Calculate a*a mod n without overflow

Question

I need to calculate a*a mod n but a is fairly large, resulting in overflow when I square it. Doing ((a % n)*(a % n)) % n doesn't work because (n-1)² can overflow. This is in C++ and I'm using int64_t

Edit:

Example value: a = 821037907258 and n = 800000000000, which overflows if you square it.

I am using DevCPP and I've already tried getting big-integer libraries working to no avail.

Edit 2:

No, there's no pattern to these numbers.

Answer 1

If you can't use a big-integer library, and you don't have a native uint128_t (or similar), you'll need to do this manually.

One option is to express a as the sum of two 32-bit quantities, i.e. a = 2³²b + c, where b contains the 32 msbs, and c contains the 32 lsbs. Squaring is then a set of four cross-multiplications; each result is guaranteed to fit into a 64-bit type. You then do the modulo operation as you recombine the individual terms (carefully taking into account the shifts needed to realign everything).

Answer 2

I know you no longer need this, and there is an alternative solution, but I want to add an alternative method to implement it. It provides two different techniques: the double and add algorithm, and the method to handle mod(a + b, n) with overflow detection.

Double and add algorithm is usually used in fields where multiplication is not possible or too costly to calculate directly (such as elliptic curves), but we could adopt it to handle it in our situation to handle overflows instead.

The following code is probably slower than the accepted solution (even when you optimize it), but if speed is not critical, you may prefer it for clarity.

unsigned addmod(unsigned x, unsigned y, unsigned n)
{
    // Precondition: x<n, y<n
    // If it will overflow, use alternative calculation
    if (x + y <= x) x = x - (n - y);
    else x = (x + y) % n;
    return x;
}

unsigned sqrmod(unsigned a, unsigned n)
{
    unsigned b;
    unsigned sum = 0;

    // Make sure original number is less than n
    a = a % n;

    // Use double and add algorithm to calculate a*a mod n
    for (b = a; b != 0; b >>= 1) {
        if (b & 1) {
            sum = addmod(sum, a, n);
        }
        a = addmod(a, a, n);
    }
    return sum;
}