斐波那契数列高效算法之快速幂法、快速倍增法与扩域法

目前网络上介绍斐波那契数列算法的文章数不胜数, 但是高效算法主要围绕矩阵快速幂法展开. 但实际上还有两种算法: 快速倍增法和扩域法, 均能达到 \(O(\log n)\) 的时间复杂度, 在一定程度的优化后甚至能做到更小的常数, 从而拥有更高的效率. 下面介绍三种算法, 同时给出必要的原理与简要推导.
快速幂法

既然要讲快速算法, 首先还是要介绍快速幂法, 因为其是快速倍增法的基础.
快速幂法基于以下递推式:

\[\begin{bmatrix}F(n)\\F(n+1)\end{bmatrix}=\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^{n}\begin{bmatrix}F(0)\\F(1)\end{bmatrix}\]
或可作

\[\begin{bmatrix}F(n-1) & F(n)\\F(n) & F(n+1)\\\end{bmatrix}=\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^{n}\]
这个式子原理比较简单, 网上也有很多推导过程, 故此处不再列出.
基于第二个递推式, 可构造一个无符号整数结构体, 为其实现乘法与快速幂运算, 即可通过求 \(n\) 次方得到 \(F(n)\).

1. use std::ops::Mul;
3. #[derive(Copy,Debug, Clone)]
4. struct Mat2x2(u64, u64, u64, u64);
6. impl Mul for Mat2x2 {
7.     type Output = Self;
8.     fn mul(self, rhs: Self) -> Self::Output {
9.         Self(
10.             self.0 * rhs.0 + self.1 * rhs.2,
11.             self.0 * rhs.1 + self.1 * rhs.3,
12.             self.2 * rhs.0 + self.3 * rhs.2,
13.             self.2 * rhs.1 + self.3 * rhs.3,
14.         )
15.     }
16. }
17. impl Mat2x2 {
18.     fn pow(mut self, mut exp: u64) -> Self {
19.         let mut ans = Self(1, 0, 0, 1);
20.         while exp > 0 {
21.             if exp & 1 == 1 {
22.                 ans = ans * self;
23.             }
24.             self = self * self;
25.             exp >>= 1;
26.         }
27.         ans
28.     }
29. }
30. fn fib(n: u64) -> u64 {
31.     let ori = Mat2x2(0, 1, 1, 1);
32.     ori.pow(n).1
33. }

复制代码

由快速幂原理可得, 时间复杂度取决于 pow 函数中指数的二进制位数, 故时间复杂度为 \(O(\log n)\); 使用的空间为常数, 故空间复杂度为 \(O(1)\).
快速倍增法

通过矩阵幂可推导出以下公式

\[\begin{cases}F(0) &= 0\\F(1) &= 1\\F(2n) &= F(n)(F(n) + 2 F(n-1))\text{ if }n>1 \\F(2n+1) &= F(n)^2 + F(n+1)^2\text{ if } n > 2\end{cases}\]
其可通过快速幂法推导得到

\[\begin{align*}\begin{bmatrix}F(2n-1) & F(2n)\\F(2n) & F(2n+1)\\\end{bmatrix}&=\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^{2n}\\&=\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^{n}\begin{bmatrix}0 & 1\\1 & 1\end{bmatrix}^{n}\\&=\begin{bmatrix}F(n-1) & F(n)\\F(n) & F(n+1)\end{bmatrix}\begin{bmatrix}F(n-1) & F(n)\\F(n) & F(n+1)\end{bmatrix}\\\end{align*}\]
通过矩阵乘法公式计算 \(F(2n)\) 和 \(F(2n+1)\) 项即可得到快速倍增法公式.
基于快速倍增法公式结合动态规划法可以得到如下代码:

1. fn fib_double(n: u64) -> u64 {
2.     let mut map = vec![1; n as usize + 1];
3.     map[0] = 0;
5.     fn inner_fib(n: u64, map: &mut Vec<u64>) -> u64 {
6.         if n == 1 || n == 2 {
7.             return 1;
8.         }
9.         match map[n as usize] {
10.             1 => {
11.                 let val = if n % 2 == 0 {
12.                     inner_fib(n / 2, map) * (inner_fib(n / 2, map) + 2 * inner_fib(n / 2 - 1, map))
13.                 } else {
14.                     inner_fib(n / 2, map).pow(2) + inner_fib(n / 2 + 1, map).pow(2)
15.                 };
16.                 map[n as usize] = val;
17.                 val
18.             }
19.             val => val,
20.         }
21.     }
22.     inner_fib(n, &mut map)
23. }

复制代码

理论上快速倍增法的时间复杂度仍为 \(O(\log n)\), 但由于递归的存在, 实际的速度不如矩阵快速幂法. 虽然可以用栈消除递归, 但这会引入额外的消耗, 所以不再进一步探索. 同时该方法的空间复杂度为 \(O(n)\), 如果使用哈希表减少额外空间的使用则空间复杂度可以减少到 \(O(\log n)\)(时间复杂度常数会上升), 这个方法的优势在于能够得到中间值, 在特定场合可能有作用. 但大部分场合效果不如快速幂与扩域法.
扩域法

终于来到了这篇文章的重点, 扩域法.
首先我们看一下斐波那契数列的通项公式:

\[F(n) = \frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n]\]
利用这个公式理论结合浮点数快速幂能够在 \(O(\log n)\) 时间内求解, 但是浮点数运算存在精度问题, 当 \(n\) 较大时很快就会精度不足.
但注意运算过程中只存在 \(\sqrt{5}\) 这唯一一个无理数, 所有数字均可表示为 \(a+b\sqrt{5}, a,b\in Q\) 形式. 如果能通过某个方法, 在运算过程中消除这个无理数, 或者用某个数据结构代替这个无理数, 就有可能避免浮点运算.
在中学阶段, 学习到复数时, 我们通过引入 \(i\) 和数对 \((a,b)\) 让复数 \(a+bi\) 运算能够用实数计算. 仿照复数的这种构造方式, 令

\[a+b\sqrt{5} = (a,b)\]
模仿虚数的实部虚部的称呼, 本文将 \(a\) 称为有理部, \(b\) 为无理部系数.
定义取有理部和无理部系数函数:

\[\begin{align*}Ra[(a,b)] &= a\\Ir[(a, b)] &= b\end{align*}\]
通过简单推算得到其四则运算的公式

\[\begin{align*}(a,b)\pm (c,d) &=(a\pm c, b\pm d) \\(a,b) \times (c,d) &= (ac + 5bd, ad+bc) \\\frac{(a, b)}{(c,d)} &= (\frac{ac-5bd}{c^2-5d^2},\frac{bc-ad}{c^2-5d^2})\end{align*}\]
可以发现计算结果仍满足 \((a,b)\) 结构, 在数学上这构成了一个 \(Q(\sqrt{5})\) 数域, 所以这个方法称为扩域法.
使用这个方法可以将斐波那契数列公式改写为

\[F(n) = \frac{1}{(0,1)}[(\frac{(1,1)}{2})^n - (\frac{(1,-1)}{2})^n]\]
但这个公式目前仍不能完全消除浮点数的存在, 因为在其中还存在以下几个非整数项

\[\begin{align*}&\frac{1}{(0,1)}\\&\frac{(1,1)}{2}\\&\frac{(1,-1)}{2}\\\end{align*}\]
但是可以发现, 第一项 \(1/(0,1)\) 乘以后续结果的作用就是将后续结果的无理部系数提取出来, 所以可以通过以下公式避免第一项的计算

\[F(n) = Ir[(\frac{(1,1)}{2})^n - (\frac{(1,-1)}{2})^n]\]
对于 \((1,1)/2\) 等导致的浮点数, 我们可以将除法放到最后计算, 得到公式

\[F(n) = \frac {Ir[(1,1)^n-(1,-1)^n]}{2^n}\]
到此我们成功消除了所有浮点数的存在, 可以得到初步的代码

1. use std::ops::{Div, Mul, Sub};
3. #[derive(Clone, Copy)]
4. pub struct QSqrt5(i64, i64);
6. impl Mul for QSqrt5 {
7.     type Output = Self;
8.     fn mul(self, rhs: Self) -> Self::Output {
9.         Self(
10.             self.0 * rhs.0 + 5 * self.1 * rhs.1,
11.             self.0 * rhs.1 + self.1 * rhs.0,
12.         )
13.     }
14. }
16. impl Sub for QSqrt5 {
17.     type Output = Self;
18.     fn sub(self, rhs: Self) -> Self::Output {
19.         Self(self.0 - rhs.0, self.1 - rhs.1)
20.     }
21. }
23. impl Div<i64> for QSqrt5 {
24.     type Output = Self;
25.     fn div(self, rhs: i64) -> Self::Output {
26.         Self(self.0 / rhs, self.1 / rhs)
27.     }
28. }
30. impl QSqrt5 {
31.     fn pow(mut self, mut exp: u64) -> Self {
32.         let mut ans = QSqrt5(1, 0);
33.         while exp > 0 {
34.             if exp & 1 == 1 {
35.                 ans = ans * self;
36.             }
37.             self = self * self;
38.             exp >>= 1;
39.         }
40.         ans
41.     }
42. }
44. fn fib(n: u64) -> u64 {
45.     (QSqrt5(1, 1).pow(n) - QSqrt5(1, -1).pow(n)).1 as u64 >> n
46. }
48. fn main() {
49.     println!("{}", fib(7));
50. }

复制代码

以上已经能够得到斐波那契数列了, 且时间复杂度为两次快速幂的指数二进制位数, 同样为 \(O(\log n)\), 使用了常数空间, 故空间复杂度为 \(O(1)\).
扩域法优化

通过之前的优化我们成功消除了公式法中所有无理数和小数, 但存在两个问题:

• 能否消除其中的负数计算, 只保留无符号整数
  
• 最后整体除 \(2^n\) 说明计算中存在大量因子没有消除, 中间项很大可能导致溢出 (虽然在模 \(2^n\) 数域中这基本不成问题)
  

通过推导我们发现

\[\begin{align*}Ra[(1,1)^n] &= Ra[(1, -1)^n]\\Ir[(1,1)^n] &= -Ir[(1, -1)^n]\end{align*}\]
所以只需要计算其中任意一个就能得到另一个, 这样就可以不必计算无理部系数为负数的项, 得到新公式

\[\begin{align*}F(n) &= \frac {Ir[(1,1)^n-(1,-1)^n]}{2^n}\\&= \frac {Ir[(1,1)^n]-Ir[(1,-1)^n]}{2^n}\\&= \frac {Ir[(1,1)^n]-(-Ir[(1,1)^n])}{2^n}\\&= \frac {2Ir[(1,1)^n]}{2^n}\end{align*}\]
这样就可以直接用无符号整数计算.
在求快速幂时, 计算中间值的同时将 \(2\) 及时约去可以避免溢出, 注意到初始项的 \((1,1)\) 不能被 \(2\) 整除, 所以我们将累乘结果 ans 的初始值设为 \((2, 0)\) 而非 \((1,0)\), 相当于 \(2Ir[(1,1)^n]/2^n\) 分子上的 \(2\), 不将其提前约去来保证每一项都能被 \(2\) 整除.
至此我们得到最终版的扩域法代码

1. use std::ops::{Div, Mul};
3. #[derive(Clone, Copy, Debug)]
4. struct QSqrt5(u64, u64);
6. impl Mul for QSqrt5 {
7.     type Output = Self;
8.     fn mul(self, rhs: Self) -> Self::Output {
9.         Self(
10.             self.0 * rhs.0 + 5 * self.1 * rhs.1,
11.             self.0 * rhs.1 + self.1 * rhs.0,
12.         )
13.     }
14. }
16. impl Div<u64> for QSqrt5 {
17.     type Output = Self;
18.     fn div(self, rhs: u64) -> Self::Output {
19.         Self(self.0 / rhs, self.1 / rhs)
20.     }
21. }
23. impl QSqrt5 {
24.     fn pow_div_2(mut self, mut exp: u64) -> Self {
25.         let mut ans = QSqrt5(2, 0);
26.         while exp != 0 {
27.             if exp & 1 == 1 {
28.                 ans = ans * self / 2;
29.             }
30.             self = self * self / 2;
31.             exp >>= 1;
32.         }
33.         ans
34.     }
35. }
37. fn fib(n: u64) -> u64 {
38.     QSqrt5(1, 1).pow_div_2(n).1
39. }
41. fn main() {
42.     println!("{}", fib(63));
43. }

复制代码

这个方法的时间复杂度仍为 \(O(\log n)\), 空间复杂度仍为 \(O(1)\), 但常数较小.
总结

以上三个方法都能实现时间复杂度在 \(O(\log n)\)的条件下求第 \(n\) 项斐波那契数. 其中快速倍增法因为递归调用性能较差, 通过充分的优化, 快速幂和扩域法的耗时其实非常相似, 在 rust 的 release 模式测试中扩域法优化后也仅比快速幂法快了不到 \(10\%\), 所以死磕性能意义已经不大.
虽然性能提升不高, 但本文的扩域算法是个很有意思的思路, 或许能给苦思算法的读者一点启发. 同时对于那些设计数学计算库的读者给出了一个结论:

• 如果需要生成 \(0\) ~ \(n\) 每一项的斐波那契数, 直接使用递推法 \(F(n) = F(n-1) + F(n-2)\) 就是最快的方案
  
• 如果只生成第 \(n\) 项斐波那契数, 则使用快速幂和扩域法都可以, 如果非常关注性能则可以选择扩域法
  

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！