需要帮助实现Lucas伪光学性测试

Question

我试图编写一个函数来确定一个数字n是素数还是复合函数，使用Lucas伪反式测试；目前，我正在使用标准测试，但是一旦我完成了这个测试，我就会编写强测试。我正在阅读贝利和瓦格斯塔夫的纸，并在trn.c文件中很好地遵循托马斯的实现。

据我所知，整个测试包括几个步骤:用小素数进行试验，检验n不是正方形，对基2进行强伪反性检验，最后再进行Lucas伪凸检验。我可以处理所有其他的部分，但我对卢卡斯的假冒伪劣测试有困难。下面是我在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，我知道这是一个卢卡斯伪粒子，isLucasPrime应该返回True。实际上，isLucasPseudoprime为我测试过的每一个n返回False。

我有几个问题：

1)我不是C/GMP方面的专家，但在我看来，从右到左(从最不重要到最重要)的(n+1)/2部分很好地运行，而其他作者则从左到右运行。上面显示的代码从左到右运行，但我也尝试过从右到左运行，结果是一样的。哪种顺序是正确的？

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

3)我做错了什么？

Answer 1

1)我不是C/GMP方面的专家，但在我看来，从右到左(从最不重要到最重要)的(n+1)/2部分很好地运行，而其他作者则从左到右运行。上面显示的代码从左到右运行，但我也尝试过从右到左运行，结果是一样的。哪种顺序是正确的？

实际上，从最不重要的一点到最重要的一点都是很好的。他在U(2^k)和mpzV2m变量中计算mpzU2m和V(2^k) (当然都是模块化N )，并将U((N+1) % 2^k) resp 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，并根据前缀p后面的下一个位是0还是1，将它们结合起来计算U(2*p+1)和V的U(2*p)或U(2*p+2) ditto。

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

我--和其他人，显然--觉得从左到右的顺序更简单。我没有做过，也没有读过一篇分析文章，可能是从右到左的计算成本平均较低，因此很好地选择了从右到左。

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

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

3)我做错了什么？

第一个小错误：

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

应该是

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

当然了。

更重要的是，您使用更新的良好的遍历顺序(从右到左)，但遍历在相反的方向。这当然会产生错误的结果。

此外，在更新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

有效(没有保证，但我认为它是正确的，并做了一些测试，所有这些都通过了)。