PrimeTest.html

Primality testing

> restart;

Since we need large primes for the public key cryptosystems, we need methods to find if a number is prime without finding factors,

A simple procedure for primality testing

In class we looked at some probabilistic tests of primality. We want to see how effective the tests are. We will simply eyeball the results to see how well the tests work.

For the first test, we test for primality by using Fermat's little theorem. We test if n is prime by taking a random number a < n and seeing if a^(n-1) = 1 mod n. We have been told that if n is not prime we have at least a 50% chance of showing it by choosing a single number.

It is worthwhile checking out the effectiveness of this test with a specific n. The first time through we will look at the case of n = 77.

We first define a vector with all the numbers from 1 through n-1.

> alphas := linalg[vector](76,i->i);

Then we see what happens when we raise these numbers to the n-1 power mod n.

> betas := map(i-> [i,Power(i,76) mod 77], alphas);

Thus the test would tell us 77 is not prime if the chosen a was anything except

1, 34, 43, and 76. Thus the test indicates that n is not prime for 72 of 76 possible values of a. When we consider that the test is always inconclusive if a is 1 or -1, we see that the test indicates 77 not prime for 72 of 74 possible values for a between 2 and n-2. That is considerably more effective than 50%.

Test some other numbers:

The following block of code is set up to test other numbers using the first test.

> n := 131;
alpha1 := linalg[vector](n-1,i->i);
beta1 := map(i -> [i, Power(i, n-1) mod n], alpha1);

Exercise:

1) Pick a nonprime number n between 100 and 150. Find out how many of the numbers between 2 and n-2 will show that n is composite.

>

2) Repeat the previous exercise for another number between 100 and 150.

>

Procedure for testing higher numbers:

Obviously, we can't use this same procedure for large numbers because we will run out of memory trying to store the arrays. A simple alternative is to test with a running from 2 to 100. The procedure is given below.

A simple prime tester

>    primetester := proc(n)
  local i, temp:
  for i from 2 to 100 do
    temp := Power(i,n-1) mod n:
    if temp <> 1 then
      print(n,` is not prime`):
      print(cat(i,` ^ `,n-1,` = `,temp)):
      break:
    fi:od:
  if i > 100 then print(cat(n,` is probably prime`)) fi:
  if i > 100 then n else 0 fi:
end:

> primetester(1122334455667789);

> primetester(131);

Finding primes

Now that we can test if a number is prime, we can try to find a prime by testing numbers in order until we find a prime

>    findprime := proc(n)
  local temp, temp2, a, b, c:
  temp := rand(n)():
  print("The random number is ", temp):
  a := 0:
  c:= 0:
  while a = 0 do:
    b := primetester(temp):
    temp := temp +1:
    c := c+1:
    if b <> 0 then a := 1:
  fi:od:
  print(`After `,c-1,` unsuccessful tries, the prime number is `, b):
  b:
end:

Let us use this procedure to find an 8 digit prime. We use the Maple command isprime to test our answer.

> findprime(10^8);
isprime(%);

Exercise:

3) Find two primes with 7 digits.

>

As the size of the prime we are looking at grows, the number of tries needed also grows. Thus we want to modify the findprime so that it only gives the answer.

>    primetester2 := proc(n)
  local i, temp:
  for i from 2 to 100 do
    temp := Power(i,n-1) mod n:
    if temp <> 1 then
      break:
    fi:od:
  if i > 100 then print(cat(n,` is probably prime`)) fi:
  if i > 100 then n else 0 fi:
end:
findprime2 := proc(n)
  local temp, temp2, a, b, c:
  temp := rand(n)():
  print("The random number is ", temp):
  a := 0:
  c:= 0:
  while a = 0 do:
    b := primetester2(temp):
    temp := temp +1:
    c := c+1:
    if b <> 0 then a := 1:
  fi:od:
  print(`After `,c-1,` unsuccessful tries, the prime number is `, b):
  b:
end:

> findprime2(10^20);

Exercise:

4) Find primes with 10, 20, 30, 40, and 50 digits.

>

Oops - Carmichael numbers

Let us try our test on 3828001 = 101*151*251 which is clearly not prime.

> primetester(3828001): isprime(3828001);ifactor(3828001);

The number 3828001 was specially chosen so that the order of the multiplicative group divides n-1 even though n is not prime. Such numbers are called Carmichael numbers. We would still spot the that n is not prime if we tested a number that is not relatively prime to n. Unfortunately, n was chosen so that its smallest prime factor is 101. Thus we need to use the more sensative Miller test.

The Miller-Rabin Test of primality.

We also want to see what happens with the more sophisticed test, the Miller-Rabin test. For that test we factor n-1 as for some odd number k. We then look at a raised to the powers k, 2*k, 2^2*k, ..., 2^t*k = n-1. If n is prime, then a^(n-1) = 1 mod n. Additionally, if a^(2j) = 1 mod n, then a^j mod n is either 1 or -1.

We illoustrate this test with n=77. We notice that . Thus we want to look at raising a to the powers 19, 38, and 76.

For the Miller-Rabin test with a to fail to show n is composite,we need the last power to be 1. We also need for the previous power, if it exists, to be 1 or -1. Recall that 76 = -1 mod 77.

>    gammas := map(i ->
       [i, Power(i,19) mod 77, Power(i,38) mod 77, Power(i,76) mod 77],
            alphas);

Thus, the test shows 77 is not a prime unless we choose a to be 1 or -1. In effect the test shows that n is composite if we choose any a between 2 and n-2.

To do the Miller-Rabin test we first need to express n-1 as an odd number times a power of 2. We do that with the procedure oddfactor. The procedure WitnessTest below performs the Miller-Rabin test and look for numbers from 2 to 100 (or n-2 if it is lower) that fails to witness to the fact that n is composite.

>    WitnessTest := proc(n)
local n1, nprimelist, primelist, inconlist, i, numprime, numnprime,
   numincon, m, testlist, testlength, a, p1, p2,top:
top := min(n-1,101):
n1 := n-1:
nprimelist := []:
primelist := []:
#   First we check that a^(n-1) = 1 mod n.
#   We break the numbers 2 to 100 into two groups,
#    those that show n is not prime, and
#    those that fail this test, but may work with more testing.
for i from 2 to (top-1) do
   if ((Power(i,n-1) mod n) = 1) then primelist := [op(primelist), i]:
   else nprimelist := [op(nprimelist), i]:
   fi:
od:
numnprime := linalg[vectdim](nprimelist):
#   Print the results up to this stage
print(cat("Testing a^",n1," shows n is not prime for ", numnprime,
   " numbers"));
print(nprimelist);
m := n1:
#   The next step check to see if a^m = 1 or -1 when a^(2m) = 1.
while (m mod 2 = 0) do
   m := m/2:
   testlist := primelist:
   testlength := linalg[vectdim](testlist);
   nprimelist := []:
   primelist := []:
   for i from 1 to testlength do
      a := testlist[i]:
      p1 := Power(a, 2*m) mod n:
      p2 := Power(a, m) mod n:
      if (p1 = 1) then
        if ((p2 <> 1) and (p2 <> n-1))
              then nprimelist := [op(nprimelist), a]:
           else primelist := [op(primelist), a]:
        fi:
      else inconlist := [op(inconlist), a]:
      fi:
   od:
   numnprime := linalg[vectdim](nprimelist):
   print(cat("Testing a^",m," shows n is not prime for ", numnprime,
      " additional numbers"));
   print(nprimelist);
od:
numprime := linalg[vectdim](primelist):
print(cat("The test fails for ", numprime," numbers"));
print(sort(primelist));
end proc:

> WitnessTest(3828001);

>

Theortic arguments say that 3 of 4 numers will witness that n is composite. Note that in our test case only 19 of the numbers from 2 through 100 fail to witness that n is not prime.

Compare the results to a few other Carmichael numbers. With two Carmichael numbers we note that the first primality test fails to identify them as composite numbers, then note that the Miller test does identify them as composites, then we factor the numbers.

> primetester(6733693): WitnessTest(6733693): ifactor(6733693);
primetester(34657141): WitnessTest(34657141): ifactor(34657141);

In these two cases respectively 86% and 82% of the numbers from 2 to 100 are witnesses that the numbers are not prime.

Exercise

5) Pick 3 random 50 digit numbers. Test if each to see if it is prime and if it is composite, note the percentage of choices of a from 2 to 100 that will show they are composite with the Miller test.

>

It is woth noting that we have yet to find a composite n that we can't show to be composite by testing the number 2. Such a number would be called a strong probable-prime base 2. It has been shown in the literature that letting the base range over the primes less than 20 will correctly identify all composites less than 10^15.

We can clean up the code and give a better prime finding routine. In particular we eliminate the printing of the unsuccessful tries

>    primetester4 := proc(n)
local a, temp, temp2, k, m:
temp2 := n:
k:= n-1:
while (k mod 2 = 0) do k:= k/2 od:
for a from 2 to 100 do
    m := k:
    temp := Power(a,k) mod n:
    if ((temp = 1) or (temp = n-1)) then break:
    else while (m < n-1) do
       m := 2*m:
       temp := Power(temp,2) mod n:
       if (temp = n-1) then break:
         elif (temp = 1) then temp2 := 0: break:
         fi:
    od:
    if ((m = n-1) and (temp <> 1)) then temp2 := 0 fi:
    fi:
    if (temp2 = 0) then break fi:
od:
temp2;
end:
findprime4 := proc(n)
local temp, a, b, c:
temp := rand(n)():
temp := temp + (1 - irem(temp,2)):
print(`the random number is `, temp):
a := 0:
c:= 0:
while a = 0 do:
      b := primetester4(temp):
      temp := temp +2:
      c := c+1:
      if b <> 0 then a := 1: fi:
   od:
print(`After `,c-1,` unsuccessful tries, the prime number is `, b):
b:
end:

> findprime4(10^100); isprime(%);

Exercise:

6) Use the findprime4 procedure to find 60, 70, and 80 digit primes.

>

Maple's method of finding primes

It is worthwhile seeing how Maple finds primes with the nextprime command.

> nextprime(rand(10^10)());

We can look at the code used with the showstat command.

> showstat(nextprime);

nextprime := proc(n)

local i, t1;

   1   if type(n,integer) then

   2     if n < 2 then

   3       2

         else

   4       if irem(n,2) = 0 then

   5         t1 := n+1

           else

   6         t1 := n+2

           end if;

   7       for i from t1 by 2 while not isprime(i) do

   8         NULL

           end do;

   9       i

         end if

       elif type(n,numeric) then

  10     error "argument must be integer"

       else

  11     'nextprime(n)'

       end if

end proc

This is essentially the same algorithm we used. From the starting point, check odd numbers going up untill a prime is found. We then want to see how Maple checks is=if a number is prime.

> showstat(isprime);

isprime := proc(n)

local btor, nr, p, r;

   1   if not type(n,'integer') then

   2     if type(n,'numeric') then

   3       error "argument must be an integer"

         else

   4       return 'isprime(n)'

         end if

       end if;

   5   if n < 2 then

   6     false

       elif member(n,`isprime/w`) then

   7     true

       elif igcd(2305567963945518424753102147331756070,n) <> 1 then

   8     false

       elif n < 10201 then

   9     true

       elif igcd(8496969489233418110532339909187349965926062586648932736611545426342203893270769390909069477309509137509786917118668028861499333825097682386722983737962963066757674131126736578936440788157186969893730633113066478620448624949257324022627395437363639038752608166758661255956834630697220447512298848222228550062683786342519960225996301315945644470064720696621750477244528915927867113,n) <> 1 then

  10     false

       elif n < 1018081 then

  11     true

       else

  12     nr := igcd(408410100000,n-1);

  13     nr := igcd(nr^5,n-1);

  14     r := iquo(n-1,nr);

  15     btor := modp(('power')(2,r),n);

  16     if `isprime/cyclotest`(n,btor,2,r) = false or irem(nr,3) = 0 and `isprime/cyclotest`(n,btor,3,r) = false or irem(nr,5) = 0 and `isprime/cyclotest`(n,btor,5,r) = false or irem(nr,7) = 0 and `isprime/cyclotest`(n,btor,7,r) = false then

  17       return false

         end if;

  18     if isqrt(n)^2 = n then

  19       return false

         end if;

  20     for p from 3 while numtheory[jacobi](p^2-4,n) <> -1 do

  21       NULL

         end do;

  22     evalb(`isprime/TraceModQF`(p,n+1,n) = [2, p])

       end if

end proc

The isprime routine first check if the test candidate n is in the list of primes below 100. If not it checks to see if any of those primes divide n by computing the gcd of n and the product of primes less than 100, thus testing if n is a prime less than (101)^2. If not, it computes the gcd of n and the product of primes less than 1000 to test for primes less than (1009)^2. Failing to find that n has any small factors, Maple then turns to techniques beyond the scope of this class

>