Modular exponentials

Modular exponentials are used in various different algorithms, mostly in cryptography. Let’s dive right in with an example: given 2 large positive integers a and p, we want to compute the following with 32-bit integers (uint32_t).

$$a^p\mod p$$

Let’s pretend that 53 and 91 are large numbers - we then have

$$53^{91}\mod 91$$

$53^{91}$ is an astronomical number, which wouldn’t fit in less than 500 or so bits, so we need to come up with an indirect way of computing its value modulo 91. The first algorithm that might come to mind is the following one, where we iteratively multiply an accumulator and take its value modulo 91 each time.

// assume mod > 1, ignore 0^0
uint32_t mod_pow(uint32_t base, uint32_t exp, uint32_t mod) {
    base %= mod;
    if (base == 0) return 0;

    uint32_t accumulator = 1;
    for (uint32_t i = 0; i < exp; ++i) {
        accumulator = (accumulator * base) % mod;
    }

    return accumulator;
}

One crucial detail to keep in mind here is that accumulator * base might not fit into a 32-bit integer. Since both accumulator and base must be less than mod - 1, all we need to do is verify at the beginning of the function that (mod - 1) * (mod - 1) fits into a uint32_t.

Can’t we do better than this naive algorithm though? The answer (luckily) is yes: you can find some pseudocode on this page, but I’ll try to explain it below.

uint32_t mod_pow(uint32_t base, uint32_t exp, uint32_t mod) {

    base %= mod;
    if (base == 0) return 0;

    uint32_t accumulator = 1;
    while (true) {
        if (exp & 1) accumulator = (accumulator * base) % mod;
        exp >>= 1;
        if (exp == 0) break;
        base = (base * base) % mod;
    }

    return accumulator;
}

The trick to understand this algorithm is to write the exponent (91) in binary form (1011011):

$$91 = 2^6 + 2^4 + 2^3 + 2^1 + 2^0$$ $$\Rightarrow 53^{91} = 53^{2^6}\cdot 53^{2^4} \cdot 53^{2^3} \cdot 53^{2^1} \cdot 53$$

We can then build up the result by multiplying the accumulator with the value corresponding to each 1 bit, while skipping all the 0s. Let’s take a closer look at the while loop:

• if (esp & 1): if the exponent is odd, we multiply the accumulator by the base. Since, in general, a * b % p == (a % p) * (b % p), we can apply the modulo operator at any time.
• exp >>= 1: move onto the next bit of the exponent.
• base = (base * base) % mod: the next bit corresponds to the square of the current base. For example, the fifth bit corresponds to:

$$53 ^ {2^4} = 53^{2\cdot 2^3} = (53^{2^3})^2$$

And that’s all, folks! The first algorithm has runtime O(exp), whereas the second one is O(log(exp)), so the winner is pretty clear.