Primality Checking, Revisited

January 26, 2010

[ There is a bug in the solution of this exercise. See the revised version of the exercise for a proper solution. ]

We examined the Miller-Rabin probabilistic primality checker in a previous exercise. Today, we examine a primality checker that combines the Miller-Rabin test with a test on Lucas pseudoprimes, devised by Robert Baillie and described by Baillie and Wagstaff in their article “Lucas Pseudoprimes” in Mathematics of Computation, Volume 35, Number 152, pages 1391-1417, October 1980; see also Thomas Nicely’s web page devoted to Baillie’s test. This is the same algorithm used in the PrimeQ function in Mathematica.

Lucas numbers are defined by a recurrence formula similar to Fibonacci numbers, where L_n = L_{n-1} + L_{n-2} with L_1 = 1 and L_2 = 3; the Lucas numbers are 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, … (Sloane’s A000204). Lucas numbers have the rather startling property that, if n is prime, L_n \equiv 1 \pmod{n}. But the converse is not true, and composite numbers n such that L_n \equiv 1 \pmod{n} are known as Lucas pseudoprimes; the first few Lucas pseudoprimes are 705, 2465, 2737, 3745, 4171, … (Sloane’s A005845).

Lucas numbers are a special case of Lucas sequences. If P and Q are integers such that the discriminant D = P^2 - 4Q, then the roots of x^2 - P x + Q = 0 are a = \frac{P + \sqrt{D}}{2} and b = \frac{P - \sqrt{D}}{2}. There are two Lucas sequences U_n(P, Q) = \frac{a^n - b^n}{a - b} and V_n(P, Q) = a^n + b^n for n \geq 0, which can be computed by the recurrence equations U_m(P,Q) = P U_{m-1}(P,Q) - Q U_{m-2}(P,Q) and V_m(P,Q) = P V_{m-1}(P,Q) - Q V_{m-2}(P,Q). The Lucas numbers are given by the sequence V(1,-1) and the Fibonacci numbers are given by the sequence U(1,-1).

We will want to compute the nth element of a Lucas sequence, mod n, for large n. Rather than computing the entire recurrence in time proportional to n, it is possible to use doubling and halving, in the same way as the exercise on Three Binary Algorithms, to compute the nth element in log n time. Such a computation is known as a Lucas chain.

Thus, to test whether an odd number n is a Lucas pseudoprime, we choose sequence parameters P and Q so that the the discriminant D is non-square and the Legendre number \left( \begin{array}{c} D \\ n\end{array} \right) = -1 (otherwise the modular arithmetic would fail). Then we construct the Lucas chain; if the nth element is zero modulo n, then n is either prime or is a Lucas pseudoprime.

Recall that the Miller-Rabin test used strong pseudoprime tests on fifty bases to check primality. Using Lucas pseudoprimes, we can reduce the number of tests substantially. It turns out that combining two strong pseudoprime tests, with bases 2 and 3, with a Lucas pseudoprime test, there are no known pseudoprimes; the test has been performed exhaustively on all numbers less than 1016, and no counter-examples have been found in over twenty-five years of use (in addition to fifteen minutes of fame, there is a monetary reward from the original authors for the finder of a counter-example).

Your task is to write a function that checks primality using the algorithm described above. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

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 600 other followers

%d bloggers like this: