Accurate modular arithmetic with double precision

Question

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

Answer 1

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.

Answer 2

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