Lucas-Lehmer素性测试的更快按位取模方法

Question

Lucas-Lehmer素性测试是一种用来检测素数的方法，特别是用来判断它们是否是梅森素数。在这个计算过程中，有一个比较耗时的步骤就是计算 (s**2 − 2) % (2**p - 1) 这个模运算。

如果使用位运算，可以大大加快这个过程（具体可以参考L-L的链接），我目前找到的最佳方法是：

def mod(n,p):
    """ Returns the value of (s**2 - 2) % (2**p -1)"""
    Mp = (1<<p) - 1
    while n.bit_length() > p: # For Python < 2.7 use len(bin(n)) - 2 > p
        n = (n & Mp) + (n >> p)
    if n == Mp:
        return 0
    else:
        return n

一个简单的测试案例是，当 p 有5到9位数字，而 s 有超过10,000位数字（或者更多，具体数字不重要）。可以通过 mod((s**2 - 2), p) == (s**2 - 2) % (2**p -1) 来测试解决方案。需要注意的是，在L-L测试中，需要进行p - 2次这样的模运算，每次的 s 都是指数级增长，因此需要进行优化。

有没有办法进一步加快这个过程，使用纯Python（包括Python 3）？有没有更好的方法？

Answer 1

当 n 的长度远远超过 2^p 时，你可以通过做一些事情来避免一些耗时的平方级别计算，像这样：

def mod1(n,p):
    while n.bit_length() > 3*p:
        k = n.bit_length() // p
        k1 = k>>1
        k1p = k1*p
        M = (1<<k1p)-1
        n = (n & M) + (n >> k1p)
    Mp = (1<<p)-1
    while n.bit_length() > p:
        n = (n&Mp) + (n>>p)
    if n==Mp: return 0
    return n

[编辑：之前格式搞错了，感谢 Benjamin 指出这个问题。教训是：不要从 Idle 窗口直接复制粘贴到 SO。抱歉！]

（注意：将 n 的长度减半而不是直接减去 p 的标准，以及 k1 的具体选择都有点不太对，但这没关系，所以我就没去修正它们。）

如果我设定 p=12345 和 n=9**200000（是的，我知道 p 不是质数，但在这里并不重要），那么这个方法大约快了 13 倍。

不幸的是，这对你没有帮助，因为在 L-L 测试中，n 的值从来不会超过大约 (2^p)^2。抱歉。

Answer 2

我找到的最佳改进方法是完全把 Mp = (1<<p) - 1 从模运算的函数中去掉，然后在开始 L-L 测试的迭代之前，在 L-L 函数中先计算好它。用 while n > Mp: 代替 while n.bit_length() > p: 也节省了一些时间。