# 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
        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
            qkd = (qkd * q) % 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 [79, 83, 89, 111, 323, 5777, 3825123056546413051, 2**89-1]:
    print n, isStrongPseudoprime(n, 2), \
             isStrongPseudoprime(n, 3), \
             isStandardLucasPrime(n), \
             isStrongLucasPrime(n), \
             isPrime(n)

79 True True True True True
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