# Accurate modular arithmetic with double precision

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".

• Do you have constraints about only using 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. Feb 4, 2018 at 21:34
• @MarcGlisse, for this particular case, I'm constrained to double. Besides, even without this constraint, I'm interested in improvements, as insights to this problem can be applied to other situations dealing precision. Feb 4, 2018 at 21:58
• Well, using one multiplication and one fma, you can rewrite a*b to c+d, so if you have already written AddBigMod... Feb 4, 2018 at 22:08
• @MarcGlisse Sadly, some fma() are not that precise: Is my fma() broken?. IMO, that case is non-compliant per the intent of the function. Feb 10, 2018 at 20:59

# 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.

• Great answer. How about 'at most, 1 less than p' -> 'at most $p - 1$'? it caught me off guard for a second. Feb 10, 2018 at 21:01
• @Coal_ OK - answer amended. Feb 10, 2018 at 21:03
• mulmodmax is a really nice function, however it is a bit esoteric. I have spent a good bit of time figuring out how it works and understand it, however I don't think it is immediately obvious to the reader. An explanation would be nice if you have time. Thanks for the thorough schooling! Feb 12, 2018 at 2:48
• @JosephWood Perhsp Avoiding overflow working modulo p would help. Feb 12, 2018 at 3:05
• I think your off-by-one is itself off by one: p <= sqrt(2^53 - 1) or p < sqrt(2^53) (if I'm right that 2^53 is the lowest value that can't be exactly represented). That said, we know sqrt(2^53) is not an integer, so I'm splitting hairs here. Feb 12, 2018 at 11:05

Looks like there is another implementation using std::fma, which looks pretty simple. It is mentioned at Modular SIMD arithmetic in Mathemagix (page 18).

// https://arxiv.org/ftp/arxiv/papers/1407/1407.3383.pdf#page=18
double ProdBigMod(double a, double b, double m) {
// a, b, m - are integers within [0; 2**53-1], a < m, b < m
double p = a * b;
double r1 = std::fma(a, b, -p);
double q = std::rint(p / m);
double r2 = -std::fma(q, m, -p);
if (r1 < 0) {
r1 += m;
}
if (r2 < 0) {
r2 += m;
}
double r = r1 - m + r2;
if (r < 0) {
r += m;
}
return r;
}


Would be nice to test efficiency against the original code.

UPDATE: the code was fixed to work within the range correctly

• Welcome to Code Review! You have presented an alternative solution, but haven't reviewed the code. Please edit to show what aspects of the question code prompted you to write this version, and in what ways it's an improvement over the original. It may be worth (re-)reading How to Answer. May 17, 2021 at 16:43