[算法] 用 Haskell 计算形式幂级数

win_hate

本贴内容来自 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)]，求和代码为：

1. 
   
2. add_series s1 [] = s1
   
3. add_series [] s2 = s2
   
4. add_series (a:s1) (b:s2) = (a+b):(add_series s1 s2)
   

复制代码

[ 本帖最后由 win_hate 于 2009-1-21 14:49 编辑 ]

win_hate

三、积分

Integrate [a(i)] = c:[a(i)/(i+1)]

按 SICP 的做法，我们给出的积分代码只算出 [a(i)/(i+1)]，前面的常数 c 可以在需要的时候加入。

1. 
   
2. import Data.Ratio
   
3. 
   
4. integrate_series::[Ratio Integer]->[Ratio Integer]
   
5. integrate_series s = f s 1 where
   
6.     f (a:s) n = (a/n):(f s (n+1))
   
7.     f [] n = []
   

复制代码

四、例子 exp(x) 和 sin(x), cos (x)

exp(x) 的导数为自身，由于 exp(0)=1，所以 exp(x) 展开式的常数项为 1。仅利用这两个条件，我们就可以形式地计算出 exp(x) 的展开式(在 0 处的幂级数展开)。

在一开始，exp(x) 对应的形式幂级数为 [1, ..]，对它积分，还是 exp(x)，但得到两项 [1,1,...]。
再积分，就能得到第三项，如此下去，就能得到"所有"项。

我们能用各种语言实现上面的描述，但用支持惰性求值的语言来实现显得特别简单。

1. 
   
2. exp_series = 1:(integrate_series exp_series)
   
3. 
   
4. 
   
5. *Main> take 10 exp_series
   
6. [1%1,1%1,1%2,1%6,1%24,1%120,1%720,1%5040,1%40320,1%362880]
   

复制代码

用 maxima 来验证一下：

1. 
   
2. (%i2) taylor (exp(x),x, 0, 10);
   
3.                   2    3    4    5     6      7      8        9        10
   
4.                  x    x    x    x     x      x      x        x        x
   
5. (%o2)/T/ 1 + x + -- + -- + -- + --- + --- + ---- + ----- + ------ + -------
   
6.                  2    6    24   120   720   5040   40320   362880   3628800
   
7.                                                                         + . . .
   

复制代码

下面我们来展开 sin 和 cos. sin 幂级数展开式的常数项为 0， 但它积分后得到的是 -cos. 刚才那套好像不灵了。但注意到 cos 的积分乃是 sin，所以我们可以同时(交互)计算这两个函数的展开。

1. 
   
2. neg_series [] = []
   
3. neg_series (x:s) = (-x):(neg_series s)
   
4. 
   
5. cos_series = 1:(integrate_series (neg_series sin_series))
   
6. sin_series = 0:(integrate_series cos_series)
   
7. 
   
8. neg_series 用于从 cos 的展开得到 -cos 的展开。
   
9. 
   
10. *Main> take 10 sin_series
    
11. [0%1,1%1,0%1,(-1)%6,0%1,1%120,0%1,(-1)%5040,0%1,1%362880]
    
12. *Main> take 11 cos_series
    
13. [1%1,0%1,(-1)%2,0%1,1%24,0%1,(-1)%720,0%1,1%40320,0%1,(-1)%3628800]
    

复制代码

同样用 maxima 验证一下：

1. 
   
2. (%i4) taylor (sin(x), x, 0, 10);
   
3.                           3    5      7       9
   
4.                          x    x      x       x
   
5. (%o4)/T/             x - -- + --- - ---- + ------ + . . .
   
6.                          6    120   5040   362880
   
7. (%i5) taylor (cos(x), x, 0, 10);
   
8.                        2    4    6      8        10
   
9.                       x    x    x      x        x
   
10. (%o5)/T/          1 - -- + -- - --- + ----- - ------- + . . .
    
11.                       2    24   720   40320   3628800
    

复制代码

[ 本帖最后由 win_hate 于 2009-1-21 01:35 编辑 ]

win_hate

五、形式级数的乘积

设有形式幂级数 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 的展开式。于是有如下代码：

1. 
   
2. scale_series k s = map (k*) s
   
3. 
   
4. mul_series (a:s1) s2 = add_series (scale_series a s2)
   
5.                        (0:(mul_series s1 s2))
   

复制代码

scale_series 用一个常数 k 乘上形式幂级数，用来实现 a(0)B(z)。

检验一下效果：sin^2+cos^2 =1

1. 
   
2. *Main> let a = mul_series sin_series sin_series
   
3. *Main> let b = mul_series cos_series cos_series
   
4. *Main> let one = add_series a b
   
5. *Main> take 100 one
   
6. [1%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1,0%1]
   

复制代码

[ 本帖最后由 win_hate 于 2009-1-21 15:04 编辑 ]

win_hate

六、形式幂级数的逆和除法

交换环 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 次项....如此计算出全部。

相应代码为：

1. 
   
2. inv_series (1:s) = inv where
   
3.     inv = 1:(neg_series (mul_series s inv))
   
4. 
   
5. div_series s1 (0:s2) = error "not a unit"
   
6. div_series s1 s2@(a:t) = (mul_series s1' s2') where
   
7.     s1'=scale_series a s1
   
8.     s2'=inv_series  (scale_series a s2)
   
9. 
   
10. *Main> let tan_series = div_series sin_series cos_series
    
11. *Main> take 10 tan_series
    
12. [0%1,1%1,0%1,1%3,0%1,2%15,0%1,17%315,0%1,62%2835]
    

复制代码

maxima 验证：

1. 
   
2. (%i6) taylor (tan(x), x, 0, 10);
   
3.                           3      5       7       9
   
4.                          x    2 x    17 x    62 x
   
5. (%o6)/T/             x + -- + ---- + ----- + ----- + . . .
   
6.                          3     15     315    2835
   

复制代码

win_hate

2003-08-03，CU 网友 loveguohuasai 在 C 版发了个帖子叫 "母牛数量算法"

见这里： http://bbs.chinaunix.net/viewthread.php?tid=130156

这个帖子很火，到目前为止，达到了 25 页共 244 楼。

在我注意到这个贴子的时候，线性递归的方法已经出来了。但用转移矩阵的好像还没出来，幂级数的方法也没有看到(组合上叫生成函数)。

我在 170 楼给出了一个链接，指出这个问题可以用形式幂级数来求解。当时用的是 maple 来求幂级数展开。有了前面的工作后，用 haskell ＋ 幂级数来求母牛数量也是可以的。

我先把题目抄过来：

若一头小母牛，从出生起第四个年头开始每年生一头母牛，按此规律，第n年有多少头母牛？

设第 i 年的母牛数量为 s(i-1)，则有 s(n)=s(n-1)+s(n-3), n>=3, s(0)=s(1)=s(2)=1.

这是个线性递归序列，这类型的序列在数学上已经被研究得很清楚了。密码学中，常用它来生成伪随机序列。多数流密码中，都包含一个这样的序列发生器。

求 s(n) 的方法有以下几种：

• 迭代/线性递归;
  
• 转移矩阵。这个可以参考 SICP 中文版第 31 页，或下面连接中的第 22，23 楼
  [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 去乘上面的级数，得到

(1-x-x^3) (s(0)+s(1)x+s(2)x^2+ ... ) =

s(0) + s(1)x + s(2)x^2 + s(3)x^3 + s(4)x^4 + ...
     - s(0)x - s(1)x^2 - s(2)x^3 - s(3)x^4 + ...
                       - s(0)x^3 - s(1)x^4 + ...

直观上看，就是把原序列 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 .....，求其逆就得到母牛幂级数。

1. 
   
2. *Main> let cows = inv_series [1,-1,0,-1]
   
3. *Main> take 20 cows
   
4. [1,1,1,2,3,4,6,9,13,19,28,41,60,88,129,189,277,406,595,872]
   
5. 
   
6. *Main> cows!!99
   
7. 16637075746565861
   

复制代码

[ 本帖最后由 win_hate 于 2009-1-21 01:32 编辑 ]

MMMIX

原帖由 win_hate 于 2009-1-21 01:18 发表
比如级数 1+x+x^2+... 就可以表示为 [1..]

[1..] 得到的就是 1, 2, 3, ...

要得到一列 1 可以用 repeat 1

[ 本帖最后由 MMMIX 于 2009-1-21 11:27 编辑 ]

MMMIX

原帖由 win_hate 于 2009-1-21 01:24 发表

scale_series k [] = []
scale_series k (a:s) = (k*a):(scale_series k s)

这个用 map 就行了，也即
scale_series k xs = map (k*) xs

如果你愿意，也可以写成
scale_series k = map (k*)

[ 本帖最后由 MMMIX 于 2009-1-21 11:49 编辑 ]

win_hate

原帖由 MMMIX 于 2009-1-21 11:26 发表

[1..] 得到的就是 1, 2, 3, ...

要得到一列 1 可以用 repeat 1

已经修改为 [1,1..] 了，谢谢。

原帖由 MMMIX 于 2009-1-21 11:42 发表

这个用 map 就行了，也即
scale_series k xs = map (k*) xs

如果你愿意，也可以写成
scale_series k = map (k*)

这个写法好。

[ 本帖最后由 win_hate 于 2009-1-21 14:54 编辑 ]

MMMIX

原帖由 win_hate 于 2009-1-21 14:51 发表


这个写法好。

这个称之为 point free style

win_hate

原帖由 MMMIX 于 2009-1-21 15:37 发表

这个称之为 point free style

point free 指的是 map (k*) 还是 (k*) 呢？

再给两个母牛的解法，其中前一个是模仿 MMMIX 处理 fib 序列的代码写的。

1. 
   
2. cows = 1:1:1:[x+y|x<-cows| y<- drop 2 cows]  -- 这个要用 ghci -XParallelListComp
   

复制代码

另一个：

1. 
   
2. cows = 1:1:1:(add_series cows (drop 2 cows))
   

复制代码

1. 
   
2. *Main> cows!!99
   
3. 16637075746565861
   

复制代码