Fast Modular Exponentiation

October 6, 2014

The Problem

… is to compute be (mod m), where b, e, and m are positive integers.

e.g., for b = 2, e = 109 and m = 109 + 7,

be (mod m) = 140625001

A Naive Solution

… would be one that first computes be and then computes its remainder mod m.

Speaking in terms of the basic arithmetic operations, this method would require O(e) multiplications, and a lot of space to store the computed value of be.

Consider the example given above where b is 2 and e is 109. Storing the computed value of be in this case will require (109 + 1) bits! Clearly, such a solution will not cut it with large values of b, e and m.

Improving the space bound

The costly O(log2 be) bound on space can be improved significantly if, instead of computing the remainder at the end, we compute remainders after each multiplication, using the property

{ p × q } (mod r) = [ p × { q (mod r) } ] (mod r)

Plugging in values for b, e and m from the example above, we begin multiplying as follows:

b (mod m) = 2 (mod m) = 2

b2 (mod m) = [ b × { b (mod m) } ] (mod m) = [ 2 × 2 ] (mod m) = 4

b3 (mod m) = [ b × { b2 (mod m) } ] (mod m) = [ 2 × 4 ] (mod m) = 8

…and so on, till

be (mod m) = [ b × { be - 1 (mod m) } ] (mod m)

This solution runs in O(e) time, using O(log2 (b × m)) space. The space improvement is significant for large values of e.

Here’s a sketch for a recursive function based on the algorithm above, taken from a Wikipedia article:

function modular_pow(b, e, m)  
    c := 1
    for e_prime = 1 to e
        c := (c * b) mod m
    return c

Improving the time bound

Now that we’ve sufficiently improved the worst-case space bound, let’s do something about the O(e) time bound, which can be troublesome for large values of e. Consider the binary representation for the exponent e:

e = ( a0 × 20 ) + ( a1 × 21 ) + … + ( an × 2n ), where n = log2 e

Thus,

be (mod m) = b( a0 × 20 ) + ( a1 × 21 ) + … + ( an × 2n ) (mod m)

be (mod m) = { b( a0 × 20 ) (mod m) } × { b( a1 × 21 ) (mod m) } × … × { b( an × 2n ) (mod m) }

be (mod m) = { (b20 )a0 (mod m) } × { (b21 )a1 (mod m) } × … × { (b2n )an (mod m) }

At this point, two observations are key:

  • When ai = 0, the term { (b2i)ai (mod m) } reduces to 1
  • When ai = 1, the same term becomes equal to b2i (mod m), which can be computed iteratively as { b2 (mod m) × b2(i - 1) (mod m) } (mod m)

This leaves us with a fast and space-efficient algorithm running in only O(log2 e) time and requiring only O(log2 (b × m)) space. Here’s a sketch of the complete algorithm:

function modular_pow(b, e, m)  
    result := 1
    b := b mod m
    while e > 0
        if (e mod 2 == 1):
           result := (result * b) mod m
        e := e >> 1
        b := (b * b) mod m
    return result
Tagged with