Lucas可能素数检验

Question

我已经尝试了几天来实现Baillie-PSW素数检验，并遇到了一些问题。当尝试使用Lucas可能素数检验的时候。我的问题不是关于Baile，而是关于如何生成正确的Lucas序列模某些数字。

对于前两个psudoprimes，我的代码给出了正确的结果，例如对于323和377。然而，对于下一次的psudoprime，标准实现和加倍版本都失败了。

试图在V_1上执行模块化操作完全破坏了Luckas序列生成器的加倍版本。

对于如何在Python中正确实现Lucas可能的素数测试，有什么建议或建议吗？

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

这是我的卢卡斯伪随机性测试；您可以在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)

注意，1159是一个已知的卢卡斯假犯罪(A217120)。