这是 Daniel Fischer 在评论中提到的 Baillie-Wagstaff 伪素性测试的伪代码实现。我们从一个简单的 Eratosthenes 筛开始，稍后我们将需要它。
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
该powerMod函数将底数b提高到指数e，所有计算均以m为模；它比先求幂，然后取结果的模要快得多，因为中间计算会很大。
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
数论中的jacobi函数告诉a是否是二次余数 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 的强伪素检验基于 Pierre de Fermat 的小定理，该定理指出如果p是素数，那么对于任何a != 0，a ^ ( p - 1) == 1 (mod p )。Miller 的检验比 Fermat 的检验要强一些，因为它不会被 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
Miller-Rabin 检验执行k个强伪素检验，其中k通常介于 10 和 25 之间。强伪素检验可能会被愚弄，但如果您执行足够多的检验，被愚弄的可能性非常小。
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
该素数测试对于大多数目的来说已经足够了，而且速度也足够快。但是如果你想要更强大一点、更快一点的东西，可以使用基于卢卡斯链的测试。这是卢卡斯链的计算。
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
由于 John Selfridge，通常使用算法初始化 Lucas 链。
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
然后卢卡斯伪素测试确定一个数字是素数还是可能是合数。像费马测试一样，它有两种风格，标准和强，和费马测试一样，它可以被愚弄，虽然费马测试的错误是合数可能被错误地报告为素数，但卢卡斯测试错误是素数可能被错误地报告为复合数。
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
那么 Baillie-Wagstaff 测试很简单。首先检查输入是否小于 2 或者是完美平方（检查平方根是否为整数）。然后用小于 100 的素数进行试除法可以快速找到大多数复合物，最后对基数 2 进行强伪素数检验（有些人在基数 3 上添加了一个强伪素数检验，以确保更加确定），然后是卢卡斯伪素数检验，最终确定.
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
Baillie-Wagstaff 检验没有已知错误。
一旦你有一个好的素数测试，你可以通过从n开始倒数找到小于n的最大素数，在第一个素数处停止。
如果你对素数编程感兴趣，我谦虚地推荐我的博客上的这篇文章，或者其他许多与素数相关的博客文章，你可以使用博客上的搜索功能找到这些文章。