I have a function ProdBigMod which computes the product of two double precise numbers $x_1$, $x_2$ (both of which are less than $2^53 - 1$) and subsequently finds the remainder $mod p$ (where $p$ is prime). I am not at liberty to use external libraries and must use double data type.

The main challenge comes into play when the product of $x_1$ and $x_2$ exceeds $2^53 - 1$.

As @TobySpeights points out, 53 is important because it is the number of mantissa bits (also called significand hence the name of the constant Significand53) for double data types (see Double-precision). We can eliminate some of these issues by first ensuring that both $x$s are less than $p$ (this is achieved by immediately applying std::fmod). In fact, when p < sqrt(2^53 - 1), we know the product prodX = x1 * x2 < 2^53 - 1. To deal with the case when p > 2^53 - 1 and prodX > 2^53 - 1, we can take advantage of some properties of modular arithmetic.

Namely:

$$(x_1 cdot x_2) pmod p = x_1 cdot (x_12 + x_22 + ... + x_n2) pmod p$$

Where $(x_12 + x_22 + ... + x_n2) = x_2$.

This gives:

$$big( x_1 cdot x_12 pmod p big) spacespace + spacespace big(x_1 cdot x_22 pmod pbig) spacespace + ... + spacespace big(x_1 cdot x_n2 pmod pbig)$$

Now, given that each of the $x_i2$ (chunkSize below) are the same except for possibly the final one (i.e. $x_n2$), we can compute the following once:

$$x_1 cdot x_12 pmod p spacespace = spacespace x_1 cdot chunkSize pmod p spacespace = spacespace chunkMod$$

To find the final $x_n2$ we first have to determine how many $x_i2$s will go into $x_2$:

numChunkMods = floor(x2 / chunkSize)

Now we can easily obtain xn2:

xn2 = (x2 - chunkSize * numChunkMods) ==>> part2 = x1 * xn2 (mod p)

Putting this altogether, Eq. (1) can be reduced to:

x1 * x2 (mod p) = (chunkMod * numChunkMods) + part2 (mod p)

This is a good start, but doesn't get us completely out of the woods, since we don't know for sure that the product numChunkMods * chunkMod < 2^53 - 1. To get around this, we continue the process above by setting x1 = chunkMod and x2 = numChunkMods, until the product is less than 2^53 - 1.

#include <iostream>
#include <cmath>

const double Significand53 = 9007199254740991.0;
const double SqrtSig53 = std::floor(std::sqrt(Significand53));

double PositiveMod(double x, double m) 
 if (x < 0)
 x = x + ceil(std::abs(x) / m) * m;
 else if (x > m)
 x = std::fmod(x, m);
 return x;


double ProdBigMod(double x1, double x2, double p) 
 double result = 0, prodX;

 x1 = PositiveMod(x1, p);
 x2 = PositiveMod(x2, p);
 prodX = x1 * x2;

 if (prodX < p) 
 result = prodX;
 else if (p < SqrtSig53 

Here is an example of how to call the function:

int main() 
 double test = ProdBigMod(914806066069, 497967734853, 732164213243);
 std::cout << std::fixed;
 std::cout << test << "n";
 return 0;


Output: 85635829849.000000

I have compared the above on random samples of numbers in the range 10^12 with a gmp analog and it seems to give correct results. However, I'm not exactly sure if my logic or implementation is bullet-proof. I'm also wondering if it could be more efficient.

Here is an online "Big Number Calculator" that you can test the results against: https://defuse.ca/big-number-calculator.htm

For example, input (914806066069 * 497967734853) % 732164213243 into the expression field and click "Calculate".

Integer-ness

It should be made more clear that ProdBigMod() is to only work with whole number values and not the full range of double including values with fractional parts, NaN and infinites.

I would expect code to detect non-whole number values and exit/complain as needed - perhaps return NaN.

Off-by-one:

"... when p < sqrt(2^53 - 1) ..." should be "... when p <= sqrt(2^53) ...".

This allows for a slightly larger p.

Recall that something mod p will be at most, p - 1.

This is reflected in the incorrect code: x > m --> x >= m

double PositiveMod(double x, double m) 
 if (x < 0)
 ...
 // else if (x > m)
 else if (x >= m)
 x = std::fmod(x, m);
 return x;

and

// } else if (p < SqrtSig53 || ...
} else if (p <= SqrtSig53 || ...

Portability

Use of Significand53 and other code relies on double as a binary64. A simple test would prevent a number of errant compilations, even if not increase portability.

#if DBL_MANT_DIG != 53
 #error TBD code
#endif

Alternative precision test.

Code does prodX = x1 * x2; and various tests when code could test FE_INEXACT after the multiplication.

feclearexcept(FE_INEXACT);
prodX = x1 * x2
if (fetestexcept(FE_INEXACT)) Handle_inexact_product();

Questionable code

I am not confident x = x + ceil(std::abs(x) / m) * m; work as expected for all x, m, especially when the math quotient is just over a whole number, but rounds down before ceil(). Clearer alternative:

double PositiveMod(double x, double m) 
 double y = std::fmod(x, m);
 if (y < 0.0) 
 y = (m < 0) ? y - m : y + m;
 
 return y;

Code alternate:

Use 64-bit or better integer math. See mulmodmax() in Modular exponentiation without range restriction. Good for at least 19 decimal digit integers vs. 12 here.

fmod() correctness.

Note: The specification of std::fmod() does not require the result to be the best possible answer, yet with IEEE Standard 754 adherence, it is. A good library would implement an exact result. Ref:
Is fmod() exact when y is an integer?
.

Conclusion

In all, using FP math to solve an integer problem has various unexpected corner concerns and so is hard to insure code is right for all x1,x2,p.