## Primality Checking, Revisited

### January 26, 2010

It is easier to begin at the end rather than the beginning:

`(define (lucas? n)`

(let loop ((a 11) (b 7))

(let ((d (- (* a a) (* 4 b))))

(cond ((square? d) (loop (+ a 2) (+ b 1)))

((not (= (gcd n (* 2 a b d)) 1))

(loop (+ a 2) (+ b 2)))

(else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))

(m (quotient (- n (legendre d n)) 2))

(f (lambda (u) (modulo (- (* u u) 2) n)))

(g (lambda (u v) (modulo (- (* u v) x1) n))))

(let-values (((xm xm1) (chain m f g 2 x1)))

(zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))

The `loop`

tries various parameters for the Lucas sequence until it finds a pair that satisfies the constraints; note some tricky math that renders the Legendre symbol as `(gcd n (* 2 a b d))`

. Then it computes the parameters for the Lucas chain *x*_{0}, *x*_{1}, …, *x*_{m}, *x*_{m-1}; note that *x*_{0} is always 2. The last line computes the modulus and determines whether *n* is prime/pseudoprime or composite.

The function that computes the Lucas chain is straight forward, `cdr`

ing down the bits of *m*:

`(define (chain m f g x0 x1)`

(let loop ((ms (digits m 2)) (u x0) (v x1))

(cond ((null? ms) (values u v))

((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))

(else (loop (cdr ms) (g u v) (f v))))))

Now we combine the Lucas pseudoprime test with Miller’s strong pseudoprime tests to bases 2 and 3, testing for divisibility by a few small primes before starting the primary algorithm:

`(define (prime? n)`

(cond ((or (not (integer? n)) (< n 2))

(error 'prime? "must be integer greater than one"))

((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))

(else (and (miller? n 2) (miller? n 3) (lucas? n)))))

We test by calculating the Mersenne primes (Sloane’s A000043) less than a thousand:

`(define (mersenne k) (- (expt 2 k) 1))`

`> (let loop ((k 1000) (ps '()))`

(cond ((= k 1) ps)

((prime? (mersenne k))

(loop (- k 1) (cons k ps)))

(else (loop (- k 1) ps))))

(2 3 5 7 13 17 19 31 61 89 107 127 521 607)

We borrowed much from our earlier work. `Expm`

, `digits`

, and `isqrt`

come from the Standard Prelude, and `square?`

is a simple use of `isqrt`

. `Jacobi`

comes from the exercise on Modular Arithmetic, and `legendre`

is just `jacobi`

for odd primes. We’ve done the modular inverse twice; the version of inverse used here derives from the exercise on Lenstra’s Algorithm. And the `miller?`

test derives from the previous Primality Checking exercise. You can see all of the code assembled and run the program at http://programmingpraxis.codepad.org/07RDsaSw.

In python:

[…] we now have a function to compute Lucas sequences, we can take another look at testing primality using Lucas […]