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.
Programming Praxis 
[…] 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. […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/02/02/programming-praxis-proving-primality/ for a version with comments):
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 qsI 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 FalseIt 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.
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.
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