Lucas 可能素数测试

Question

我已经尝试了几天来实现Baillie-PSW素性测试，但遇到了一些问题。特别是在尝试使用Lucas probable prime test时。我的问题不是关于Baile，而是如何生成模某个数字的正确Lucas序列 对于前两个伪质数，我的代码给出了正确的结果，例如323和377。然而对于下一个伪质数，标准实现和加倍版本都失败了。

尝试对V_1进行模运算会完全破坏Luckas序列生成器的加倍版本。

有没有关于如何在Python中正确实现Lucas probable prime test的提示或建议？

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)

- N3buchadnezzar

---

从您提供的文章中可以看到：“如果其中任意一个分子是奇数，我们可以通过将其增加n来使其变为偶数，因为所有这些计算都是在模n下进行的。” 您尝试过这个方法吗？ - Lynn

谢谢！现在luckas_sequence_doubling返回与luckas_sequence_standard相同的值，但它们仍然显示不正确的值。例如，说1159不是伪素数。我应该更新我的问题来修复错误吗？ - N3buchadnezzar

当然你应该 :) - Lynn

Answer 1

这是我的Lucas伪素性测试；您可以在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)。

- user448810