## 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 x0, x1, …, xm, xm-1; note that x0 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.

Pages: 1 2

### 2 Responses to “Primality Checking, Revisited”

1. Mike said

In python:

```from math import sqrt

def bits_of( n ):
'''generator that returns the bits in n, from MSB to LSB.'''

bit = 1
while bit < n:
bit <<= 8
while bit > n:
bit >>= 1

while bit:
yield 1 if ( bit & n ) else 0
bit >>= 1

def euclid(a, b):
'''Extended Euclid's algorithm.'''
lastx, x = 1, 0
lasty, y = 0, 1

while b:
q = a/b
a, b = b, a%b
lastx, x = x, lastx - q*x
lasty, y = y, lasty - q*y

return lastx, lasty, a

def gcd(m, n):
'''Return greatest common divisor or m and n.'''

if m < n:
m,n = n,m

while n:
m, n = n, m % n

return m

def inverse(x, m):
a, b, g = euclid(x, m)
if g != 1:
raise ValueError("x and m must be coprime.")

return a % m

def is_square( xx ):
if xx&7==1 or xx&31==4 or xx&127==16 or xx&191==0:
x = sqrt( xx )

return x == int( x )

return False

def jacobi(a, n):
if a == 0:
return 0
elif a == 1:
return 1
elif a == 2:
return 1 if (n%8) in (1,7) else -1
elif a&1 == 0:
return jacobi(2,n) * jacobi(a/2, n)
elif n < a:
return jacobi(a%n, n)
elif a%4 == 3 and n%4 == 3:
return -jacobi(n, a)
else:
return jacobi(n, a)

legendre = jacobi       # well, for odd primes anyway

def lucas_test( n ):
a, b = 11, 7
while True:
d = a * a - 4 * b
if is_square(d):
a += 2,
b += 1

elif gcd( n, 2 * a * b * d ) != 1:
a += 2
b += 2

else:
break

x1 = ( a * a * inverse(b, n) - 2 ) % n
m = ( n - legendre(d, n) ) / 2

u = 2
v = x1
for bit in bits_of( m ):
if bit:
u = ( u * v - x1 ) %n
v = ( v * v - 2 ) % n
else:
v = ( u * v - x1 ) %n
u = ( u * u - 2 ) % n

return ( x1 * u - 2 * v ) % n == 0

def miller_rabin_test(n, a):
'''checks if n is prime base a.

returns True is a is likely a composite.
'''
r, s = 0, n-1
while s&1 == 0:
r += 1
s >>= 1

if pow(a, s, n) == 1:
return True      # possibly prime

for j in range(r):
if pow(a, s, n) == n-1:
return True  # possibly prime
s *= 2

return False  # definitely composite

small_primes = set([  2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ])

def is_prime(n):
if n < 100:
return n in small_primes

return all(n%p for p in small_primes) and \
miller_rabin_test( n, 2 ) and \
miller_rabin_test( n, 3 ) and \
lucas_test( n )

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