Here is a pseudocode implementation of the Baillie-Wagstaff Pseudoprimality Test that Daniel Fischer mentioned in his comment. We begin with a simple Sieve of Eratosthenes that we will need later.
function primes(n)
ps := []
sieve := makeArray(2..n, True)
for p from 2 to n step 1
if sieve(p)
ps.append(p)
for i from p * p to n step p
sieve[i] := False
return ps
The powerMod function raises the base b to the exponent e with all calculations done modulo m; it is much faster than performing the exponentiation first, then taking the modulus of the result, because the intermediate calculation will be huge.
function powerMod(b, e, m)
x := 1
while e > 0
if e % 2 == 1
x := (b * x) % m
b := (b * b) % m
e := floor(e / 2)
return x
The jacobi function from number theory tells if a is a quadratic residue mod p.
function jacobi(a, p)
a := a % p
t := 1
while a != 0
while a % 2 == 0
a := a / 2
if p % 8 == 3 or p % 8 == 5
t := -t
a, p := p , a # swap
if a % 4 == 3 and p % 4 == 3
t := -t
a := a % p
if p == 1 return t else return 0
Gary Miller's strong pseudoprime test is based on the Little Theorem of Pierre de Fermat, which states that if p is a prime number, then for any a != 0, a ^ (p - 1) == 1 (mod p). Miller's test is somewhat stronger than Fermat's because it can't be fooled by Carmichael Numbers.
function isStrongPseudoprime(n, a)
d := n - 1; s := 0
while d % 2 == 0
d := d / 2; s := s + 1
t = powerMod(a, d, n)
if t == 1 return ProbablyPrime
while s > 0
if t == n - 1 return ProbablyPrime
t := (t * t) % n; s := s - 1
return Composite
The Miller-Rabin test performs k strong-pseudoprime tests, where k is typically somewhere between 10 and 25. The strong pseudoprime test can be fooled, but if you perform enough of them, the likelihood of being fooled is very small.
function isPrime(n) # Miller-Rabin
for i from 1 to k
a := randInt(2 .. n-1)
if not isStrongPseudoprime(n, a)
return Composite
return ProbablyPrime
That primality test is sufficient for most purposes, and fast enough. But if you want something a little bit stronger and a little bit faster, a test based on Lucas chains can be used. Here is the calculation of the Lucas chain.
function chain(n, u, v, u2, v2, d, q, m)
k := q
while m > 0
u2 := (u2 * v2) % n; v2 := (v2 * v2 - 2 * q) % n
q := (q * q) % n
if m % 2 == 1
t1 := u2 * v; t2 := u * v2
t3 := v2 * v; t4 := u2 * u * d
u, v := t1 + t2, t3 + t4
if u % 2 == 1 u := u + n
if v % 2 == 1 v := v + n
u, v, k := (u / 2) % n, (v / 2) % n), (q * k) % n
m := floor(m / 2)
return u, v, k
It is common initialize the Lucas chain using an algorithm due to John Selfridge.
function selfridge(n)
d, s := 5, 1; ds := d * s
repeat
if gcd(ds, n) > 1 return ds, 0, 0
if jacobi(ds, n) == 1 return ds, 1, (1 - ds) / 4
d, s := d + 2, s * -1; ds := d * s
Then a Lucas pseudoprime test determines if a number is prime or probably composite. Like the Fermat test, it comes in two flavors, both standard and strong, and like the Fermat test it can be fooled, though with the Fermat test the error is that a composite number could be improperly reported prime but with the Lucas test the error is that a prime number could be improperly reported composite.
function isLucasPseudoprime(n) # standard
d, p, q := selfridge(n)
if p == 0 return n == d
u, v, k := chain(n, 0, 2, 1, p, d, q, (n + 1) / 2)
return u == 0
function isLucasPseudoprime(n) # strong
d, p, q := selfridge(n)
if p == 0 return n == d
s, t := 0, n + 1
while t % 2 == 0
s, t := s + 1, t / 2
u, v, k := chain(n, 1, p, 1, p, d, q, t // 2
if u == 0 or v == 0 return Prime
r := 1
while r < s
v := (v * v - 2 * k) % n; k := (K * k) % n
if v == 0 return Prime
return ProbablyComposite
Then the Baillie-Wagstaff test is simple. First check if the input is less than 2 or is a perfect square (check if the square root is an integer). Then trial division by the primes less than 100 finds most composites quickly, and finally a strong pseudoprime test to base 2 (some people add a strong pseudoprime test to base 3, to be extra sure) followed by a Lucas pseudoprime test makes the final determination.
function isPrime(n) # Baillie-Wagstaff
if n < 2 or isSquare(n) return False
for p in primes(100)
if n % p == 0 return n == p
return isStrongPseudoprime(n, 2) \
and isLucasPseudoprime(n) # standard or strong
The Baillie-Wagstaff test has no known errors.
Once you have a good primality test, you can find the largest prime less than n by counting down from n, stopping at the first prime number.
If you are interested in programming with prime numbers, I modestly recommend this essay at my blog, or many of the other blog entries related to prime numbers, which you can find by using the search function at the blog.
