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, cdring 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

Posted by programmingpraxis
Filed in Exercises
2 Comments »

2 Responses to “Primality Checking, Revisited”

  1. Mike said
    April 12, 2010 at 11:55 PM

    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. Lucas Pseudoprimes | Programming Praxis said
    February 4, 2014 at 9:00 AM

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

Leave a comment