Modular Exponentiation and Primality Testing

Math 242 Modern Computational Math

Modular Arithmetic

The notation a≡b mod n means that a and b have the same remainder when divided by n, or equivalently, that n divides a-b. For example, 2≡23 mod 7, since 7 divides (23-2).

Mathematica can compute modular arithmetic:

Mod[23,7]

Given a modulus m (say 7), the value a mod m is in the range 0,1,...,m-1:

m=7;
Mod[#,m]&/@Range[0,27]

Modular Powers: The Simple Version

In what follows, it will be very important to compute the value of (b^e) mod m for large values of b, e, and m. Mathematica can, of course, just compute (b^e) and then reduce mod m. However, for large values of b, e, and m, this may be inefficient or impossible.

Here is a first version, written in pseudocode:

result = 1
   repeat e times:
        result = (result * b) mod m
    return result

This code accumulates the result of repeated multiplication by b. Since the modulus is computed at each iteration of the loop, the numbers won’t get too big.

Write your own Module to implement this algorithm, according to the following specification.

Module: modPow
    Arguments: base b, exponent e, modulus m
    Return: (b^e) mod m

How efficient is modPow[b, e, m], in terms of computing time, compared to just computing Mod[(b^e), m]? Check it out on different values of b, e, and m. Push these up high enough that you can determine clear patterns in the runtimes.

Repeated Squaring

We can greatly improve the efficiency using the base two expansion of the exponents. For example, if e=11=8+2+1=(2^3)+(2^1)+(2^0), then (b^e)= (b^(((2^3)+ (2^1)+ (2^0))))=(b^(2^3))(b^(2^1))(b^(2^0)). The following pseudocode exploits this idea:

result = 1
   b = b mod m
   while e > 0
        if (e mod 2 == 1):
           result = (result * b) mod m
        e = floor[e / 2]
        b = (b * b) mod m
    return result

Write your own Module to implement this algorithm, according to the following specification:

Module: modPow2
    Arguments: base b, exponent e, modulus m
    Return: (b^e) mod m

How efficient is modPow2[b, e, m], time-wise, compared to just computing Mod[(b^e), m] and modPow[b, e, m]? As before, try it out on different values of b, e, and m. This new function modPow2 should work efficiently even for b and e on the order of (10^20).

Fermat’s Little Theorem and Primality Testing

Fermat’s Little Theorem says that if p is prime and a is not divisible by p, then (a^(p-1))≡1 mod p.

This simple observation provides a way to test if an arbitrary number n is prime. To do so, pick some values (Subscript[a, 1]), (Subscript[a, 2]), ..., (Subscript[a, m]) (not divisible by n) and compute

((Subscript[a, i])^(n-1))mod n    for    i=1, 2, ..., m.

If all of these are 1, then n is probably prime.

Here is a brief demonstration. Note that 97 is prime.

Mod[#^96,97]&/@Range[2,96]

However, 93 is not prime, and the modular powers (mostly) are not equal to 1.

Mod[#^92,93]&/@Range[2,92]

Write your own Module to implement the following prime-testing algorithm:

Module: fermatPrime
Arguments: n, m
Return: Select m random values of a in the range 2, ..., n-1. If  (a^(n-1))≡1 mod n for each of these values of a, then return True. Otherwise return False.

Implement this and test it well. One way to test this is to use your best sieve algorithm to generate a large list of primes. Then apply fermatPrime to some or all of these. Does it correctly identify them as prime? Can you find any values for which fermatPrime does not work?

Known primes:

fermatPrime[104393,10]

fermatPrime[282174488599599500573849980909,10]

Try it with composites:

fermatPrime[282174488598,10]

What' s next?

For our encryption activities starting on Friday, you will need two large prime numbers that no one else knows. Use your fermatPrime algorithm to find two primes of at least 30 decimal digits. Keep these in a safe place.