1

我正在尝试编写一个函数，该函数使用卢卡斯伪素数测试确定数字n是素数还是复合数；目前，我正在使用标准测试，但是一旦我开始工作，我就会编写强测试。我正在阅读 Baillie 和 Wagstaff 的论文，并遵循Thomas Nicely 在trn.c文件中的实现。

我知道完整的测试涉及几个步骤：用小素数进行除法，检查n是否不是正方形，对基数 2 执行强伪素数测试，最后是卢卡斯伪素数测试。我可以处理所有其他部分，但卢卡斯伪素测试遇到了麻烦。这是我在 Python 中的实现：

def gcd(a, b):
    while b != 0:
        a, b = b, a % b
    return a

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):
    dAbs, sign, d = 5, 1, 5
    while 1:
        if 1 < gcd(d, n) > n:
            return False
        if jacobi(d, n) == -1:
            break
        dAbs, sign = dAbs + 2, sign * -1
        d = dAbs * sign
    p, q = 1, (1 - d) / 4
    print "p, q, d =", p, q, d
    u, v, u2, v2, q, q2 = 0, 2, 1, p, q, 2 * q
    bits = []
    t = (n + 1) / 2
    while t > 0:
        bits.append(t % 2)
        t = t // 2
    h = -1
    while -1 * len(bits) <= h:
        print "u, u2, v, v2, q, q2, bits, bits[h] = ",\
               u, u2, v, v2, q, q2, bits, bits[h]
        u2 = (u2 * v2) % n
        v2 = (v2 * v2 - q2) % n
        if bits[h] == 1:
            u = u2 * v + u * v2
            u = u if u % 2 == 0 else u + n
            u = (u / 2) % n
            v = (v2 * v) + (u2 * u * d)
            v = v if v % 2 == 0 else v + n
            v = (v / 2) % n
        if -1 * len(bits) < h:
            q = (q * q) % n
            q2 = q + q
        h = h - 1
    return u == 0

当我运行它时，isLucasPrime返回False83 和 89 等素数，这是不正确的。它还返回False复合 111，这是正确的。它返回False复合 323，我知道它是一个isLucasPrime应该返回的 Lucas 伪素数True。事实上，我测试过的每个nisLucasPseudoprime都会返回。False

我有几个问题：

1）我不是 C/GMP 方面的专家，但在我看来，Nicely 似乎(n+1)/2从右到左（从最不重要到最重要）的位运行，而其他作者从左到右运行位。上面显示的代码从左到右运行位，但我也尝试从右到左运行位，结果相同。哪个顺序是正确的？

2) 我觉得 Nicely 只更新1 位的u和v变量看起来很奇怪。它是否正确？我希望每次通过循环更新所有四个卢卡斯链变量，因为链的索引在每一步都会增加。

3）我做错了什么？

4

1 回答 1

3

1）我不是 C/GMP 方面的专家，但在我看来，Nicely 似乎(n+1)/2从右到左（从最不重要到最重要）的位运行，而其他作者从左到右运行位。上面显示的代码从左到右运行位，但我也尝试从右到左运行位，结果相同。哪个顺序是正确的？

事实上，Nicely 从最低有效位到最高有效位。他在and变量中计算U(2^k)and V(2^k)（和Q^(2^k)；当然都是模数N），并将resp存储在and中。当遇到 1 位时，余数会发生变化，并相应地更新 和。mpzU2mmpzV2mU((N+1) % 2^k)V((N+1) % 2^k)mpzUmpzV(N+1) % 2^kmpzUmpzV

另一种方法是计算U(p), U(p+1),V(p)和 ( 可选地)的V(p+1)前缀并将它们组合起来计算and或[ditto for ] 取决于前缀后面的下一位是 0 还是 1。pN+1U(2*p+1)U(2*p)U(2*p+2)Vp

这两种方法都是正确的，就像您可以计算x^N从左到右，具有x^p和x^(p+1)作为状态的功率，或者从右到左具有x^(2^k)和x^(N % 2^k)作为状态的功率[并且，计算U(n)和U(n+1)基本上是计算ζ^n位置ζ = (1 + sqrt(D))/2]。

我——显然还有其他人——发现从左到右的顺序更简单。我没有完成或阅读分析，可能是从右到左的平均计算成本较低，因此很好地选择了从右到左。

2) 我觉得很奇怪，Nicely 只更新1 位的u和变量。v它是否正确？我希望每次通过循环更新所有四个卢卡斯链变量，因为链的索引在每一步都会增加。

是的，这是正确的，因为(N+1) % 2^k == (N+1) % 2^(k-1)如果该2^k位为 0，则余数。

3）我做错了什么？

首先是一个小错字：

if 1 < gcd(d, n) > n:

应该

if 1 < gcd(d, n) < n:

当然。

更重要的是，您将更新用于 Nicely 的遍历顺序（从右到左），但沿另一个方向遍历。这当然会产生错误的结果。

此外，更新时v

    if bits[h] == 1:
        u = u2 * v + u * v2
        u = u if u % 2 == 0 else u + n
        u = (u / 2) % n
        v = (v2 * v) + (u2 * u * d)
        v = v if v % 2 == 0 else v + n
        v = (v / 2) % n

你使用 的新值u，但你应该使用旧值。

def isLucasPrime(n):
    dAbs, sign, d = 5, 1, 5
    while 1:
        if 1 < gcd(d, n) < n:
            return False
        if jacobi(d, n) == -1:
            break
        dAbs, sign = dAbs + 2, sign * -1
        d = dAbs * sign
    p, q = 1, (1 - d) // 4
    u, v, u2, v2, q, q2 = 0, 2, 1, p, q, 2 * q
    bits = []
    t = (n + 1) // 2
    while t > 0:
        bits.append(t % 2)
        t = t // 2
    h = 0
    while h < len(bits):
        u2 = (u2 * v2) % n
        v2 = (v2 * v2 - q2) % n
        if bits[h] == 1:
            uold = u
            u = u2 * v + u * v2
            u = u if u % 2 == 0 else u + n
            u = (u // 2) % n
            v = (v2 * v) + (u2 * uold * d)
            v = v if v % 2 == 0 else v + n
            v = (v // 2) % n
        if h < len(bits) - 1:
            q = (q * q) % n
            q2 = q + q
        h = h + 1
    return u == 0

有效（不能保证，但我认为它是正确的，并且已经做了一些测试，所有这些都通过了）。

于 2013-02-24T17:08:00.690 回答