programmingpraxis
-
Python,
pasted
on Feb 21:
|
# lucas pseudoprimality tester
# based on trn.c from http://www.trnicely.net/misc/bpsw.html
def gcd(a, b):
while b != 0:
a, b = b, a % b
return a
# jacobi a m -- jacobi symbol
def jacobi(a, m):
a = a % m; t = 1
while a != 0:
while a % 2 == 0:
a = a / 2
if m % 8 == 3 or m % 8 == 5:
t = -1 * t
a, m = m, a # swap a and m
if a % 4 == 3 and m % 4 == 3:
t = -1 * t
a = a % m
if m == 1:
return t
return 0
def isLucasPrime(n):
d = 5
while gcd(d, n) > 1 or jacobi(d, n) > -1:
d = (abs(d) + 2) * (-1 if d > 0 else 1)
p = 1; q = (1 - d) / 4
print "d, p, q =", d, p, q
u = 0; u2 = 1; v = 2; v2 = p; bits = (n + 1) / 2
while bits > 0:
print "u, u2, v, v2, q, bits =", u, u2, v, v2, q, bits
u2 = (u2 * v2) % n
v2 = (v2 * v2 - 2 * q) % n
q = (q * q) % n
if bits % 2 == 1:
u = u2 * v + u * v2
u = u + n if u % 2 == 1 else u
u = (u / 2) % n
v = v2 * v + d * u2 * u
v = v + n if v % 2 == 1 else v
v = (v / 2) % n
bits = bits // 2 # integer division
return u == 0
print isLucasPrime(97)
print ""
print isLucasPrime(99)
|
Output:
|
d, p, q = 5 1 -1
u, u2, v, v2, q, bits = 0 1 2 1 -1 49
u, u2, v, v2, q, bits = 1 1 54 3 1 24
u, u2, v, v2, q, bits = 1 3 54 7 1 12
u, u2, v, v2, q, bits = 1 21 54 47 1 6
u, u2, v, v2, q, bits = 1 17 54 73 1 3
u, u2, v, v2, q, bits = 38 77 18 89 1 1
False
d, p, q = 13 1 -3
u, u2, v, v2, q, bits = 0 1 2 1 -3 50
u, u2, v, v2, q, bits = 0 1 2 7 9 25
u, u2, v, v2, q, bits = 7 7 3 31 81 12
u, u2, v, v2, q, bits = 7 19 3 7 27 6
u, u2, v, v2, q, bits = 7 34 3 94 36 3
u, u2, v, v2, q, bits = 26 28 58 52 9 1
False
|
|