CommonLisp函数调用的CPS变换

之前在新浪博客上已经发表了一篇同样主题的文章了,不过心血来潮想重写一份发到这边的博客上来,顺便看看自己能不能够在自己重写这篇的东西的过程中再发掘出一些有价值的东西,充实一下内容。

什么是CPS

首先双手奉上来自维基百科的Continuation-passing-style条目。CPS一般译作续延传递风格,这里最不好理解的词语莫过于续延了。维基百科上关于续延也有详细的介绍,不过我也可以简单地介绍一下。续延是一种表示接下来所要进行的计算的对象,也就是你在计算了这一步之后需要继续做的事情的一个描述。举个例子,对于Common Lisp代码(+ (+ 1 2) 3)而言,最左边的参数首先会被计算,也就是Lisp实现先计算的是(+ 1 2)这个表达式,于是可以得到返回值为3。那么计算完了最左边的参数之后,Lisp实现接下来要做什么呢?显然,接下来会计算作为直接值的3得到3本身,然后再将第一次计算得到的3和第二次求值得到的3进行相加。而上面的这两个步骤,其实可以用Lisp中的对象进行描述,只需要用一个匿名函数就可以了,写成

(lambda (v) (+ v 3))。

那么我们就可以思考,在这一个匿名函数被计算之后,再一次的接下来的计算又是什么呢?对于在REPL中输入的代码而言,Lisp实现会把表达式最终计算得到的返回值打印在终端上,所以我们可以认为上面的匿名函数,也就是一个计算过程,所对应的续延是

(lambda (u) (print u))。

那么原来的表达式(+ (+ 1 2) 3)就可以写成下面这样的风格的代码:

(+& 1 2 (lambda (v)
          (+& v 3 (lambda (u)
                    (print u)))))

你或许会感到奇怪:这里的+&是个什么东西啊?+&是一个和加法类似的函数,只是它除了会把前两个参数相加之外,还会将这个参数传递给作为第三个参数传入的匿名函数,并进行函数调用。按照维基百科上的说明,CPS代码有两个规范:

  1. 每一个函数都要接受一个续延作为额外的参数;
  2. 每一个函数调用中的参数都只能是一个变量或者λ表达式。

可以说,+&这样的函数是为了将代码变换成CPS才创造出来的。当然了,除了+函数需要进行变换之外,其它的普通函数(regular function)也是需要进行变换的。例如对代码(\ (+ 1 2) 3)进行变换的话,就需要为\函数提供相应的变体,结果应该是:

(+& 1 2 (lambda (v)
          (*& v 3 (lambda (u)
                    (print u)))))。

上面的代码中有个可以稍微琢磨的地方:为什么在两个不同的匿名函数的参数列表中偏偏要使用不同的参数名呢?事实上,在上面的两个变换出来的结果中,完全可以统一使用v或者u来作为参数名,那是因为参数v不需要在第二层嵌套的匿名函数中被使用。可是有时候一开始就进行求值的表达式不一定是马上就被使用的,例如代码(- (+ 1 2) (\ 3 4))*的变换结果就是

(+& 1 2 (lambda (v)
          (*& 3 4 (lambda (u)
                    (-& v u)))))。

可以说,当参数列表表达式中不止一个复杂表达式时,就需要至少两个不重名的参数了。讲了这么多手动变换的事情,现在可以讨论最重要的问题了:

如何对Lisp代码进行CPS变换

如果我们可以对一类简单的代码进行CPS变换了,那么如何在这个基础上对复杂的代码进行CPS变换呢?例如我们有一段简单的代码(op3 (op2 (op1 a b) c) d),并且我们已经有了一个通用的功能F可以对其中的嵌套的表达式(op2 (op1 a b) c)进行变换,那么要如何将其和外层的对op3函数的调用结合在一起呢?显然,我们可以先把(op2 (op1 a b) c)换掉,换成一个普通的符号v,以表示一个变量。那么接下来我们的处理对象就是一个稍微简单的表达式(op3 v c)了,对于这种参数列表均由变量所构成的表达式,要将其变换成CPS代码只需要修改其函数为加&的形式,并且添上一个额外的参数——续延就可以了。至于这个续延应该是由调用者提供,而又可以使用默认的(lambda (XXX) (print XXX))的。所以,上面的代码最终会变成

(经过F处理的代码 (lambda (v)
                (op3& v c (lambda (XXX)
                            (print XXX)))))。

对内层的表达式应用F处理完之后还要再计算外层表达式的续延,而计算续延正是函数F所做的事情,所以实际上对一个复杂的表达式做CPS变换,是一个递归的过程。用两段代码进行变换可以有更直观的体验:假设我们有一个正确实现的变换函数F,它接受一个代码和一个续延,并产生出对代码进行CPS变换所得到的结果。那么如果应用在代码(+ (+ 1 1) 3)上,就是

(F (+ (+ 1 1) 3) (lambda (v) (print v)))。

显然,这段代码应该和下面的这段代码是等价的:

(F (+ 1 1) (lambda (u)
             (+& u 3 (lambda (v)
                       (print v)))))。

并且上面这段代码的续延参数部分,和下面的这段代码是等价的:

(F (+ u 3) (lambda (v) (print v)))。

是不是就是说,对于一个复杂的表达式expr0,如果想要在提供一个额外续延cont0的情况下,将其正确变换成CPS形式,那么可以先提取出它的第一个被求值的子表达式expr1,将原来子表达式出现的地方替换为一个变量,那么原本的表达式就变成了expr1对应的续延cont1。此时先将cont1和cont0传入函数F进行CPS变换,得到的结果cont将会成为和expr1一起传入函数F进行变换,并最终得到对expr0和cont0进行变换的正确结果呢?显然这是对的,因为cont0紧随expr0求值,而expr0可以拆分为两部分:expr1以及紧随其求值的cont1。那么显然cont0紧随cont1求值,所以cont0是cont1的续延,同样需要进行CPS变换,然后得到新的代码作为expr1的续延继续进行变换,那么最终就会得到代码,它反映的求值顺序正是expr1 -> cont1 -> expr0。

那么重点就在于怎么把原本的expr0给拆分成expr1和cont1,以及具体结合代码和续延的技术细节了。假设实现拆分功能的函数叫做split-funcall,那么它需要在什么情况下做什么事情呢?一般而言,split-funcall处理的是有复杂的子表达式的代码,那么它需要能够在整个表达式expr0中找出尽量左的一个复杂子表达式,也就是expr1。如果找到这样的expr1了,那么就需要同时找出对应的续延cont1并将两者都返回;如果没有这样的expr1,那么直接返回逻辑假和整个expr0就可以了。说到需要返回多个值的情况,正是Common Lisp的values函数派上大用场的时候了。

对于含有复杂子表达式的情况,split-funcall函数会返回提取到的子表达式和对应的续延,那么这里的续延应该用什么样的形式来表示呢?为了方便使用从split-funcall内部所gensym出来的符号,所以直接返回一个经过构造的lambda表达式就可以了,但是在递归调用函数F时,不能直接对其进行变换处理,因为lambda表达式不是一个会对参数求值的操作符,所以真正需要变换的是lambda表达式的函数体。因此必须先拆开lambda表达式,对其函数体进行变换并再次构造成lambda表达式,以传递给函数F使用。

最后,贴一下自己放在Gist上的代码:

连分数与ProjectEuler的第66题

可耻的蛮力法

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和闭包中,我更喜欢这种用闭包来实现的风格;-)

PE的第26题

今天上百度贴吧的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,重新进行递归。