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

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