Lucas probable prime test

Question

I have been trying to implement the Baillie-PSW primality test for a few days, and have ran into some problems. Sepcifically when trying to use the Lucas probable prime test. My question is not about Baile, but on how to generate the correct Lucas sequence modulo some number

For the first two psudoprimes my code gives the correct result, eg for 323 and 377. However for the next psudoprime, both the standard implementation and the doubling version fails.

Trying to do modulo operations on V_1 completely breaks the doubling version of the Luckas sequence generator.

Any tips or suggestions on how to correctly implement the Lucas probable prime test in Python?

from fractions import gcd
from math import log

def luckas_sequence_standard(num, D=0):
    if D == 0: 
        D = smallest_D(num) 

    P = 1
    Q = (1-D)/4

    V0 = 2
    V1 = P

    U0 = 0
    U1 = 1

    for _ in range(num):
        U2 = (P*U1 - Q*U0) % num
        U1, U0 = U2, U1

        V2 = (P*V1 - Q*V0) % num
        V1, V0 = V2, V1  

    return U2%num, V2%num


def luckas_sequence_doubling(num, D=0):
    if D == 0:
        D = smallest_D(num) 
    P = 1
    Q = (1 - D)/4

    V0 = P
    U0 = 1

    temp_num = num + 1
    double = []
    while temp_num > 1:
        if temp_num % 2 == 0:
            double.append(True)
            temp_num //= 2
        else:
            double.append(False)
            temp_num  += -1

    k = 1
    double.reverse()
    for is_double in double:
        if is_double:

            U1 = (U0*V0) % num
            V1 = V0**2 - 2*Q**k 

            U0 = U1
            V0 = V1

            k *= 2

        elif not is_double:

            U1 = ((P*U0 + V0)/2) % num
            V1 = (D*U0 + P*V0)/2

            U0 = U1
            V0 = V1

            k += 1
    return U1%num, V1%num


def jacobi(a, m):
    if a in [0, 1]:
        return a
    elif gcd(a, m) != 1:
        return 0
    elif a == 2:
        if m % 8 in [3, 5]:
            return -1
        elif m % 8 in [1, 7]:
            return 1
    if a % 2 == 0:
        return jacobi(2,m)*jacobi(a/2, m)
    elif a >= m or a < 0:
        return jacobi(a % m, m)
    elif a % 4 == 3 and m % 4 == 3:
        return -jacobi(m, a)
    return jacobi(m, a)


def smallest_D(num):
    D = 5
    k = 1
    while k > 0 and jacobi(k*D, num) != -1:
        D += 2
        k *= -1
    return k*D


if __name__ == '__main__':

    print luckas_sequence_standard(323)
    print luckas_sequence_doubling(323)
    print 
    print luckas_sequence_standard(377)
    print luckas_sequence_doubling(377)
    print 
    print luckas_sequence_standard(1159)
    print luckas_sequence_doubling(1159)

Answer 1

Here is my Lucas pseudoprimality test; you can run it at ideone.com/57Iayq.

# lucas pseudoprimality test

def gcd(a,b): # euclid's algorithm
    if b == 0: return a
    return gcd(b, a%b)

def jacobi(a, m):
    # assumes a an integer and
    # m an odd positive integer
    a, t = a % m, 1
    while a <> 0:
        z = -1 if m % 8 in [3,5] else 1
        while a % 2 == 0:
            a, t = a / 2, t * z
        if a%4 == 3 and m%4 == 3: t = -t
        a, m = m % a, a
    return t if m == 1 else 0

def selfridge(n):
    d, s = 5, 1
    while True:
        ds = d * s
        if gcd(ds, n) > 1:
            return ds, 0, 0
        if jacobi(ds, n) == -1:
            return ds, 1, (1 - ds) / 4
        d, s = d + 2, s * -1

def lucasPQ(p, q, m, n):
    # nth element of lucas sequence with
    # parameters p and q (mod m); ignore
    # modulus operation when m is zero
    def mod(x):
        if m == 0: return x
        return x % m
    def half(x):
        if x % 2 == 1: x = x + m
        return mod(x / 2)
    un, vn, qn = 1, p, q
    u = 0 if n % 2 == 0 else 1
    v = 2 if n % 2 == 0 else p
    k = 1 if n % 2 == 0 else q
    n, d = n // 2, p * p - 4 * q
    while n > 0:
        u2 = mod(un * vn)
        v2 = mod(vn * vn - 2 * qn)
        q2 = mod(qn * qn)
        n2 = n // 2
        if n % 2 == 1:
            uu = half(u * v2 + u2 * v)
            vv = half(v * v2 + d * u * u2)
            u, v, k = uu, vv, k * q2
        un, vn, qn, n = u2, v2, q2, n2
    return u, v, k

def isLucasPseudoprime(n):
    d, p, q = selfridge(n)
    if p == 0: return n == d
    u, v, k = lucasPQ(p, q, n, n+1)
    return u == 0

print isLucasPseudoprime(1159)

Note that 1159 is a known Lucas pseudoprime (A217120).