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) \space\space + \space\space \big(x_1 \cdot x_{22} \pmod p\big) \space\space + ... + \space\space \big(x_1 \cdot x_{n2} \pmod p\big)$$
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 \space\space = \space\space x_1 \cdot chunkSize \pmod p \space\space = \space\space 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 || prodX < Significand53) {
result = std::fmod(prodX, p);
} else {
double numChunkMods, part1 = Significand53;
double chunkSize, chunkMod, part2;
while (part1 >= Significand53) {
chunkSize = std::floor(Significand53 / x1); // Ensures chunkMod < 2^53 - 1
chunkMod = std::fmod(x1 * chunkSize, p);
numChunkMods = std::floor(x2 / chunkSize);
part2 = std::fmod((x2 - chunkSize * numChunkMods) * x1, p);
part1 = numChunkMods * chunkMod;
x1 = chunkMod;
x2 = numChunkMods;
result = std::fmod(result + part2, p);
}
result = std::fmod(part1 + result, p);
}
return result;
}
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".
double
? If you could convert to uint128_t... This reminds me a bit of double-double libraries, where before FMA (are you allowed to use fma?) multiplication would split each operand into a sum of short-mantissa numbers. \$\endgroup\$double
. Besides, even without this constraint, I'm interested in improvements, as insights to this problem can be applied to other situations dealing precision. \$\endgroup\$fma()
are not that precise: Is my fma() broken?. IMO, that case is non-compliant per the intent of the function. \$\endgroup\$