Dev Notes

Software Development Resources by David Egan.

Fast Modular Exponentiation


Arithmetic, Modular
David Egan

Modular exponentiation is used in public key cryptography.

It involves computing b to the power e (mod m):

cbe (mod m)

You could brute-force this problem by multiplying b by itself e - 1 times and taking the answer mod m, but it is important to have fast (efficient) algorithms for this process to have any practical application.

In cryptography, the numbers involved are usually very large. Without an efficient algorithm, the process would take too long.

This article is educational - it is a summary of what I have learned about the process of modular exponentiation, with a few code implementations of a possible algorithm rather than a presentation of the most efficient methods.

Naive Exponentiation

THe most basic approach would be to multiply b by itself e - 1 times, and performing floor division on the result to obtain the result modulo m (which will be the remainder or residue of floor division of the result by the modulus). There are a few problems with this approach:

  1. If b and e are large numbers, be will be enormous - which causes problems representing the resultant data as a native type in many languages/systems.
  2. If b and e are large, a lot of multiplications are required. The universe might not last long enough for us to carry out the computation…

Slightly Less Naive Approach

For multiplication (mod m) congruence is maintained. In other words:

if a ≡ x (mod m) then a ∙ k ≡ x (mod m)

It follows that if we’re just concerned with congruences (i.e. the residue mod m), multiplying the congruences provides the same result as multiplying the factors and then taking the result modulo m:

If a ∙ b ≡ x (mod m), then a (mod m) ∙ a (mod m) ≡ x (mod m)

In terms of an exponentiation algorithm, multiplying the result modulo m at each step leads to much smaller numbers which spares computational resources.

Slightly Better Algorithm

For cbe (mod m)

Start with 1, multiply by b, take the result mod(m), repeat e times.

In other words:

  1. Start with c ← 1
  2. Repeat e times: ccb mod m

Exponent is a Power of Two

If cbe (mod m) and

e = 2k

We can compute c using the “squares” method - this allows for fast computation of large positive integer powers of a number.

From rules of indices:

(be)f = bef

For example, this allows a⁸, can be represented as ((a²)²)².

If you calculate a⁸ naively:

a⁸ = a ∙ a ∙ a ∙ a ∙ a ∙ a ∙ a ∙ a

…7 multiplications are required (the exponent - 1).

Alternatively, computing a⁸ as ((a²)²)² requires three multiplications:

a ∙ a = s₁ s₁ ∙ s₁ = s₂, where s₂ is equivalent to (a²)² s₂ ∙ s₂ = s₃, where s₃ is equivalent to ((a²)²)²

In this way, aⁿ requires no more than 2 log₂(e) multiplications, where e is the exponent.

So long as our exponent is a power of 2, and we multiply congruences at each stage, we have an efficient algorithm that can be converted to code:

Example Code: Exponent is a Power of 2

#!/usr/bin/env python3

def exponent_power_of_2(a, e_power2, mod):
    for i in range (0, e_power2):
        a = (a * a) % mod
    return a

Exponent is Not Necessarily a Power of Two

The previous algorithm gets us on the right track, but it has a major limitation - it only works if the exponent is a power of two.

We can however take the principle and generalise it so that it works for any number.

The idea is to take any number represented as the summation of powers of two - which as luck would have it, is exactly the way that modern computers represent numbers - and to create a running total of the required squares.

Example:

a11 = a8 ∙ a2 ∙ a1

…Notice that 8, 2 and one are powers of 2 (3, 1 and 0 respectively).

We can get the answer by working through each power of two up to the maximum possible given the size of e, squaring the base at each stage and only adding the squares to the result if the given power of two is a factor in the exponent.

Algorithm

Given a base b, an exponent e and modulo m, compute be (mod m):

  1. Create an integer (or long) variable called result and set this result equal to 1.
  2. Check the least significant bit (2⁰) of the exponent e. If it is 1, set result equal to base.
  3. Check each bit in the exponent by iteratively bitshifting and masking against 1 - this checks each position in order, starting from the second-least-significant bit (we have already considered the least most significant bit in stage 2.
  4. Start a loop
  5. At each iteration, set base equal to the value of the previous base squared, modulo m
  6. At each stage, if the LSB of e is set, set result equal to the product of the previous result and the current base (which is the previous base squared, as described in stage 3), all modulo m
  7. When the value of e is NULL, end the loop
  8. The value of result is the product of all b to the power of two for all powers of two that constitute the exponent.

In summary, set the result to be 1. Starting from the least significant bit of the exponent, iteratively check each bit - which denotes that the particular power of two is a component of the exponent. Square the base modulo m at each stage. If the exponent bit is set, multiply the base with the current result, modulo m. The final result is the base raised to the exponent mod m - the product of a set of base raised to exponents that constitute the original exponent broken into powers of two.

Example Code

Examples applying the above algorithm in C and Python are shown below:

Example: C

// See: https://github.com/DavidCWebs/input-integer-C for integer inputfunction
#include <stdio.h>
#include "integer-input.h"

int fastExp(int b, int e, int m)
{
	int result = 1;
	if (1 & e)
		result = b;
	while (1) {
		if (!e) break;
		e >>= 1;
		b = (b * b) % m;
		if (e & 1)
			result = (result * b) % m;
	}
	return result;
}

int main()
{
	puts("Compute b to the power e modulo m");
	int b = 0, e = 0, m = 0;
	puts("Enter b:");
	getIntegerFromStdin(&b);
	puts("Enter e:");
	getIntegerFromStdin(&e);
	puts("Enter m:");
	getIntegerFromStdin(&m);
	printf("%d to the power %d ≡ %d (mod%d)\n", b, e, fastExp(b, e, m), m);
	return 0;
}

Example: Python

#!/usr/bin/env python3
def fast_exp(b, e, m):
    r = 1
    if 1 & e:
        r = b
    while e:
        e >>= 1
        b = (b * b) % m
        if e & 1: r = (r * b) % m
    return r

def main():
    b = int(input("Enter b:"))
    e = int(input("Enter e:"))
    m = int(input("Enter m:"))
    r = fast_exp(b, e, m)
    print("{} ^ {} ≡ {} (mod {})".format(b, e, r, m))

if __name__ == '__main__':
    main()

References


comments powered by Disqus