Haskell: 基于方阵快速幂，求解斐波那契数列项

Haskell: 基于方阵快速幂，求解斐波那契数列项

原理

• 斐波那契：$F(n+1)=F(n)+F(n-1)$, $F(-1)=0,F(0)=1,F(1)=1,F(2)=2,...$

• 使用矩阵的幂求解斐波那契：
  $$\left[\begin{array}{c}F(n+1)\\ F(n)\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\ast \left[\begin{array}{c}F(n)\\ F(n-1)\end{array}\right]$$
  $$\left[\begin{array}{c}F(n+1)\\ F(n)\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n\ast \left[\begin{array}{c}F(1)\\ F(0)\end{array}\right]={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n\ast \left[\begin{array}{c}1\\ 1\end{array}\right]$$

• 矩阵的快速幂方法：基于公式$A^{2m}=A^m\ast A^m$和分治法，可以通过$O(logn)$的对数时间计算$A^n$

• 可以通过快速幂来快速求解 ${\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^n$，从而在对数时间内求解斐波那契的第n项。

代码实现

实现方阵的乘法

my_product :: Num a => [[a]] -> [[a]] -> [[a]]
my_product a b = [[sum [a !! i !! k * b !! k !! j | k <- [0..dim - 1]] | j <- [0..dim - 1]] | i <- [0..dim - 1]]
                 where dim = 2

定义中dim=2是方阵的维数。因为要用到的是2*2方阵的乘法，所以dim=2.

然后，因为此处是自乘，写一个自乘函数比较方便：

self_product :: Num a => [[a]] -> [[a]]
self_product = \a -> my_product a a

这样做，比直接定义自乘函数泛用性要强，而且并不会使计算变慢。

Haskell有惰性求值的特性，在上述函数中，参数a被函数体中的my_product第一次引用时求值，但是只会求一次值，第二次引用时不会重新求一遍。

实现分治法快速求幂

以下是第一版代码：

fast_exp :: (Num a, Integral int) => [[a]] -> int -> [[a]]
fast_exp _ 0 = [[1, 0], [0, 1]]
fast_exp a num = if num `mod` 2 == 0 
                 then self_product $ fast_exp a $ halfnum
                 else my_product a $ self_product $ fast_exp a $ halfnum
                 where halfnum = floor $ fromIntegral num / 2

• 脑子里时刻想明白求值顺序,通过检查现在写出来的
• 这里主要练习了一下$的使用。$的使用其实就是改变求值顺序。
  • 编程时，应该先想明白求值顺序，然后灵活使用$，而不是用$等价替换括号。
  • 编程时，永远不要想着用$完全代替括号，但是在括号嵌套很深的时候，可以通过使用$减少括号的数量。
• 注意halfnum的实现。为了类型正确，中间需要一步fromIntegral
• 以上提到的所有函数，如果用到，一定通过:help来查看函数的具体函数签名，来保证写出来的程序正确。

在此基础上，试着不用if-else，而是使用模式匹配写出来：

fast_exp :: (Num a, Integral int) => [[a]] -> int -> [[a]]
fast_exp a num | num == 0 = [[1, 0], [0, 1]]
               | num `mod` 2 == 0 = self_product $ fast_exp a $ halfnum
               | otherwise = my_product a $ self_product $ fast_exp a $ halfnum
               where halfnum = floor $ fromIntegral num / 2

这样写更有Haskell风范。

实现斐波那契数列项的求解

fibo :: Integral int => int -> int
fibo num = fast_exp fibo_base num !! 0 !! 0 where fibo_base = [[1, 1], [1, 0]]

通过快速幂求解矩阵，然后再取出对应位置的元素，即为数列的项值。

完成

在终端里面运行ghci打开Haskell交互平台，用:l加载脚本，然后可以调用函数求解。

*Main> fibo 100
573147844013817084101
(0.00 secs, 159,264 bytes)

测试时间

用:set +s显示函数执行状态，从而可以测试指令执行时间。
在此基础上，用下面这个函数封装一下，使得计算过程不变，但是不会被显示数字的时间干扰。

my_zero :: Integral int => int -> int
my_zero a = if a == 0 then 0 else 1

通过下面的执行结果可见，显示一长串数字确实影响对测试结果的计算。用my_zero封装一下可以排除显示数字所需时间的干扰。

*Main> fibo 2000
6835...26(many numbers)
(0.05 secs, 569,920 bytes)

*Main> my_zero $ fibo 2000
1
(0.02 secs, 208,080 bytes)

*Main> my_zero $ fibo 200000
1
(0.00 secs, 509,608 bytes)

那么我们测试一下速度（从1e7开始倍增）：

*Main> my_zero $ fibo 10000000
1
(0.16 secs, 12,081,960 bytes)
*Main> my_zero $ fibo 20000000
1
(0.33 secs, 23,806,616 bytes)
*Main> my_zero $ fibo 40000000
1
(0.73 secs, 47,246,816 bytes)
*Main> my_zero $ fibo 80000000
1
(1.61 secs, 94,117,912 bytes)
*Main> my_zero $ fibo 160000000
1
(3.53 secs, 187,849,912 bytes)

果然是一个对数时间复杂度。
（在数字达到1e9之后，可能因为内存不够用了而时间上涨。但整体没什么问题。）

完整代码

-- 方阵快速幂.hs

my_product :: Num a => [[a]] -> [[a]] -> [[a]]
my_product a b = [[sum [a !! i !! k * b !! k !! j | k <- [0..dim - 1]] | j <- [0..dim - 1]] | i <- [0..dim - 1]]
                 where dim = 2

{-
self_product :: Num a => [[a]] -> [[a]]
self_product a = [[sum [a !! i !! k * a !! k !! j | k <- [0..1]] | j <- [0..1]] | i <- [0..1]]
-}

self_product :: Num a => [[a]] -> [[a]]
self_product = \a -> my_product a a

{-
fast_exp :: (Num a, Integral int) => [[a]] -> int -> [[a]]
fast_exp _ 0 = [[1, 0], [0, 1]]
fast_exp a num = if num `mod` 2 == 0 
                 then self_product $ fast_exp a $ halfnum
                 else my_product a $ self_product $ fast_exp a $ halfnum
                 where halfnum = floor $ fromIntegral num / 2
                 -}

fast_exp :: (Num a, Integral int) => [[a]] -> int -> [[a]]
fast_exp a num | num == 0 = [[1, 0], [0, 1]]
               | num `mod` 2 == 0 = self_product $ fast_exp a $ halfnum
               | otherwise = my_product a $ self_product $ fast_exp a $ halfnum
               where halfnum = floor $ fromIntegral num / 2

fibo :: Integral int => int -> int
fibo num = fast_exp fibo_base num !! 0 !! 0 where fibo_base = [[1, 1], [1, 0]]

my_zero :: Integral int => int -> int
my_zero a = if a == 0 then 0 else 1