# 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)

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