[算法] 用 Haskell 计算形式幂级数
本贴内容来自 SICP 习题 3.59 - 3.62. 在本帖中，我将使用 Haskell，因为它的惰性求值是内建的，而且我可以顺便熟悉一下 Haskell. 对 Haskell 了解有限，如果有些地方写得不够“Haskell”，请帮忙指出。谢谢。
一、形式幂级数的表示
a(0)+a(1)x+ ... + a(n)x^n + ...
在 Haskell 中，我们可以用一个无穷列表来表示它: [a(0), a(1), .... ]。当然，我们要给出一种规则，从有限项生成无限项。
比如级数 1+x+x^2+... 就可以表示为 [1,1..]
而 1+2x+3x^2+4x^3 +... 就可以表示为 [1,2..]
如果形式幂级数的系数只有有限个是非零的，也可以用一个有穷列表来表示。
二、求和
[a(i)] + [b(i)] = [a(i)+b(i)]，求和代码为：
- add_series s1 [] = s1
- add_series [] s2 = s2
- add_series (a:s1) (b:s2) = (a+b):(add_series s1 s2)
复制代码
[ 本帖最后由 win_hate 于 2009-1-21 14:49 编辑 ]
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（9）
point free 指的是 map (k*) 还是 (k*) 呢？
再给两个母牛的解法，其中前一个是模仿 MMMIX 处理 fib 序列的代码写的。
复制代码
另一个：
复制代码
复制代码
这个称之为 point free style
已经修改为 [1,1..] 了，谢谢。
这个写法好。
[ 本帖最后由 win_hate 于 2009-1-21 14:54 编辑 ]
这个用 map 就行了，也即
scale_series k xs = map (k*) xs
如果你愿意，也可以写成
scale_series k = map (k*)
[ 本帖最后由 MMMIX 于 2009-1-21 11:49 编辑 ]
[1..] 得到的就是 1, 2, 3, ...
要得到一列 1 可以用 repeat 1
[ 本帖最后由 MMMIX 于 2009-1-21 11:27 编辑 ]
2003-08-03，CU 网友 loveguohuasai 在 C 版发了个帖子叫 "母牛数量算法"
见这里： http://bbs.chinaunix.net/viewthread.php?tid=130156
这个帖子很火，到目前为止，达到了 25 页共 244 楼。
在我注意到这个贴子的时候，线性递归的方法已经出来了。但用转移矩阵的好像还没出来，幂级数的方法也没有看到(组合上叫生成函数)。
我在 170 楼给出了一个链接，指出这个问题可以用形式幂级数来求解。当时用的是 maple 来求幂级数展开。有了前面的工作后，用 haskell ＋ 幂级数来求母牛数量也是可以的。
我先把题目抄过来：
设第 i 年的母牛数量为 s(i-1)，则有 s(n)=s(n-1)+s(n-3), n>=3, s(0)=s(1)=s(2)=1.
这是个线性递归序列，这类型的序列在数学上已经被研究得很清楚了。密码学中，常用它来生成伪随机序列。多数流密码中，都包含一个这样的序列发生器。
求 s(n) 的方法有以下几种：
[http://bbs.chinaunix.net/thread-1263002-3-1.html]
这个方法见：
http://blog.chinaunix.net/u/20/showart_215141.html
这个方法见下面以及原贴第 170 楼。
各年母牛数量为 s(0), s(1), s(2), .... 可以对应为一个形式幂级数:
s(0)+s(1)x+s(2)x^2+ ...
注意到 s(n)-s(n-1)-s(n-3)=0 的条件，我们用 1-x-x^3 去乘上面的级数，得到
直观上看，就是把原序列 S， 右移一列的序列 Sx 和 右移 3列的 Sx^3 进行线性组合。容易看出，从第四列起，全都是 0，因为有 s(n)-s(n-1)-s(n-2)=0。
不难验证，第二，三列也是 0，只有第一列为 1。所以我们有：
(1-x-x^3) (s(0)+s(1)x+s(2)x^2+ ... ) = 1
或 s(0)+s(1)x + s(2) x^2+ ... = 1/(1-x-x^3)
把 1-x-x^3 看成幂级数 1-x-0x^2-x^3+ 0x^4 +0x^5 .....，求其逆就得到母牛幂级数。
复制代码
[ 本帖最后由 win_hate 于 2009-1-21 01:32 编辑 ]
六、形式幂级数的逆和除法
交换环 A 上的形式幂级数构成一个环 A[[x]]，有逆的充分必要条件是常数项可逆。由于我们这里只处理实数，所以有逆相当于常数项非零。说明 a(0)+a(1)x+a(2)x^2 + ... 可逆的常见方法是设其逆为 b(0)+b(1)x+b(2)x^2+....然后声称所有的 b(i) 可以从
(a(0)+a(1)x+...) (b(0)+b(1)x+...) = 1
求出。
SICP 提供了另一种描述，本质上是一样的。令 S = 1+ S' 是一个常数项为 1 的形式幂级数，其逆设为 X。有 S X = 1 => (1+S')X=1 => X=1-S'X
X 的常数项为 1， 从 X=1-S'X 可以计算出 X 的一次项。代入右边，可以计算出 2 次项....如此计算出全部。
相应代码为：
复制代码
maxima 验证：
复制代码
五、形式级数的乘积
设有形式幂级数 A(z)=a(0)+A'(z)z, B(z), 则
A(z)B(z) = a(0) B(z)+ A'(z)B(z)z
从 a(0) B(z) 中已经可以得到 A(z)B(z) 的第一项，A'(z)B(z)z 等于把 A'(z)B(z) 往右推了一步，若 A'(z)B(z) 的展开式为 s，则 0:s 是 A'(z)B(z)z 的展开式。于是有如下代码：
复制代码
scale_series 用一个常数 k 乘上形式幂级数，用来实现 a(0)B(z)。
检验一下效果：sin^2+cos^2 =1
复制代码
[ 本帖最后由 win_hate 于 2009-1-21 15:04 编辑 ]
三、积分
Integrate [a(i)] = c:[a(i)/(i+1)]
按 SICP 的做法，我们给出的积分代码只算出 [a(i)/(i+1)]，前面的常数 c 可以在需要的时候加入。
复制代码
四、例子 exp(x) 和 sin(x), cos (x)
exp(x) 的导数为自身，由于 exp(0)=1，所以 exp(x) 展开式的常数项为 1。仅利用这两个条件，我们就可以形式地计算出 exp(x) 的展开式(在 0 处的幂级数展开)。
在一开始，exp(x) 对应的形式幂级数为 [1, ..]，对它积分，还是 exp(x)，但得到两项 [1,1,...]。
再积分，就能得到第三项，如此下去，就能得到"所有"项。
我们能用各种语言实现上面的描述，但用支持惰性求值的语言来实现显得特别简单。
复制代码
用 maxima 来验证一下：
复制代码
下面我们来展开 sin 和 cos. sin 幂级数展开式的常数项为 0， 但它积分后得到的是 -cos. 刚才那套好像不灵了。但注意到 cos 的积分乃是 sin，所以我们可以同时(交互)计算这两个函数的展开。
复制代码
同样用 maxima 验证一下：
复制代码
[ 本帖最后由 win_hate 于 2009-1-21 01:35 编辑 ]