Essays/Primality Tests

Various methods are presented for testing whether a positive integer n is prime. The result is 1 if the number is prime or is highly probably prime, 0 if it is not prime, or _1 if the method failed to determine its status.

Reference: Eric Bach and Jeffrey Shallit, Algorithmic Number Theory, Volume I: Efficient Algorithms, MIT Press, 1996, Chapter 9.

p:

1 p: n is 1 if n is prime or is highly probably prime, and 0 if it is not prime.

   1 p: i.20
0 0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1
   I. 1 p: i.20
2 3 5 7 11 13 17 19

   1 p: _1+2^67x
0

   m=: 366384x * */ p: i.9x
   1 p: 6171054912832631x + m * i.25
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Trial Division

n is prime if it is not divisible by any prime less than or equal to its square root. Trial division can be applied to eliminate "easy" composites before more expensive methods are used.

TrialDiv=: 3 : 0 " 0
 for_p. i.&.(p:^:_1) (%:+2&<) y<.2^31 do.
  if. 0=p|y do. 0 return. end.
 end.
 (1<y)*_1^y>2^31
)

   TrialDiv i.20
0 0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1
   1 p: i.20
0 0 1 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1

   TrialDiv 2147483659x
_1
   TrialDiv *: 46349x
_1

Miller-Rabin

MillerRabin is deterministic when n<3.4155e14 and probabilistic if n>:3.4155e14 .  m MillerRabin n is 0 or 1.  If 0 , then n is not prime; if 1 , then n is prime with an error probability of at most 0.25^m (the error probability is 0 if n<3.4155e14 ). Deterministic MillerRabin is faster than trial division for all but small values of n .

The magic numbers in witnesses are from http://primes.utm.edu/prove/prove2_3.html .

huo=: <.@-:^:(0=2&|)^:a:  NB. halve until odd

witnesses=: 4 : 0
 r=. 9080191 4759123141 2152302898747 3474749660383  341550071728321x
 s=. 31 73;  2 7 61;    2 3 5 7 11;   2 3 5 7 11 13; 2 3 5 7 11 13 17
 if. y>:{:r do.
  1+x ?@$ y-1
 else.
  ((r-1) I. y){:: s
 end.
)

MillerRabin=: 100&$: : (4 : 0) " 0
 if. 0=2|y do. 2=y return. end.
 w=. x witnesses y
 if. (100>y)*.y e. w do. 1 return. end.
 e=. huo y-1
 for_a. w do. if. (+./c=y-1) +: 1={:c=. a y&|@^ e do. 0 return. end. end.
 1
)

For example:

   (MillerRabin -: 1&p:) 1000 ?@$ 2e9
1
   MillerRabin n=: */ p: 1e8+0 100x
0
   MillerRabin _1+2^127 131x
1 0

Lucas-Lehmer

The Lucas-Lehmer test for Mersenne numbers states that if p is an odd prime, then m=._1+2x^p is prime if and only if  (_2+*:)^:(p-2) 4x is divisible by m . Thus:

   LucasLehmer=: 3 : '0 = (<:2x^y)&|@(_2+*:)^:(y-2) 4x' " 0

   ] b=. LucasLehmer p: i.50
0 1 1 1 0 1 1 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
   p: I. b
3 5 7 13 17 19 31 61 89 107 127
   ] m=. _1 + 2x ^ p: I. b
7 31 127 8191 131071 524287 2147483647 2305843009213693951 618970019642690137449562111
      162259276829213363391578010288127 170141183460469231731687303715884105727
   1 p: m
1 1 1 1 1 1 1 1 1 1 1

   b -:&}. 1 p: _1 + 2x ^ p: i.50
1

Regarding the Mersenne number _1+2^67 , Bach and Shallit relate the following anecdote (page 7):

In 1903 F.N. Cole factored $2^{67}-1$. When E.T. Bell, the historian of mathematics, asked Cole how long it had taken him to find the factorization, he replied, "Three years of Sundays." Bell goes on to say

At the October, 1903 meeting in New York of the American Mathematical Society, Cole had a paper on the program with the modest title On the factorization of large numbers. When the chairman called on him for his paper, Cole -- who was always a man of very few words -- walked to the board, and, saying nothing, proceeded to chalk up the arithmetic for raising 2 to the sixty-seventh power. Then he carefully subtracted 1. Without a word he moved over to a clear space on the board and multiplied out, by longhand,

    193,707,721 x 761,838,257,287.

The two calculations agreed ... For the first and only time on record, an audience of the American Mathematical Society vigorously applauded the author of a paper delivered before it. Cole took his seat without having uttered a word. Nobody asked him a question.

   ] n=: _1+2^67x
147573952589676412927
   q: n
193707721 761838257287
   */ q: n
147573952589676412927

   LucasLehmer 67
0

Collected Definitions

TrialDiv=: 3 : 0 " 0
 for_p. i.&.(p:^:_1) (%:+2&<) y<.2^31 do.
  if. 0=p|y do. 0 return. end.
 end.
 (1<y)*_1^y>2^31
)

huo=: <.@-:^:(0=2&|)^:a:  NB. halve until odd

witnesses=: 4 : 0
 r=. 9080191 4759123141 2152302898747 3474749660383  341550071728321x
 s=. 31 73;  2 7 61;    2 3 5 7 11;   2 3 5 7 11 13; 2 3 5 7 11 13 17
 if. y>:{:r do.
  1+x ?@$ y-1
 else.
  ((r-1) I. y){:: s
 end.
)

MillerRabin=: 100&$: : (4 : 0) " 0
 if. 0=2|y do. 2=y return. end.
 w=. x witnesses y
 if. (100>y)*.y e. w do. 1 return. end.
 e=. huo y-1
 for_a. w do. if. (+./c=y-1) +: 1={:c=. a y&|@^ e do. 0 return. end. end.
 1
)

LucasLehmer=: 3 : '0 = (<:2x^y)&|@(_2+*:)^:(y-2) 4x' " 0

---

Contributed by Roger Hui, with improvements to MillerRabin suggested by John Randall.