## Proving Primality

### February 2, 2010

[ Updated 14 April 2010 to fix a bug. ]

We use the trial-division function `td-factors` from the exercise on wheel factorization and `expm` from the Standard Prelude:

```(define (prime? n)   (if (even? n) (= n 2)     (let* ((n-1 (- n 1)) (fs (td-factors n-1)))       (let loop1 ((b 2))         (cond ((= b n) #f)               ((= (expm b n-1 n) 1)                 (let loop2 ((qs fs))                   (cond ((null? qs) #t)                         ((= (expm b (/ n-1 (car qs)) n) 1)                           (loop1 (+ b 1)))                         (else (loop2 (cdr qs))))))               (else (loop1 (+ b 1))))))))```

As an example, we compute the prime factors of 289-1. The prime factors of n-1 are 2, 3, 5, 17, 23, 89, 353, 397, 683, 2113, and 2931542417, and it is 2(n-1)/89 ≡ 4096 (mod n) that provides the counterexample that proves 289-1 is prime:

```> (prime? (- (expt 2 89) 1)) #t```

Since the first version of this program had a bug, we provide the following test program, which uses the `primes` function from a previous exercise:

```(define (test-prime? n)   (let ((ps (primes n)))     (do ((i 2 (+ i 1))) ((= i n))       (let ((p? (prime? i)))         (when (and p? (not (member i ps)))           (display i) (display " found to be prime") (newline))         (when (and (not p?) (member i ps))           (display i) (display " found to be composite") (newline))))))```

Calling `(test-prime? 1000)` shows no errors.

The speed of this algorithm depends primarily on the speed by which n-1 can be factored, so it is better to replace the factorization by trial division with one of the other factorization functions that we have studied. But if you do that, be sure that if you use one of the algorithms that can return a non-prime factor, you don’t use a probabilistic prime checker to determine whether or not the factor is prime!

You can run this program at http://programmingpraxis.codepad.org/1ttTLsZF.

Pages: 1 2

### 6 Responses to “Proving Primality”

1. […] Praxis – Proving Primality By Remco Niemeijer In today’s Programming Praxis exercise we have to implement an algorithm to prove the primality of a number. […]

2. Remco Niemeijer said

```import Data.Bits
import Data.List

tdFactors :: Integer -> [Integer]
tdFactors n = f n [2..floor . sqrt \$ fromIntegral n] where
f _ []     = []
f r (x:xs) | mod r x == 0 = x : f (div r x) (x:xs)
| otherwise    = f r xs

expm :: Integer -> Integer -> Integer -> Integer
expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
filter (flip testBit 0 . snd) .
zip (iterate (flip mod m . (^ 2)) b) .
takeWhile (> 0) \$ iterate (`shiftR` 1) e

isPrime :: Integer -> Bool
isPrime n = f 2 \$ tdFactors (n - 1) where
f _ []     = False
f b (q:qs) | expm b (n - 1)         n /= 1 = f (b + 1) (q:qs)
| expm b (n - 1 `div` q) n /= 1 = True
| otherwise                     = f b qs
```
3. Mike said

I found the description of the algorithm confusing. When I implemented it, the routine would report some multiples of 5 (e.g., 15 and 35) to be prime.

Studying the proof at the mathpages.com link, given above, I think the algorithm should be to find a b so that b^(x-1) == 1 (mod n) then, if b^((n-1)/q) != 1 (mod n) for *all* prime divisors q of n-1, n is prime. If b^((n-1)/q) == 1 (mod n) for some q, then try a new value for b. This matches the pseudocode at wikipedia.

My python code is below:

```from itertools import chain, cycle

def factors_of(n):
f = 2
for step in chain([0,1,2,2], cycle([4,2,4,2,4,6,2,6])):
f += step

if f*f > n:
if n != 1:
yield n
break

if n%f == 0:
yield f
while n%f == 0:
n /= f

def is_prime(n):
'''proves primality of n using Lucas test.'''
x = n-1

if n == 2 or n%2 == 0:
return False

factors = list( factors_of(x))

b = 2
while b < n:
if pow( b, x, n ) == 1:

for q in factors:
if pow( b, x/q, n ) == 1:
break
else:
return True

b += 1

return False
```
4. programmingpraxis said

It wasn’t confusing, it was wrong. Thanks for pointing out the bug. I’ve updated the solution and added a test function.

By the way, your version of is_prime incorrectly reports that 2 is composite. You should change the test at lines 22 and 23 to read:

```if n%2 == 0:     return n == 2```

Don’t feel bad; Jon Bentley got that wrong in one of his books, so you’re in good company.

5. Mike said

D’oh…guess that’ll teach me not to edit a routine and not test it again before posting.
Is there a way to edit my earlier post fix lines 22 and 23? Or should I just repost a correct version?

My confusion was over the word “any” in

“… if b^((n-1)/q) != 1 (mod n) for *any* prime divisor q of n-1, n is prime.”

At first, I read the phrase to mean if I found any one q for which the non-congruence was true, then n is prime. However, the correct reading is that the non-congruence holds no matter which q you try. Therefore, you have to try every q to prove n is prime. Perhaps using “every” instead of “any” would be clearer.

6. Mike said

Here’s a corrected version of the python code for proving a number is prime using the Lucas algorithm.

```
from itertools import chain, cycle

def factors_of(n):
f = 2
for step in chain([0,1,2,2], cycle([4,2,4,2,4,6,2,6])):
f += step

if f*f > n:
if n != 1:
yield n
break

if n%f == 0:
yield f
while n%f == 0:
n /= f

def is_prime(n):
'''proves primality of n using Lucas test.'''
x = n-1

if n%2 == 0:
return n == 2

factors = list( factors_of(x))

b = 2
while b < n:
if pow( b, x, n ) == 1:

for q in factors:
if pow( b, x/q, n ) == 1:
break
else:
return True

b += 1

return False
```