Application Center - Maplesoft

App Preview:

Test for primes

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

Primality testing

Mike May, S. J., 1998

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 see 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);

alphas := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...
alphas := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...

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);

betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...
betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...
betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...
betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...
betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...
betas := vector([[1, 1], [2, 9], [3, 25], [4, 4], [...

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 not prime for 72 of 74 possible a between 2 and n-2.

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);

n := 131

alpha1 := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...
alpha1 := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...
alpha1 := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...
alpha1 := vector([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11...

beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...
beta1 := vector([[1, 1], [2, 1], [3, 1], [4, 1], [5...

Exercise:

1) Pick 3 nonprime numbers between 100 and 150. For each one, find out how many of the numbers between 2 and n-2 will shouw that each number is composite.

>

A 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(i,` ^ `,n-1,` = `,temp):
break:
fi:od:
if i > 100 then print(n,` is probably prime`,i) fi:
if i > 100 then n else 0 fi:
end:

> primetester(1122334455667789);

1122334455667789, ` is not prime`

2, ` ^ `, 1122334455667788, ` = `, 210589670002713

0

> primetester(131);

131, ` is probably prime`, 101

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(");

>

Warning, incomplete string;  use " to end the string

Exercise:

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

>

Oops

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

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

3828001, ` is probably prime`, 101

false

``(101)*``(151)*``(251)

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 Test of primality.

We also want to see what happens with the more sophisticed test, the Miller test. For that test I factor n-1 as k*2^t 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.

I notice that 76 = 2^2*19 . Thus I want to look at raising a to the powers 19, 38, and 76.

For a test to show prime, I need the last power to be 1. I 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);

gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...
gammas := vector([[1, 1, 1, 1], [2, 72, 25, 9], [3,...

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 not prime if we choose any a between 2 and n-2.

To do the Miller 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 primetester2 below performs the Miller test and look for number that the Miller fail to witness to the fact that n is not prime.

> oddfactor := proc(n)
local temp, count:
temp := n:
count := 0:
while (temp mod 2) = 0 do
temp := temp/2:
count := count + 1:
od:
[temp, count];
end:
primetester2 := proc(n)
local i, temp, nums, count, temp2:
temp2 := n:
nums := oddfactor(n-1):
for i from 2 to 100 do
temp := Power(i,n-1) mod n:
if temp <> 1 then
if temp2 <> 0 then
print(n,` is not prime`):
print(i,` ^ `,n-1,` = `,temp):
fi:
temp2 := 0:
break:
else
for count from 1 to nums[2] do
temp := Power(i,(n-1)/2^count) mod n:
if (temp <> 1) and ((n-temp) <> 1) then
if temp2 <> 0 then
print(n,` is not prime`):
### WARNING: note that `I` is no longer of type `^`
print(i,`^`,(n-1)/2^count,` = `,temp,` and `,
### WARNING: note that `I` is no longer of type `^`
i,`^`,(n-1)/2^(count-1),`=1`):
fi:
temp2 := 0:
break:
fi:
if temp = -1 then
print(i,` ^ `,(n-1)/2^count,` = `,temp):
break:
fi:
od:
if count > nums[2] then print(i,` is not a witness`):fi
fi:
od:
if temp2 = n then print(n,` is probably prime`,i) fi:
temp2;
end:

> primetester2(3828001);

3828001, ` is not prime`

2, proc () option builtin; 85 end proc, 239250, ` =...

5, ` is not a witness`

6, ` is not a witness`

14, ` is not a witness`

16, ` is not a witness`

25, ` is not a witness`

30, ` is not a witness`

31, ` is not a witness`

33, ` is not a witness`

36, ` is not a witness`

58, ` is not a witness`

68, ` is not a witness`

70, ` is not a witness`

77, ` is not a witness`

80, ` is not a witness`

81, ` is not a witness`

82, ` is not a witness`

84, ` is not a witness`

88, ` is not a witness`

96, ` is not a witness`

0

Note that only 18 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): primetester2(6733693): ifactor(6733693);
primetester(34657141): primetester2(34657141): ifactor(34657141);

6733693, ` is probably prime`, 101

6733693, ` is not prime`

2, proc () option builtin; 85 end proc, 3366846, ` ...

9, ` is not a witness`

12, ` is not a witness`

16, ` is not a witness`

21, ` is not a witness`

22, ` is not a witness`

25, ` is not a witness`

26, ` is not a witness`

28, ` is not a witness`

29, ` is not a witness`

31, ` is not a witness`

49, ` is not a witness`

81, ` is not a witness`

82, ` is not a witness`

97, ` is not a witness`

``(109)*``(163)*``(379)

34657141, ` is probably prime`, 101

34657141, ` is not prime`

2, proc () option builtin; 85 end proc, 17328570, `...

3, ` is not a witness`

7, ` is not a witness`

9, ` is not a witness`

16, ` is not a witness`

20, ` is not a witness`

21, ` is not a witness`

25, ` is not a witness`

27, ` is not a witness`

46, ` is not a witness`

48, ` is not a witness`

49, ` is not a witness`

60, ` is not a witness`

63, ` is not a witness`

74, ` is not a witness`

75, ` is not a witness`

81, ` is not a witness`

94, ` is not a witness`

97, ` is not a witness`

``(191)*``(421)*``(431)

Exercise:

3) 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 a rnage over the primes less than 20 will correctly identify all composites less than 10^15.

We can clean up the code and fine a better prime finding routine

> oddfactor := proc(n)
local temp, count:
temp := n:
count := 0:
while (temp mod 2) = 0 do
temp := temp/2:
count := count + 1:
od:
[temp, count];
end:
primetester3 := proc(n)
local i, temp, nums, count, temp2:
temp2 := n:
nums := oddfactor(n-1):
for i from 2 to 100 do
temp := Power(i,n-1) mod n:
if temp <> 1 then
if temp2 <> 0 then
print(n,` is not prime`):
print(i,` ^ `,n-1,` = `,temp):
fi:
temp2 := 0:
break:
else
for count from 1 to nums[2] do
temp := Power(i,(n-1)/2^count) mod n:
if (temp <> 1) and ((n-temp) <> 1) then
if temp2 <> 0 then
print(n,` is not prime`):
### WARNING: note that `I` is no longer of type `^`
print(i,`^`,(n-1)/2^count,` = `,temp,` and `,
### WARNING: note that `I` is no longer of type `^`
i,`^`,(n-1)/2^(count-1),`=1`):
fi:
temp2 := 0:
break:
fi:
if temp = -1 then
print(i,` ^ `,(n-1)/2^count,` = `,temp):
break:
fi:
od:
fi:
if temp2 = 0 then break: fi:
od:
if temp2 = n then print(n,` is probably prime`,i) fi:
temp2;
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 := primetester3(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); isprime(");

>

Warning, incomplete string;  use " to end the string

Exercise:

4) Use the findprime2 procedure to find 60, 70, and 80 digit primes.

>