Computing a^x mod n, where x is not a power of 2, is only slightly harder. Binary notation expresses x as a sum of powers of 2: 25 is 11001 in binary, so 25 = 2⁴ + 2³ + 2⁰. So

a²⁵ mod n = (a*a²⁴) mod n = (a*a⁸*a¹⁶) mod n

= (a*((a²)²)²*(((a²)²)²)²) mod n = ((((a²*a)²)²)²*a) mod n

With judicious storing of intermediate results, you only need six multiplications:

(((((((a² mod n)*a) mod n)² mod n)² mod n)² mod n)*a) mod n

This is called addition chaining [863], or the binary square and multiply method. It uses a simple and obvious addition chain based on the binary representation. In C, it looks like:

    unsigned long qe2(unsigned long x, unsigned long y, unsigned long n) {
    	unsigned long s,t,u;
    	int i;
    
    	s = 1; t = x; u = y;
    	while(u) {
    		if(u&1) s = (s* t)%n;
    		u>>=1;
    		t = (t* t)%n;
    	}
    	return(s);
    }

Another, recursive, algorithm is:

    unsigned long fast_exp(unsigned long x, unsigned long y, unsigned long N) {
    	unsigned long tmp;
         if(y==1) return(x % N);
         if ((y&1)==0) {
    	       tmp = fast_exp(x,y/2,N);
    	       return ((tmp* tmp)%N);
         }
         else {
    	       tmp = fast_exp(x,(y-1)/2,N);
    	       tmp = (tmp* tmp)%N;
    	       tmp = (tmp* x)%N;
    	       return (tmp);
         }
    }

This technique reduces the operation to, on the average, 1.5*k operations, if k is the length of the number x in bits. Finding the calculation with the fewest operations is a hard problem (it has been proven that the sequence must contain at least k- 1 operations), but it is not too hard to get the number of operations down to 1.1*k or better, as k grows.

An efficient way to do modular reductions many times using the same n is Montgomery’s method [1111]. Another method is called Barrett’s algorithm [87]. The software performance of these two algorithms and the algorithm previously discussed is in [210]: The algorithm I’ve discussed is the best choice for singular modular reductions; Barrett’s algorithm is the best choice for small arguments; and Montgomery’s method is the best choice for general modular exponentiations. (Montgomery’s method can also take advantage of small exponents, using something called mixed arithmetic.)

The inverse of exponentiation modulo n is calculating a discrete logarithm . I’ll discuss this shortly.

Prime Numbers

A prime number is an integer greater than 1 whose only factors are 1 and itself: No other number evenly divides it. Two is a prime number. So are 73, 2521, 2365347734339, and 2⁷⁵⁶⁸³⁹ - 1. There are an infinite number of primes. Cryptography, especially public-key cryptography, uses large primes (512 bits and even larger) often.

Evangelos Kranakis wrote an excellent book on number theory, prime numbers, and their applications to cryptography [896]. Paulo Ribenboim wrote two excellent references on prime numbers in general [1307, 1308].

Greatest Common Divisor

Two numbers are relatively prime when they share no factors in common other than 1. In other words, if the greatest common divisor of a and n is equal to 1. This is written:

gcd(a,n) = 1

The numbers 15 and 28 are relatively prime, 15 and 27 are not, and 13 and 500 are. A prime number is relatively prime to all other numbers except its multiples.

One way to compute the greatest common divisor of two numbers is with Euclid’s algorithm . Euclid described the algorithm in his book, Elements, written around 300 B.C. He didn’t invent it. Historians believe the algorithm could be 200 years older. It is the oldest nontrivial algorithm that has survived to the present day, and it is still a good one. Knuth describes the algorithm and some modern modifications [863].

In C:

    /* returns gcd of x and y */
    int gcd (int x, int y)
    {
         int g;
         if (x < 0)
    	  x = -x;
         if (y < 0)
    	  y = -y;
         if (x + y == 0)
    	  ERROR;
         g = y;
         while (x > 0) {
    	  g = x;
    	  x = y % x;
    	  y = g;
         }
         return g;
    }

This algorithm can be generalized to return the gcd of an array of m numbers:

    /* returns the gcd of x1, x2...xm */
    int multiple_gcd (int m, int *x)
    {
         size_t i;
         int g;

         if (m < 1)
    	       return 0;
         g = x[0];
         for (i=1; i<m; ++i) {
    	       g = gcd(g, x[i]);
    /* optimization, since for random x[i], g==1 60% of the time: */
    	        if (g == 1)
    	             return 1;
         }
         return g;
    }

Inverses Modulo a Number

Remember inverses? The multiplicative inverse of 4 is 1/4, because 4*1/4 = 1. In the modulo world, the problem is more complicated:

4*x ≡ 1 (mod 7)

This equation is equivalent to finding an x and k such that

4x = 7k + 1

where both x and k are integers.

The general problem is finding an x such that

1 = (a*x) mod n

This is also written as

a ^-1 ≡ x (mod n)

The modular inverse problem is a lot more difficult to solve. Sometimes it has a solution, sometimes not. For example, the inverse of 5, modulo 14, is 3. On the other hand, 2 has no inverse modulo 14.

In general, a^-1 ≡ x (mod n) has a unique solution if a and n are relatively prime. If a and n are not relatively prime, then a-1 ≡ x (mod n) has no solution. If n is a prime number, then every number from 1 to n- 1 is relatively prime to n and has exactly one inverse modulo n in that range.