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 10¹⁶, 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.

Posted by programmingpraxis

Filed in Exercises

2 Comments »

2 Responses to “Primality Checking, Revisited”

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 )

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

Programming Praxis