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.

About these ads

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 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 643 other followers

%d bloggers like this: