标题说了两道比较好玩的编程题,如果读者觉得标题绕的让人眩晕的话,那么让我再说得清晰一点:

两百万前素数之和指的是所有不超过两百万的素数的和;
前两百万素数之和指的是前两百万个素数的和。

我是从子谋的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://kexue.fm/archives/2612

更详细的转载事宜请参考:《科学空间FAQ》

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (Jun. 10, 2014). 《两百万前素数之和与前两百万素数之和 》[Blog post]. Retrieved from https://kexue.fm/archives/2612

@online{kexuefm-2612,
        title={两百万前素数之和与前两百万素数之和},
        author={苏剑林},
        year={2014},
        month={Jun},
        url={\url{https://kexue.fm/archives/2612}},
}