N+1 Primality Prover
February 7, 2014
We use the lucas
and selfridge
functions from the previous exercise:
(define (nplus1-prime? n)
(let-values (((d p q) (selfridge n)))
(if (< 1 (gcd d n)) #f
(if (positive? (u p q n (+ n 1))) #f
(let f-loop ((fs (facts (+ n 1))))
(if (null? fs) #t
(let d-loop ((p p) (q q))
(if (= (gcd (u p q n (/ (+ n 1) (car fs))) n) 1)
(f-loop (cdr fs))
(d-loop (+ p 2) (+ p q 1))))))))))
Our facts
function finds some of the factors of n using trial division, stopping when the factored portion is greater than the remaining unfactored portion. Here, f is the current trial factor, fs is the current list of factors, without duplicates, fp is the current product of the factors, and r is the remaining unfactored portion of n:
(define (facts n)
(let loop ((f 2) (fs (list)) (fp 1) (r n))
(cond ((< n (* fp fp)) (reverse fs))
((< r (* f f)) (reverse (cons r fs)))
((zero? (modulo r f))
(loop f (if (member f fs) fs (cons f fs))
(* fp f) (/ r f)))
(else (loop (+ f 1) fs fp r)))))
Facts
can be made to run faster, but if you use a method other than trial division, you will have to recursively prove prime all of factors that you find.
Here are some examples:
> (nplus1-prime? (+ #e1e12 39))
#t
> (nplus1-prime? (+ #e1e12 61))
#t
You can run the program at http://programmingpraxis.codepad.org/J9oc7Q3B.