```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 ``` ```# 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) / 2 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) ```
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ``` ```d, p, q = 5 1 -2 u, u2, v, v2, q, bits = 0 1 2 1 -2 49 u, u2, v, v2, q, bits = 1 1 56 5 4 24 u, u2, v, v2, q, bits = 1 5 56 17 16 12 u, u2, v, v2, q, bits = 1 85 56 63 62 6 u, u2, v, v2, q, bits = 1 20 56 62 61 3 u, u2, v, v2, q, bits = 12 76 87 36 35 1 False d, p, q = 13 1 -6 u, u2, v, v2, q, bits = 0 1 2 1 -6 50 u, u2, v, v2, q, bits = 0 1 2 13 36 25 u, u2, v, v2, q, bits = 13 13 57 97 9 12 u, u2, v, v2, q, bits = 13 73 57 85 81 6 u, u2, v, v2, q, bits = 13 67 57 34 27 3 u, u2, v, v2, q, bits = 14 1 16 13 36 1 False ```