# lucas pseudoprimality test
def iSqrt(n):
if type(n) not in (int, long): raise TypeError, "must be integer"
if n < 1: raise ValueError, "must be positive"
x = n
while 1:
y = (x + n // x) // 2
if x <= y: return x
x = y
def isSquare(n):
n2 = iSqrt(n)
n == n2 * n2
def gcd(a, b):
while b != 0:
a, b = b, a % b
return a
def jacobi(a, m):
a, t = a % m, 1
while a != 0:
while a % 2 == 0:
a = a / 2
if m % 8 in (3, 5):
t = -t
a, m = m, a
if a % 4 == 3 and m % 4 == 3:
t = -t
a = a % m
if m == 1:
return t
return 0
def isStrongPseudoprime(n, a):
d, s = n - 1, 0
while d % 2 == 0:
d, s = d / 2, s + 1
t = pow(a, d, n)
if t == 1: return True
while s > 0:
if t == n - 1: return True
t, s = (t * t) % n, s - 1
return False
def selfridge(n):
dAbs, sign, d = 5, 1, 5
while 1:
if 1 < gcd(d, n): return d, 0, 0
if jacobi(d, n) == -1: break
dAbs, sign = dAbs + 2, -sign
d = dAbs * sign
return d, 1, (1 - d) / 4
def chain(n, u, v, u2, v2, d, q, m):
qkd = q
while m > 0:
u2, v2 = (u2 * v2) % n, (v2 * v2 - 2 * q) % n
q = (q * q) % n; qkd = (qkd * q) % n
if m % 2 == 1:
t1, t2, t3, t4 = u2 * v, u * v2, v2 * v, u2 * u * d
u = t1 + t2; u = u + n if u % 2 == 1 else u; u = (u / 2) % n
v = t3 + t4; v = v + n if v % 2 == 1 else v; v = (v / 2) % n
m = m // 2
return u, v, qkd
def isStandardLucasPrime(n):
d, p, q = selfridge(n)
if p == 0: return n == d
u, v, qkd = chain(n, 0, 2, 1, p, d, q, (n + 1) / 2)
return u == 0
def isStrongLucasPrime(n):
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, qkd = chain(n, 1, p, 1, p, d, q, t // 2)
if u == 0 or v == 0: return True
for r in xrange(1, s):
v = (v * v - 2 * qkd) % n
if v == 0: return True
qkd = (qkd * qkd) % n
return False
def isPrime(n):
if type(n) not in (int, long): raise TypeError, "must be integer"
if n < 2: return False
if isSquare(n): return False
for p in [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]:
if n % p == 0: return n == p
return isStrongPseudoprime(n, 2) \
and isStrongPseudoprime(n, 3) \
and isStrongLucasPrime(n)
for n in [83, 89, 111, 323, 5777, 3825123056546413051, 2**89-1]:
print n, isStrongPseudoprime(n, 2), \
isStrongPseudoprime(n, 3), \
isStandardLucasPrime(n), \
isStrongLucasPrime(n), \
isPrime(n)