需要帮助实施 Lucas Pseudoprimality 测试

Question

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

我知道完整的测试涉及几个步骤：用小素数进行除法，检查 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 会为 83 和 89 等素数返回 False，这是不正确的。它还为复合 111 返回False，这是正确的。它为复合 323 返回 False，我知道它是 Lucas 伪素数，isLucasPrime 应该返回 True。事实上，isLucasPseudoprime 为我测试过的每个 n 返回 False。

我有几个问题：

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

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

3) 我做错了什么？

Answer 1

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

确实，Nicely 从最低有效位到最高有效位。他在mpzU2m 和mpzV2m 变量中计算U(2^k) 和V(2^k)（还有Q^(2^k)；当然都是模N），并在U((N+1) % 2^k) 和V((N+1) % 2^k) 中存储了mpzU 和mpzV。当遇到 1 位时，余数 (N+1) % 2^k 会发生变化，mpzU 和 mpzV 也会相应更新。

另一种方法是计算 U(p)、U(p+1)、V(p) 和（可选）V(p+1) 以获得 N+1 的前缀 p，并将它们组合起来计算 U(2*p+1) 和 U(2*p)或 U(2*p+2) [V 同上] 取决于前缀 p 之后的下一位是 0 还是 1。

这两种方法都是正确的，就像你可以计算从左到右的幂 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 变量看起来很奇怪。它是否正确？我希望每次循环都更新所有四个 Lucas 链变量，因为链的索引在每一步都会增加。

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

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

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