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.

Solution 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).

Solution 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;
}

Solution 3:

Here is a double-and-add implementation of a multiplication a * b % m, without overflows, whatever the size of a, b and m.

(Note that the res -= m and temp_b -= m lines rely on 64-bit unsigned integer overflow in order to give the expected results. This should be fine since unsigned integer overflow is well-defined in C and C++. For this reason it's important to use unsigned integer types.)

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t res = 0;
    uint64_t temp_b;

    /* Only needed if b may be >= m */
    if (b >= m) {
        if (m > UINT64_MAX / 2u)
            b -= m;
        else
            b %= m;
    }

    while (a != 0) {
        if (a & 1) {
            /* Add b to res, modulo m, without overflow */
            if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
                res -= m;
            res += b;
        }
        a >>= 1;

        /* Double b, modulo m */
        temp_b = b;
        if (b >= m - b)       /* Equiv to if (2 * b >= m), without overflow */
            temp_b -= m;
        b += temp_b;
    }
    return res;
}

This is my modification of another answer to another similar question.