Modular Exponentiation for high numbers in C++

Question

So I've been working recently on an implementation of the Miller-Rabin primality test. I am limiting it to a scope of all 32-bit numbers, because this is a just-for-fun project that I am doing to familiarize myself with c++, and I don't want to have to work with anything 64-bits for awhile. An added bonus is that the algorithm is deterministic for all 32-bit numbers, so I can significantly increase efficiency because I know exactly what witnesses to test for.

So for low numbers, the algorithm works exceptionally well. However, part of the process relies upon modular exponentiation, that is (num ^ pow) % mod. so, for example,

3 ^ 2 % 5 = 
9 % 5 = 
4

here is the code I have been using for this modular exponentiation:

unsigned mod_pow(unsigned num, unsigned pow, unsigned mod)
{
    unsigned test;
    for(test = 1; pow; pow >>= 1)
    {
        if (pow & 1)
            test = (test * num) % mod;
        num = (num * num) % mod;
    }

    return test;

}

As you might have already guessed, problems arise when the arguments are all exceptionally large numbers. For example, if I want to test the number 673109 for primality, I will at one point have to find:

(2 ^ 168277) % 673109

now 2 ^ 168277 is an exceptionally large number, and somewhere in the process it overflows test, which results in an incorrect evaluation.

on the reverse side, arguments such as

4000111222 ^ 3 % 1608

also evaluate incorrectly, for much the same reason.

Does anyone have suggestions for modular exponentiation in a way that can prevent this overflow and/or manipulate it to produce the correct result? (the way I see it, overflow is just another form of modulo, that is num % (UINT_MAX+1))

Answer 1

Exponentiation by squaring still "works" for modulo exponentiation. Your problem isn't that 2 ^ 168277 is an exceptionally large number, it's that one of your intermediate results is a fairly large number (bigger than 2^32), because 673109 is bigger than 2^16.

So I think the following will do. It's possible I've missed a detail, but the basic idea works, and this is how "real" crypto code might do large mod-exponentiation (although not with 32 and 64 bit numbers, rather with bignums that never have to get bigger than 2 * log (modulus)):

• Start with exponentiation by squaring, as you have.
• Perform the actual squaring in a 64-bit unsigned integer.
• Reduce modulo 673109 at each step to get back within the 32-bit range, as you do.

Obviously that's a bit awkward if your C++ implementation doesn't have a 64 bit integer, although you can always fake one.

There's an example on slide 22 here: http://www.cs.princeton.edu/courses/archive/spr05/cos126/lectures/22.pdf, although it uses very small numbers (less than 2^16), so it may not illustrate anything you don't already know.

Your other example, 4000111222 ^ 3 % 1608 would work in your current code if you just reduce 4000111222 modulo 1608 before you start. 1608 is small enough that you can safely multiply any two mod-1608 numbers in a 32 bit int.