Lucas-Lehmer素性检验的快速位模算法

Question

Lucas-Lehmer primality test测试质数以确定它们是否也是Mersenne primes。其中一个瓶颈是(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 (包括Python3)？有没有更好的方法？

Answer 1

我能找到的最好的改进是从模函数中完全删除Mp = (1<<p) - 1，并在开始L-L测试的迭代之前在L-L函数中预先计算它。使用while n > Mp:代替while n.bit_length() > p:也节省了一些时间。

Answer 2

在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

编辑是因为我之前搞砸了格式；感谢本杰明指出这一点。启示:不要从空闲窗口复制粘贴到SO中。抱歉的!

(注意:将n的长度减半而不是去掉p的标准，以及k1的确切选择，都有一点错误，但这并不重要，所以我没有费心修改它们。)

如果我使用p=12345和n=9**200000 (是的，我知道p不是素数，但这在这里无关紧要)，那么速度大约快13倍。

遗憾的是，这对没有帮助，因为在L-L测试中，n永远不会大于(2^p)^2。对不起。