在Python中快速确定小于10亿的数是否为质数

Question

我的当前python质数检查算法对于在1000万和10亿之间的数字来说太慢了。我希望它能得到改进，尽管我知道我永远不会得到比10亿更大的数字。

背景是我无法得到一个足够快的实现来解决项目欧拉第60题：我在75秒内得到问题的答案，而我需要在60秒内得到它。http://projecteuler.net/index.php?section=problems&id=60

我手头可用的内存非常少，因此无法存储1亿以下的所有素数。

我目前正在使用调整为6k ± 1的标准试除法。还有比这更好的方法吗？我是否需要为这么大的数字使用Rabin-Miller方法。

primes_under_100 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
def isprime(n):
    if n <= 100:
        return n in primes_under_100
    if n % 2 == 0 or n % 3 == 0:
        return False

    for f in range(5, int(n ** .5), 6):
        if n % f == 0 or n % (f + 2) == 0:
            return False
    return True

如何改进这个算法？

需要说明的是，我是Python新手，想要只使用Python 3+。

---

最终代码

对于那些感兴趣的人，使用MAK的想法，我生成了以下代码，速度约快1/3，使我能够在不到60秒的时间内得到欧拉问题的结果！

from bisect import bisect_left
# sqrt(1000000000) = 31622
__primes = sieve(31622)
def is_prime(n):
    # if prime is already in the list, just pick it
    if n <= 31622:
        i = bisect_left(__primes, n)
        return i != len(__primes) and __primes[i] == n
    # Divide by each known prime
    limit = int(n ** .5)
    for p in __primes:
        if p > limit: return True
        if n % p == 0: return False
    # fall back on trial division if n > 1 billion
    for f in range(31627, limit, 6): # 31627 is the next prime
        if n % f == 0 or n % (f + 4) == 0:
            return False
    return True

- Olivier Grégoire

---

警告：纯Python示例（他的第一个代码片段）并不能适用于所有素数。应将for f in range(5, int(n ** .5), 6):这一行修改为for f in range(5, int(n ** .5) + 1, 6):，因为它在能够显示该数字可被其平方根整除之前就退出了（也就是退出得过早）。 - deceleratedcaviar

第二个例子有效，问题已被证明是有用的。没有理由仅因此而投反对票。我明确要求“改进算法”，并声明我是 Python 编程新手。这意味着我在这个问题上遇到了困难。仅基于这个原因投反对票违反了 SO 的宗旨（自你发表评论后两个月内有 2 次反对票，这是有关系的）。当时我提出这个问题显然适合 SO。无论如何，我会修复第一个片段中的算法。欢迎提供示例数字。 - Olivier Grégoire

@ogregoire：我不知道Daniel是否真的给你点了踩，但我觉得他的警告很有用。我几乎使用了第一个片段，因为我只需要一个快速而简单的isprime函数，而第二个片段在Python 2.x上不能直接运行 。 - Joseph Garvin

获取gmpy2，并使用gmpy2.is_prime(n)。 - xylon97

相关：列出小于N的所有质数的最快方法 - jfs

Answer 1

针对10⁹这样大的数字，一种方法是生成所有小于sqrt(10⁹)的质数列表，然后仅检查输入数字对该列表中的数字是否整除。如果一个数字不可被其平方根以下的任何质数整除，则它本身一定是一个质数（它必须至少有一个因子<=sqrt及另一个>= sqrt才能不是质数）。注意，您无需测试所有数字的可除性，只需测试到平方根即可（大约为32,000-我认为相当容易处理）。您可以使用埃氏筛法来生成质数列表。

您还可以选择概率质数测试。但它们可能更难理解，在此问题上，仅使用生成的质数列表就足够了。

- MAK

---

是的，32k个数字确实是我可以存储的。好主意。 - Olivier Grégoire

如果数字小于32k，请使用二分查找，并使用 bisect 模块。 - st0le

Answer 2

为了解决Project Euler问题，我按照您在问题中提到的建议进行操作：实现Miller Rabin测试（使用C#实现，但我认为在Python中也会很快）。该算法并不难。对于小于4,759,123,141的数字，只需检查一个数是否是2、7、61这些底数的强伪素数即可。再结合试除法和小质数即可。

我不知道您已经解决了多少问题，但拥有一个快速的素性测试将对许多问题非常有价值。

- heijp06

---

好的，那么你称小质数为什么？我应该设置什么限制？ - Olivier Grégoire

@Frór：你需要进行实验来找到最优值，但我建议从100以下的所有质数开始尝试。如果我没记错的话，甚至可能会跳过对除基数（在这种情况下为2、7、61）以外的所有值进行试除法。 - heijp06

Python：已证明在大N范围内正确 - P i

好的，4,759,123,141非常小... 可以很快地通过将其除以奇数直到平方根来检查。但感谢@Pi提供的链接——我仍然不明白为什么没有np.miller_rabin函数（或者如果这太科学了，可以使用scipy）。 - Tomasz Gandor

Answer 3

好的，我对Peter Van Der Heijden在回答中提到“流行”的Python库中没有针对真正大质数（或数字）的好方法的评论有一个后续。结果我错了 - 在sympy（众多符号代数库之一）中有一个。

https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.primetest.isprime

当然，它可能会产生误报高于10**16，但这已经比我什么都不做（除了pip install sympy ;)）得到的任何结果要好得多。

- Tomasz Gandor

---

SymPy 1.1（2017年7月）切换到BPSW，因此对于任何64位输入都没有错误的阳性结果。在某些情况下，它将使用确定性Miller-Rabin算法，但它们也已经验证了2^64。对于64位输入，唯一比这更好的方法是优化预测试以使其更快。对于更大的输入，在大多数情况下，做更多并没有令人信服的好处（需要长时间讨论）。 - DanaJ

谢谢更新！我先阅读了 sympy 0.x 的旧资料，然后链接到最新的文档。这并不改变我的观点，sympy 很棒，只是比我想象的更好 ;) - Tomasz Gandor

实际上，对于大多数情况，我认为这是正确的答案。OP正在做欧拉计划，所以这可能不适用于早期的问题，但对于后面的问题和任何其他实际用途可能会很好。 - DanaJ

Answer 4

你可以先将你的n通过primes_under_100进行除法运算。

同时，预先计算出更多的素数。

此外，你实际上是将你的range()结果存储在内存中 - 使用irange()代替并利用这个内存来运行Eratosthenes筛法算法。

- crazylammer

---

好的，我记性不差 ;) 我正在使用Python 3。我从未在Python 3中看到过xrange。 - Olivier Grégoire

xrange 在 py3k 中变成了简单的 range。 - user3850

Answer 5

def isprime(num):
if (num==3)or(num==2):
    return(True)
elif (num%2 == 0)or(num%5 == 0):
    return (False)
elif ((((num+1)%6 ==0) or ((num-1)%6 ==0)) and (num>1)):
    return (True)
else:
    return (False)

我认为这段代码是最快的...

- Priyadarsan Priyadarsan

---

每个质数（除了2和3）都可以表示为6n（+/-）1的形式。 - Priyadarsan Priyadarsan

在我的机器上，检查所有小于1000000的质数只需0.42194461822509766秒。函数体内没有循环或迭代。 - Priyadarsan Priyadarsan

是的，但反过来并不一定成立 - 可以表示为6n+-1的每个数字都不一定是质数。例如，25不是质数（6 * 4 + 1）。 - JJJ

感谢您发现了这个错误！很抱歉浪费了您的时间。我会修改代码的。不过，我只是一个自学的业余程序员。 - Priyadarsan Priyadarsan

我感谢您的指导。 - Priyadarsan Priyadarsan