两百万前素数之和与前两百万素数之和
By 苏剑林 | 2014-06-10 | 70628位读者 |标题说了两道比较好玩的编程题,如果读者觉得标题绕的让人眩晕的话,那么让我再说得清晰一点:
两百万前素数之和指的是所有不超过两百万的素数的和;
前两百万素数之和指的是前两百万个素数的和。
我是从子谋的blog中看到这道题目的,前一道题目是Project Euler的第10题,后一道则是我跟子谋探索着玩的。关于子谋的研究和代码,大家可以去他的blog上学习。本文分享一下我自己的想法。
两道题目本质上都是素数表的构建问题,前一题稍微简单,后一题则稍微复杂,计算量也大一点。构建素数表和判断素数都是用基本的“埃拉托斯特尼筛法”,即用2到$\sqrt{n}$的素数去除$n$。为了编程上的方便,我们通常都是2到$\sqrt{n}$的所有整数去除。但是本文所探讨的题目设计的数都很大,提高效率非常有必要。因此,我们得实现只用2到$\sqrt{n}$的素数去除$n$。
需要注意的是,不少读者会对筛法有个误解。一般来说,我们会用2到$\sqrt{n}$的素数去除$n$的方法来判断$n$的素性,然后对1到200万的数字进行判断,然后把所有素数加起来。这样是可行的,但是效率会大打折扣!要知道,最原始的埃拉托斯特尼筛法并不是这样除的,它是像下面那样表述的
埃拉托斯特尼筛法
要找到不超过n的素数,先找到2到$\sqrt{n}$的所有素数$2,3,5,\dots,p_m$,然后从1到n中依次删除2的倍数(2本身除外)、3的倍数(3本身除外)、...、$p_m$的倍数($p_m$本身除外)
也就是说,用埃拉托斯特尼筛法构建素数表,只需要乘法以及删除操作!要知道乘法的效率会比除法高很多,而且这样的算法比一个个判断所进行的运算次数也会少很多。
在编程实现上面算法的思路是,先构建一个包含1到$n$的全体整数的数组,然后依次让2的倍数的那些项为0、3的倍数的那些项为0、...、$p_m$的倍数那些项为0。可是如何得到1到$\sqrt{n}$的素数呢?其实也简单,就从前面的数组拿,也就是边判断、边收集。素数的生成速度会比需要的素数个数增加的速度要快。(这句话的意思是,比如在20以内,只需要删除2、3除本身外的倍数,就可以得到所有的素数2、3、5、7、...、19,最大的素数到了19,这对于判断400以内的素数都够用了。我们要得到2到$\sqrt{n}$的各素数$p_i$,就是要在进行数组前面的非0项取,因为合数已经被设为0了。读者只需要简单思考,就可以发现这个过程是逐步推进的。)
说到这里,编程的基本思路已经出来了。然而,还需要考虑的是这两道题目已经初步涉及到了大数运算,尤其是后一道题,前两百万个素数也就意味着大约3500万以内的素数。这时我们得考虑编程工具的问题了。我一开始想用C++,毕竟效率高,但是C++的数组大小有限制,而且整数型大小也有限制,容易溢出;于是我考虑在C++上使用GMP大数库,但是我的C++基础太弱,始终弄不好(据说GMP的作者始终不肯支持VC和VS,网上的方法都是群众修改出来的。同时想过自己编写简单的大数类,但是没有心思折腾了);也想过用matlab,但是matlab不支持高精度的大数运算,算出来的是近似值,不大适合;还想过用Mathematica,它本身就是符号运算工具,支持超大整数的运算,但问题是我发现mathematica处理几百万个元素的数组时速度也很慢。后来就锁定了Python。
很早前就尝试过Python,但是搁下了很长很长时间。今天发现Python本身支持大整数运算!虽然效率比不上C++,但是它可以方便完成本文的编程任务,而且经过优化,所用的时间也大大缩短了。在此把代码附在文后跟大家分享。
测试结果:
在Python 3.3中
两百万前的素数之和:
142913828922
time: 2.4048174478605646前两百万的素数之和:
31381137530481
time: 46.75734807838953
据说用python 2.7会更快,但在我的机器上实测,倒没有明显的区别。计算前两百万的素数之和平均时间在50秒左右,比子谋的算法快一些,原因应该就是我使用了构造的方式而不是逐个判断的方式来得到素数表,相同的是,子谋跟我都是只用2到$\sqrt{n}$之间的所有素数而不是所有整数去判断。
期待看到在Python上更高效率的算法~
200万以内的素数之和
#200万以内的素数之和
import time
start=time.clock()
import math
n=2000000 #定义上限
prime=[i for i in range(1,n+1)] #定义整数表
r=int(math.sqrt(n))
#下面是用删除式的方法把整数表中的合数删除掉
for j in range(2,r+1):
if prime[j-1] != 0:
s=j*j
while s <= n:
prime[s-1]=0
s=s+j
print(sum(prime)-1) #求和
end=time.clock()
print("time:",end-start)
前200万个素数之和
#前200万个素数之和
import time
start=time.clock()
import math
n=35000000
#定义上限,两百万个素数大约是前3500万的素数
#这是根据公式pi(n)约等于n/ln(n)得到的。
prime=[i for i in range(1,n+1)] #定义整数表
r=int(math.sqrt(n))
#下面是用删除式的方法把整数表中的合数删除掉
for j in range(2,r+1):
if prime[j-1] != 0:
s=j*j
while s <= n:
prime[s-1]=0
s=s+j
#下面开始确定素数个数到200万个。
prime.sort() #从小到大重新排列(让0位于前面)
z=prime.count(0) #统计0的个数
print(sum(prime[z+1:z+1+2000000])) #计算前两百万个素数之和。
end=time.clock()
print("time:",end-start)
转载到请包括本文地址:https://spaces.ac.cn/archives/2612
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Jun. 10, 2014). 《两百万前素数之和与前两百万素数之和 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/2612
@online{kexuefm-2612,
title={两百万前素数之和与前两百万素数之和},
author={苏剑林},
year={2014},
month={Jun},
url={\url{https://spaces.ac.cn/archives/2612}},
}
June 12th, 2014
你的计算前200万个素数和的代码在我的python2.7上用时21秒,比我的方法快太多。。
我觉得我们的方法本质是一样的,但我的用的是除法,导致计算更慢。不知是否是这个原因。
我很好奇你的python2.7为什么会这么快...我在win8上的python2.7和3.3一样慢
June 12th, 2014
我知道了,比如对于一个数49,你的方法是它作为7的倍数而直接从素数表里删除,而我的方法却是分别用2,3,5,7去除它,不但增加了计算量,用的还是除法,所以速度更慢了。
你说的应该就是主要原因,直接通过乘法构造这个表出来,会比除法省心好多。另外还有一些小原因,比如你用了两个表而我只用一个表等,但是这些都是小影响了。
我突然想到,虽然对于诸如49这样的数,你要除以2、3、5、7而导致计算量增加,但是对于30这样的数,我的方法却要将它分别当作2、3、5的倍数各删除一次,从这点来说,判断次数估计是扯平的。所以主要的差别就是乘除法的差别了。
August 22nd, 2015
第10题的thread里,第五页Lucy_Hedgehog 给了一个很快的算法
请看这里:http://kexue.fm/index.php/archives/2991/
这里已经达到0.13s了,估计已经是python下最快的了。
都一年多了,我都忘了这里有留言了。那个算法10亿以内素数的和用Python才0.5s,比你的快多了。这是知乎上对那个算法的解释
https://www.zhihu.com/question/29580448/answer/45218281
直接递归求和,不求素数表,着实巧妙~~
我觉得吧,那本文的算法作为造素数表算法来说,应该还是不错的^_^
July 16th, 2016
其实Mathematica算这两个不到1秒就够了,Python使用numpy的话也可以不到1秒:
primes[n_]:=Module[{p=Range[1,n,2]},p[[1]]=2;
Do[If[p[[(k+1)/2]]!=0,p[[(k^2+1)/2;; ;;k]]=0],{k,3,n^.5,2}];
SparseArray[p]["NonzeroValues"]];
primes[2*^6]//Tr//AbsoluteTiming
primes[InverseFunction[PrimePi][2*^6]]//Tr//AbsoluteTiming
更短的代码,但没有上面的速度快:
AbsoluteTiming@Tr@Prime@Range@PrimePi@2000000
AbsoluteTiming@Total@Prime@Range@2000000
January 23rd, 2018
用C++试了一下,
2,000,000以前的素数和只用了0.015秒,
2,000,000个素数和算法比较差,但只用了3.8秒
December 11th, 2019
筛素数表可以使用Sieve of Atkin算法,求素数和可以使用min 25筛或者Meissel-Lehmer算法应该可以更快。
December 11th, 2019
判断素数可以用Miller_Rabin素性测试法,或者是用这种方法:
素数总是形如 $6x-1$ 或者 $6x+1$,其中 $x$ 是正整数(除$2,3$外)所以循环的步长可以设为 $6$,然后每次只判断 $6$ 两侧的数即可。
```cpp
inline bool is_prime(const unsigned num)
{
if (num1;
if (num%6!=1&&num%6!=5) return false;
const unsigned lim=sqrt(num);
for (unsigned i=5;i
评论区暂不支持代码,不过我在后台能看到,理解你的意思了,谢谢~
March 27th, 2021
c++的boost.multiprecision可以做任意精度的计算