Project:
 ```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 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 ``` ```# 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) ```
 ```1 2 3 4 5 6 7 ``` ```83 True True True True True 89 True True True True True 111 False False False False False 323 False False True False False 5777 False False True True False 3825123056546413051 True True False False False 618970019642690137449562111 True True True True True ```