algorithm - 如何实现Frobenius伪素算法？

Question

有人告诉我，Frobenius伪素算法的运行时间比Miller-Rabin素数测试的时间长三倍，但分辨率却高出七倍。因此，如果一个地方要运行前十次，然后再运行三十次，那么两者都将花费相同的时间运行，但是前者将提供大约233％的分析能力。为了找出执行测试的方法，发现了以下论文，最后给出了该算法：

A Simple Derivation for the Frobenius Pseudoprime Test

尝试实现以下算法，但该程序从不打印数字。可以让更熟悉数学符号或算法的人验证发生了什么吗？

编辑1：下面的代码已添加更正，但缺少compute_wm_wm1的实现。有人可以从算法的角度解释递归定义吗？对我来说，这不是“点击”。

编辑2：错误的代码已被删除，并在下面添加了compute_wm_wm1函数的实现。它似乎可行，但可能需要进一步优化才能实用。

from random import SystemRandom
from fractions import gcd
random = SystemRandom().randrange

def find_prime_number(bits, test):
    number = random((1 << bits - 1) + 1, 1 << bits, 2)
    while True:
        for _ in range(test):
            if not frobenius_pseudoprime(number):
                break
        else:
            return number
        number += 2

def frobenius_pseudoprime(integer):
    assert integer & 1 and integer >= 3
    a, b, d = choose_ab(integer)
    w1 = (a ** 2 * extended_gcd(b, integer)[0] - 2) % integer
    m = (integer - jacobi_symbol(d, integer)) >> 1
    wm, wm1 = compute_wm_wm1(w1, m, integer)
    if w1 * wm != 2 * wm1 % integer:
        return False
    b = pow(b, (integer - 1) >> 1, integer)
    return b * wm % integer == 2

def choose_ab(integer):
    a, b = random(1, integer), random(1, integer)
    d = a ** 2 - 4 * b
    while is_square(d) or gcd(2 * d * a * b, integer) != 1:
        a, b = random(1, integer), random(1, integer)
        d = a ** 2 - 4 * b
    return a, b, d

def is_square(integer):
    if integer < 0:
        return False
    if integer < 2:
        return True
    x = integer >> 1
    seen = set([x])
    while x * x != integer:
        x = (x + integer // x) >> 1
        if x in seen:
            return False
        seen.add(x)
    return True

def extended_gcd(n, d):
    x1, x2, y1, y2 = 0, 1, 1, 0
    while d:
        n, (q, d) = d, divmod(n, d)
        x1, x2, y1, y2 = x2 - q * x1, x1, y2 - q * y1, y1
    return x2, y2

def jacobi_symbol(n, d):
    j = 1
    while n:
        while not n & 1:
            n >>= 1
            if d & 7 in {3, 5}:
                j = -j
        n, d = d, n
        if n & 3 == 3 == d & 3:
            j = -j
        n %= d
    return j if d == 1 else 0

def compute_wm_wm1(w1, m, n):
    a, b = 2, w1
    for shift in range(m.bit_length() - 1, -1, -1):
        if m >> shift & 1:
            a, b = (a * b - w1) % n, (b * b - 2) % n
        else:
            a, b = (a * a - 2) % n, (a * b - w1) % n
    return a, b

print('Probably prime:\n', find_prime_number(300, 10))

Answer 1

由于不熟悉该符号，您似乎完全误解了该算法。

def frobenius_pseudoprime(integer):
    assert integer & 1 and integer >= 3
    a, b, d = choose_ab(integer)
    w1 = (a ** 2 // b - 2) % integer

从线来

W0≡2（mod n）和W1≡a2b−1 − 2（mod n）

但是b-1在这里不是 1/b的意思，而是 b模 n的模逆，即带有 c的整数 b·c ≡ 1 (mod n)。通过 c的连续小数展开，或者通过 extended Euclidean algorithm等效地进行计算，您可以最容易地找到这样的 b/n。由于您可能不熟悉连续分数，因此建议您使用后者。

    m = (integer - d // integer) // 2

来自

n −（∆ / n）= 2m

并将 Jacobi symbol误认为是分数/分区（诚然，我在这里将其显示为甚至是分数），但是由于该站点不支持LaTeX渲染，因此我们必须这样做。
Jacobi符号是Legendre符号的概括-用相同的符号表示-表示数字是否为模为奇数质数的二次余数（如果 n为模为 p的二次余数，即存在 k为< cc>和 k^2 ≡ n (mod p)不是 n的整数倍，然后是 p，如果 (n/p) = 1是 n的整数倍，那么 p，否则是 (n/p) = 0）。 Jacobi符号解除了对“分母”为奇数质数的限制，并允许任意奇数作为“分母”。它的值是Legendre符号的乘积，所有除 (n/p) = -1的质数均按相同的“分子”（根据多重性）。除此之外，还有如何在链接的文章中有效地计算Jacobi符号。
该行应正确读取

m = (integer - jacobi_symbol(d,integer)) // 2

我完全无法理解以下几行，从逻辑上讲，这里应该遵循
Wm和Wm + 1使用递归

W2j≡Wj2-2（mod n）
W2j + 1≡WjWj + 1-W1（mod n）

围绕PDF的公式（11）给出了使用该递归计算所需值的有效方法。

    w_m0 = w1 * 2 // m % integer
    w_m1 = w1 * 2 // (m + 1) % integer
    w_m2 = (w_m0 * w_m1 - w1) % integer

该函数的其余部分几乎是正确的，当然，由于先前的误解，它现在会得到错误的数据。

    if w1 * w_m0 != 2 * w_m2:

这里的（不）等式应为模 n，即 integer。

        return False
    b = pow(b, (integer - 1) // 2, integer)
    return b * w_m0 % integer == 2

但是请注意，如果 if (w1*w_m0 - 2*w_m2) % integer != 0是质数，则

b^((n-1)/2) ≡ (b/n) (mod n)

其中， n是Legendre（或Jacobi）符号（对于素数“分母”，Jacobi符号是Legendre符号），因此是 (b/n)。因此，可以将其用作额外的检查，如果Wm不是2或 b^((n-1)/2) ≡ ±1 (mod n)，则 n-2不能是素数，如果 n不是1或 b^((n-1)/2) (mod n)也不能。
大概先计算 n-1并检查它是1还是 b^((n-1)/2) (mod n)是个好主意，因为如果该检查失败（顺便说一下，这就是 Euler pseudoprime测试），您将不需要另一个，同样也很昂贵，计算，如果成功，很可能仍然需要计算。
关于更正，它们似乎是正确的，除了一个使我先前忽略的故障可能更糟的故障：

if w1 * wm != 2 * wm1 % integer:

这仅将模数应用于 n-1。
关于Wj的递归，我认为最好以可行的实现方式进行解释，首先是为了方便复制和粘贴：

def compute_wm_wm1(w1,m,n):
    a, b = 2, w1
    bits = int(log(m,2)) - 2
    if bits < 0:
        bits = 0
    mask = 1 << bits
    while mask <= m:
        mask <<= 1
    mask >>= 1
    while mask > 0:
        if (mask & m) != 0:
            a, b = (a*b-w1)%n, (b*b-2)%n
        else:
            a, b = (a*a-2)%n, (a*b-w1)%n
        mask >>= 1
    return a, b

然后在两者之间进行解释：

def compute_wm_wm1(w1,m,n):

我们需要W1的值，所需数字的索引以及将模数用作输入的数字。值W0始终为2，因此我们不需要将该值作为参数。
称其为

wm, wm1 = compute_wm_wm1(w1,m,integer)

在 2 * wm1中（除了：不是一个好名字，大多数返回 frobenius_pseudoprime的数字都是实质数）。

    a, b = 2, w1

我们分别将 True和 a初始化为W0和W1。在每个点上， b保持Wj的值，而 a保持Wj + 1的值，其中 b是到目前为止消耗的 j的位的值。例如，对于 m， m = 13， j和 a的值如下所示：

consumed remaining  j    a    b
           1101     0   w_0  w_1
   1        101     1   w_1  w_2
   11        01     3   w_3  w_4
   110        1     6   w_6  w_7
   1101            13  w_13  w_14

这些位从左到右被消耗，因此我们必须找到 b的第一个设置位并将“指针”放在它的前面

    bits = int(log(m,2)) - 2
    if bits < 0:
        bits = 0
    mask = 1 << bits

我从计算的对数中减去了一点，以确保我们不会被浮点错误所迷惑（顺便说一句，使用 m将您限制为最多1024位，大约308个十进制数字；如果如果您想处理更大的数字，则必须以其他方式找到 log的以2为底的对数，使用 m是最简单的方式，并且这只是概念的证明，因此我在这里使用了它。

    while mask <= m:
        mask <<= 1

移动遮罩，直到它大于 log，这样设置位就指向 m的第一个设置位之前。然后向后移动一个位置，所以我们指向钻头。

    mask >>= 1
    while mask > 0:
        if (mask & m) != 0:
            a, b = (a*b-w1)%n, (b*b-2)%n

如果设置了下一位，则 m的已消耗位的初始部分的值将从 m变为 j，因此我们需要的W序列的下一个值是 2*j+1和W2j的W2j + 1 a +2。通过以上递归公式，

W_{2j+1} = W_j * W_{j+1} - W_1 (mod n)
W_{2j+2} = W_{j+1}^2 - 2 (mod n)

由于 b为Wj， a为Wj + 1，因此 b变为 a，而 (a*b - W_1) % n变为 b。

        else:
            a, b = (a*a-2)%n, (a*b-w1)%n

如果未设置下一位，则 (b * b - 2) % n的已消耗位的初始部分的值将从 m变为 j，因此 2*j变为W2j =（Wj2-2）（mod n），并且< cc>变成
W2j + 1 =（Wj * Wj + 1-W1）（mod n）。

        mask >>= 1

将指针移到下一位。当我们经过最后一位时， a变为0，循环结束。现在， b的已消耗位的初始部分是 mask的所有位，因此该值当然是 m。
那我们可以

    return a, b

一些补充说明：

def find_prime_number(bits, test):
    while True:
        number = random(3, 1 << bits, 2)
        for _ in range(test):
            if not frobenius_pseudoprime(number):
                break
        else:
            return number

在较大的数字中，素数不太常见，因此仅选择随机数很可能会尝试很多次。如果您选择一个随机数并按顺序检查候选者，则可能会更快地找到素数（或可能的素数）。
另一个要点是，像Frobenius检验这样的检验对于发现例如3的倍数是复合的。在使用这样的测试（或Miller-Rabin测试，Lucas测试或Euler测试，...）之前，您绝对应该做一些试验划分以淘汰大部分复合材料，并且仅在以下情况下进行工作：它有值得一战的机会。
哦， m函数还没有准备好处理小于2的参数，那里潜伏着零除错误，

def is_square(integer):
    if integer < 0:
        return False
    if integer < 2:
        return True
    x = integer // 2

应该有所帮助。