乍听之下,不无道理;仔细揣摩,胡说八道

0%

可耻的蛮力法

ProjectEuler的第66题的传送门。话说今天挺累的,去实习单位那边做了一下小小的阶段性总结,然后走回了学校的另一个门,吃完饭又走回学校另一头的我自己的宿舍,别提有多累了。看来今天这样的运动量也够了,今晚的跑步就先搁置了吧。好了,回到正题,也就是ProjectEuler的第66题。第66题乍看起来似乎很简单,无非就是查找一个满足条件的所有x当中的最大值而已。于是像几乎所有的其它题目一样,我尝试使用蛮力法来解决这个题目,而为了保险起见,我为自己的查找单个d所对应的最小的x的函数的查找过程设置了上界,因此保证了每一个函数调用都会结束。然后就开始对所有小于1000的d进行查找了,一开始判断一个数是不是平方数的方法有误,计算出了结果却是错的。待我修改了之后重新运行了一次程序,结果更令我大吃一惊——对于d为61的情况,在上限之下找不到结果——尼玛我设的上界可是一亿啊!(最近对一亿特别有感情)

抛弃蛮力法重新上路

再试了几次之后,我确定自己的方法是没有错的,所以真正的答案就是对于有些d而言结果确实无法在上界之下找到(最后计算到了真正的答案之后我才发现自己用蛮力法是多么的天真无邪可爱不懂事……)。于是我放弃了自己的方法,转而在网上寻找别人的求助。在Google中输入关键字后,很快就有英文的博客出现在了候选项中。为了不放过任何的学习机会,我点击了第一个。进去之后才发现,这个问题原来还真不简单……

数学的逆袭

ProjectEuler的第66题是一道不折不扣的数学题。在第一个博客中,我看到了博主提到了丢番图方程以及两个我不太懂的单词convergent和continue fraction。这时候还没啥感觉,只是对博主说他以0.x级别的速度做出来这道题目感到震惊和难以置信,然后又转而去看了Google到的第二个链接。说实在的,真的非常感谢这个链接的博主,因为正是他的指引才让我知道了一个关于连分数的极好的讲解,看完了博主的文章后大致知道了这道题目是怎么一回事,然后就一头扎进了那个讲解里面去看了。好吧,为了做题我又回去研究数学了。

大致的数学内容

首先题目中出现的方程是所谓的“丢番图方程”,并且可以归为更具体的类别叫做“佩尔方程”,幸运的是,对于佩尔方程的解是一个已经研究过的问题,并且根据维基百科上的说明,当d不为平方数时,佩尔方程对每一个d都有解,这是由拉格朗日证明了的。所以即使是使用蛮力法,也尽可以放心地去从1开始遍历所有的自然数——你总会找到那个心仪的它的!除了证明了有解之外,佩尔方程的解法也是被人们研究过的,并且有几种不同的解法。可惜的是本人的英语水平和数学水平都不是很强,只能看懂中文维基百科中关于佩尔方程的用连分数求解方法,因此也决定了使用这个方法来解决第66题。

重点来了:怎么计算一个数的算术平方根的连分数形式?

方法来源在前面提及的关于连分数的极好的讲解那里,不过重点其实在于用程序实现这个东西。一开始的时候我觉得这个东西可能要用到符号计算,然后就开始在那里纠结要用怎样的方式来表达一个数的算术平方根才好。到后来我盯了那个表格很久之后我才发现我错了,根本不是酱紫的!压根就不需要什么符号计算,只是纯粹的数值计算就足够了。这一切的一切的关键,就在于怎么把一个数拆分成一个整数加上一个未知数的倒数的形式。

(为了方便起见用了原文中的例子)首先,如果要计算√14,那么按照算法(请自己看上面的链接里的英文原文)的第一步,要找到一个数m,使得m的平方是一个最接近14的完全平方数。显然这个数是9,但是怎么算呢?可以一个个地从1开始找(太慢了一边去!),另一个方法是用计算器求平方根的功能直接计算得到3.7416575,然后对这个结果取不大于它的最大整数(也就是所谓的floor),得到了3。那么就可以把√14表示成3+1/x1的形式。

这样我们得到了√14的一个等级的连分数近似形式3/1(也就是3啊)。显然3*3-14*1*1得到-5,所以连分数3/1不是所希望的解,因此需要继续计算第二级的连分数近似,那么要怎么算呢?我们只要把上式中的x1也表示成连分数的形式就可以了,但是必须先求出x1的某个数值表示形式,不然无法应用上面的第一步。方法是把x1用√14和其它数来表示。于是有

x1 = 1/(√14-3) = (√14+3)/(√14-3)(√14+3) = (√14+3)/5

按照前面的方法,floor((√14+3)/5)的结果为1,所以x1可以写成1+1/x2的形式,那么√14就可以写成3+1/(1+1/x2)的形式,忽略未知的1/x2就得到了连分数4/1,仍然不是方程x^2-Dy^2=1的解。总之继续这么进行下去,就可以得到满足上面的方程的整数解了,并且其中的x必定是所有解当中最小的一个。

规律在哪?紧盯表格!

表格的第一列是Finding,显然,每一行的Finding的值决定了在Step 1中计算得到的红色的数字是多少。如果知道了每一行的Finding是多少,就有办法计算出那些红色的数字了。在Finding中出现的x1、x2和x3,其实就是前一行最后所计算出来的那些值罢了。不难发现x1、x2和x3其实都有相同的结构,那就是它们都是(√14-a)/b的形式。而a和b分别为3、2、2和5、2、5。如果再想想就会知道,其实一开始的√14也满足这样的形式。对于√14而言,a为0而b为1。既然有了规律和相应的符号表示,那么就可以从代数的角度来找规律了。加入要计算算术平方根的数是n,红色的要计算的数为m,那么在第i次计算中

(√n+a_(i-1))/b_(i-1) = m_i + 1/x_i

显然要计算的就是m_i和x_i,而x_i其实就是由a_i和b_i所构成的,因此变换x_i就可以得到

x_i = (√n-(a_(i-1) - m_ib_(i-1)))/((n-(a_(i-1) - m_ib_(i-1))^2)/b_(i-1))

上面的公式太丑了,LaTeX版本在这里。这样也就知道了对应于这个xi的a和b分别为多少了,参见前面给出的LaTeX所描述的公式的图片的链接。这些公式实际上是递推公式,也就是说你必须有了前一个值之后才能计算出下一个值。如果是希望快速计算出指定的连分数的一个部分的话,那么这种方式不适合,不过用在这里却刚好,因为我不知道第几个连分数的部分才是解,所以一定是从第一个近似连分数一直迭代过去的。

知道了计算方法,那么剩下的编码过程也就简单了很多了,详情请参见我的解法

后记

这一次的收获还真不少。首先是关于题目本身的。以前我做ProjectEuler的题目,几乎都是直接使用蛮力法来解决的,不过这一次蛮力法彻底失败了,因为消耗的时间实在是太久了所以没有必要完全计算出来了(实际上我没有完全计算,因为设置了上限在中途就发现问题了)。而上网搜索才发现这其实是一道彻头彻尾的数学题——当然了,还是需要写程序来算的。也强烈地让我认识到了,PE的题目看中的,真的如之前在豆瓣上的讨论所看到的,注重的是数学能力,而不是运用一些经典的算法的能力。第二个,是关于所谓的“过早优化是万恶之源”的说法的感想。

一开始我并不知道用蛮力法来计算是不实际的,因此我一直在尝试着改进蛮力法的性能——是的,我还没有确定这个东西是否能够正确地计算出结果的时候我已经开始考虑优化的事情了。问题出在哪里呢?就在于我对于所遇到的很多题目都是尝试使用蛮力法来解决的,以致于我有了这样的思维定势:每当我发现运行时间过长时,我就会觉得一定是我的代码性能不好,因此会尝试对其进行优化,例如加上类型声明或者使用记忆化之类的。可是偏偏这样走上了歪路,因为正确的做法是分析一下问题,找到正确的解法,而不是在错误的解法上不断地尝试优化。如果我一直优化下去,那么永远不会得到正确的答案——过早优化是万恶之源~

第三条——糟糕,洗完澡忘记了第三条是什么了!>_<(总算想起来了,于23:36)在这一次的解法中,对于连分数的构造是一个递推的过程,实际上就有种可以用递归来解决的味道。刚想着怎么表示这么一个构造连分数的过程的时候,我忽然灵光一闪,想到了只看过却从来没有使用过的技术——把闭包当作数据结构来使用。于是根据自己看《SICP》的经验,我弄了个叫做make-cf的函数,只要传递一个整数给它,那么就可以得到相应的一个闭包,对这个闭包做不同的操作就有不同的结果,而且那些状态信息还可以由闭包自己保持。等我用起来之后才发现,用闭包好像是一个很有用的想法,因为

这样子相当具有复用的可能!

代码中的make-cf、update-cf、to-cf和get-cf-list,完全可以写入一个package中然后作为一个独立的模块来提供,这样子以后的一些程序也可以使用了——PE上好像还有不少用到连分数的题目呢^_^虽然用CLOS完全可以实现这里提及的闭包的功能,可是比起CLOS提供的面向对象功能,闭包反而更有封装的感觉。总之在CLOS和闭包中,我更喜欢这种用闭包来实现的风格;-)

今天上百度贴吧的ProjectEuler吧时发现大家在讨论第26题的一个解法的原理[求值]求解释原理,然后话题很快又从对解法的原理的解释转移到了对性能的攀比之中——这也是人之常情啦,哪个做PE题目的人不希望自己的方法是尽量快的呢^_^看了一下吧主的解法的运行时间,发现居然不是毫秒级别的,顿时对自己以前的解法需要运行多久心生好奇,于是在已经打开的SLIME中加载了以前的源文件,一运行,发现相当地快,速度上是吧主给出解法的快一百倍,连我自己都吓了一跳,不过还是掩盖不了自己内心的欢喜。

我的思路

我不知道吧主是怎么做的,他在贴吧里面没有给出代码,不过我倒是可以给一下我的解法的链接第26题。刚翻出来看的时候我自己也看蒙了,完全不知道自己当初究竟是怎样的思路的,也不知道自己当初是在怎样的情况下想出这样的代码的。不过再仔细看了一会儿之后,总算有点明白是什么意思了。这其中的做法,还得从小学的时候学习除法来说起……

哦,不好意思,说太远了,没必要追溯到小学那里去,不过思路本身可能确实就是小学的除法,只不过是找了一下其中的规律,然后用计算机代替了手动计算而已。考虑一个简单的情况,例如一除以二,要怎么计算呢?嗯,因为1比2要小,所以计算结果应该是个小数,所以商必然有0.(零以及小数点),然后在1后面添一个0,再继续计算。这时候我们知道了2×5=10,所以在商的部分再添上一个5,就得到了余数为0的情况。因为余数为0,所以我们知道2是可以整除1的,所以计算就到此结束了,因此1/2=0.5。

上面的例子是可以整除的情况,因此在小数点后面并没有循环的部分,所以循环部分的长度理所当然的就是0。那如果是一除以六呢?计算方法是一样的,因为1比6小,所以商添上0.,而给1添上一个0得到10,而10除以6的商为1余数为4。接下来给4添上一个0得到40,再进行40除以6的计算过程。最后得到的结果是一个无限循环小数,按照第26题的写法就是0.1(6)。不过且慢,这里有个地方值得注意,那就是对于得到余数为4的情况,我们为其添上了一个0再继续进行计算,这个步骤和一开始为1添上0再进行计算难道有什么本质的差别吗?

显然没有啊!而且计算40除以6的结果是6.(6),如果直接把这个计算1除以6时的第一步得到的商进行拼接,就得到了0.1(6)啊。也就是说,小学的时候学到的手动进行除法计算的方法,实际上可以分解成一个递归的算法。具体如何分解不是这里的讨论重点,这里的重点在于怎样从递归的过程中发现可以得到循环部分的长度的方法。

余数的奥秘

如果在手动计算的过程中,遇到了某个曾经在余数的位置上出现过的数,那么还有必要继续计算下去吗?显然没有,40除以6得到的第一个余数是4,而这个余数比6小(事实上我们都知道,所有的余数都比除数要小),就需要先将其乘上10再继续计算下去。结果显而易见,一乘上10,就又得到了一开始的40了,这样的话要是继续计算下去,那么结果一定是在重复前面的那些动作,并且永无休止。结论就是,当你遇到了一个余数是之前出现的数的时候,计算就不需要继续下去了,此时根据历史积累的信息,就可以知道循环部分的长度是多少。这也就是我的代码中第一个函数cycle-length的内部辅助函数rec的作用。

rec接收三个参数,第一个参数就是被除数,第二个数是除数,第三个数是这一次计算得到的商在原来的小数部分上是位于哪一个位置上的一个整数。对于每一个出现过的被除数,都将它们记录在了一个内部的哈希表中(为了查找速度的尽可能快),如果一个被除数已经被计算过了,那么就直接取出它所对应的商所在的位置,然后用现在的商所在的位置减掉前者,就知道了循环部分的长度了;如果没有出现过,那么就计算出这一次计算的余数,如果余数为0,意味着除数可以整除被除数,所以没有循环部分可言,长度为0;如果不能整除,那么这个余数的10倍就是新的被除数(因为余数必定小于除数,所以要乘以10直到大于除数为止),位置计数器加上1,重新进行递归。