本专栏会在以后的文章中专门更新欧拉计划(Project Euler)的算法题目。
欧拉计划
欧拉计划是一个发布一系列算法问题的网站,题目的类型会比较偏重于数学,以及算法效率。欧拉计划网站由Colin Hughes在2001年创立,并持续更新问题,目前网站共有700+问题,注册用户超过一百万。欧拉计划的
问题的难度从5%(最简单)到100%(最困难)不等,但是无论是什么难度的题目,欧拉计划中的问题都有一个核心的标准:当前CPU的技术条件下,单核CPU计算出每一个问题的时间均小于一分钟。因此欧拉计划中的题目主要目的不是用程序如何实现,而是如何更快。并且欧拉计划中的问题不仅需要从语法和逻辑需要优化,还需要从数学的角度进行优化。漂亮并且高效的解答,不仅仅需要编程水平,还需要数学水平。
关于欧拉计划的系列文章会涉及程序语言,比如Python,Wolfram|mathematica等,也会涉及算法的知识,以及数学知识(主要是组合与数论)。
作为正式开始的第一篇文章,这里挑选了欧拉计划中最简单的5%难度的两道题目来进行解答,一方面读者可以熟悉Python以及Wolfram语言的基本语法,第二可以感受以下欧拉计划的题目风格。在今后的文章中会更新难度更大的题目,也会就题目中涉及的数学问题以及算法的效率进行分析。
PE1问题:
第一题要求寻找所有1000以下,含有因子3或者5的数字。因为这个问题比较简单,而设立时间又大致二十年前,因此我们可以试着把问题难度加大,变成:
寻找以内,所有含有素因子,3或者5的数的总和。
要用代码实现它并不困难,最简单的想法是,对于1到
Python3:
sum
Mathematica:
Total
(对于Mathematica的代码千万要注意,mma列表的计数是从1开始,而python是从0开始。而题目要求的是10^10以下的数字,而10^10本身是被5整除的,因此在计算的时候,mma的代码要去掉最后的10^10)
Mathematica还有一种写法,稍微更好一点,就是把3,5的倍数都生成出来,用集合的求并运算生成列表,最后求和。因为这里取得是倍数,稍微列表会短一些。
Mathematica:
Union
不过如果运行这些串代码,几乎都会出现一个问题:内存不足的报错。
这是因为,为了计算满足条件的数字的总和的时候先生成了这些数字的一个列表,而这样的列表占用了过多的内存导致内存不足。不过我们想要的问题的结果只是这些数的总和是多少,因此可以只是求和而没有必要存储究竟是哪些数字满足条件。于是可以这样写:
Python 3:
ans
然而这样的作法虽然解决了内存的问题,但时间上仍然无法符合一分钟原则(事实上运算时间远远超出一分钟。实际时间是:3123.14秒,大致一个小时)。
但如果,仔细分析问题并从数学上进行优化,那么这个问题的时间复杂度可以达到
首先所有能被3整除的自然数,自然是3的倍数,因此具有
因此不超过n的所有3的倍数的总和是:
对于n而言,这个数字必然大于等于不超过n的最大的3的倍数
同样的道理,不超过n的5的倍数,的总和为
但是这样做还有一个漏洞需要修复,重复计算问题。因为例如15这样的,15的倍数,在计算3的倍数的时候被算入其中,而计算5的倍数的时候再次重复计算,因此我们必须要减去重复计算的部分。道理很简单,则来自于集合论的事实:
假设其中一个集合代表不超过n的所有3的倍数的集合,而另一个是所有不超过n5的倍数的集合。而这两个集合的交集就是不超过n所有15的数的集合。要正确的计算3,5倍数的集合,就应该是两个集合的并减去两个集合的交集。(因为直接相加实际上,把交集计算了两次。)
在集合论上,一般的如果由m个集合,想要正确的计算元素的个数,这个叫做容斥原理(Inclusion–exclusion principle)。
因此计算结果就应该是:
其中:
不过上述结果的前提是不超过n,但是问题要求是n以下,如果
因此最后的结果是:,。
而这个结果好算到,就算是用纸和笔也能得到想要的结果。
在Python里面可以用 //整除符号代替取下整floor函数。但是注意在Mathematica里面 // 表示是在 // 后面加上函数名字,可以嵌套一个函数。比如Table[ ]//Total 对列表求和。
Python 3:
def
Mathmatica:
f
PE问题3:
这个问题非常好概括,就是求一个数的最大素因子。
这个问题可以归结为一个一般的问题:
对于任意给定的自然数,如何分解这个自然数为素因子的乘积。
这个问题本身是一个相当困难的问题,因为如果这个问题能够高效率的解决,那么相当于证明了:
但是,该题目中给出的数字600851475143本身算不得什么大数(按照当前计算机的运算能力至少要是155位的十进制数,也就是一个512比特的数字,才可以称之为大数),因此我们并不需要真的去解决一般的问题,只要对较小的数字有一些最基本的提升运算效率的技巧即可。
其次对于Mathematica而言已经有很好的内置函数去解决分解较小的自然数,所以对于该问题只需要:
Mathematica:
FactorInteger
(//在MMA中的功能和Python有很大的区别。)
就可以知道最大的素因子,但是这样没有任何意思,这里我们自己写一个算法来解决问题。
这个问题所涉及的两个要素:素数,因子。
因此最朴素的思想莫过于,判断600851475143以下的所有自然数,判断其是否为素数,如果是素数是否能整除600851475143,如果能,就留下来,最后选出最大的素因子。
但这样实际上有很多判断是不必要的。
因为对于任意自然数n,如果n是合数,假设
也就是说,只需要判断一个数n,小于等于
比如,分解100:
由于5是属于{2,3,5,7}中的一个,所以,100的分解就是:
因此想要分解600851475143我们首先要知道,小于等于
这个不妨提前生成,而生成素数所需要的算法无非是素数判定算法,而这个素数判定算法效率很高,大概
from
这样便生成了我们所需要的素数。
对于上述的寻找100的最大素因子的情况做一个分析,
整个过程可以表述为(此处假设primes = {2,3,5,7}):
1.从primes中最小的素数2开始,取素因子p,比如
2.把
上述的算法为了最快解决问题,我们可以把它写成一个递归的函数。
不过这个时候如果我们需要使得我们的函数能够针对于更佳一般的情况,我们要考虑到一些特殊情形。
调整递归深度上限!
为什么需要调整递归深度上限?毕竟不是所有正整数总是有不均匀的素因子的。比如
2 x 224737这种素因子大小分布非常不均匀的情况,恰好是非常好计算的,因为只需要剥离2,剩下的立刻就是最大素因子。然而如果随机给我们的递归函数输入一个数,可能会出现,50508067591 =224737*24743的数。这个数并不算大,甚至没有该题目要求的数字大。但是两个素因子都是
所以如果要剥离像上面例子中的自然数的第一个素因子需要递归深度大于20000,而Python的递归深度默认10000,mathematica更是只有1024。所以我们需要调整系统默认的递归深度。
from
做成递归的定义,这样函数运行起来的算法效率更高。实际上对于该问题的运行时间是,大约0.04s。并且该函数考虑了一般情形的n,而不仅仅只是该问题中的n。
Mathematica:
(完全一样的算法,但是函数不同,风格也不同)
(*MMA不需要生成primes列表,因为自带有第i的素数的函数Prime[i]*)
不过这样的递归算法是有明显的缺陷的,对于即便不太大的特殊的数字,但是两个素因子都不小,比如:50508067591 =224737 * 224743 ,运行上述代码会出现递归超过上限。
最后感谢@早生华发 提供的部分代码以及对整体,包括Python3以及Mathematica的代码的修改建议。