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.
[…] 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):
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:
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.
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.