Modular Exponentiation and Primality Testing

Math 242 Modern Computational Math, spring 2021

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 the difference 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]

1234232890895376789452390234^23489329084230982342436556

Modular Exponentiation

We will need to compute the value of (b^e) (mod m) for large values of b, e, and m. We could try to simply compute (b^e) and then reduce mod m. However, for large values of b, e, and m, this will be inefficient or impossible.

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^((8+ 2+ 1)))=(b^8)(b^2)(b^1). 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

Complete the module below to efficiently compute (b^e) (mod m).

modPow[base_,exponent_,modulus_]:=Module[{b=base,e=exponent,result},
result=1;
b=Mod[b,modulus];
While[e > 0,
If[Mod[e,2]==1,result = Mod[(result*b),modulus];
];
e=Floor[e/2];
b=Mod[b*b,modulus];
];
result   (* return the result *)
]

How efficient is modPow[b, e, m] compared to just computing Mod[(b^e), m]? Your function modPow should work efficiently even when b and e have thirty or more digits each!

Fermat’s Little Theorem and Primality Testing -- we’ll get to this on Friday

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]) less than n and compute

((Subscript[a, i])^(n-1)) (mod n)    for each    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:

fermatPrime[n_,m_]:=Module[{},
(* 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?

Try it with some known primes:

fermatPrime[104393,10]

fermatPrime[282174488599599500573849980909,10]

Also try it with composites:

fermatPrime[282174488598,10]

fermatPrime[1853020188851841,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.