P30 - 矩阵快速幂计算斐波那契数
Fibonacci numbers with matrix exponentiation
官方模块:
Problems.P30核心函数:fibonacci'
核心公式
利用转移矩阵
问题变成整数矩阵的快速幂。
函数签名
1fibonacci' :: Integral a => a -> a
实现
方法一:矩阵快速幂(O(log n))
1type Matrix a = ((a, a), (a, a))
2
3fibonacci' :: Integral a => a -> a
4fibonacci' n
5 | n < 1 = error "fibonacci': index must be positive"
6 | otherwise = fst (apply (power q (n - 1)) (1, 0))
7 where q = ((1, 1), (1, 0))
8
9apply :: Num a => Matrix a -> (a, a) -> (a, a)
10apply ((a,b),(c,d)) (x,y) = (a*x + b*y, c*x + d*y)
11
12multiply :: Num a => Matrix a -> Matrix a -> Matrix a
13multiply ((a,b),(c,d)) ((e,f),(g,h)) =
14 ((a*e+b*g, a*f+b*h), (c*e+d*g, c*f+d*h))
15
16identity :: Num a => Matrix a
17identity = ((1,0),(0,1))
18
19power :: Num a => Matrix a -> a -> Matrix a
20power _ 0 = identity
21power m n
22 | even n = power (multiply m m) (n `div` 2)
23 | otherwise = multiply m (power m (n - 1))
复杂度 O(log n),适合计算超大索引(如 fibonacci' 10^6)。矩阵乘法用元组表示,清晰但代码量大。
C++ 对照实现(Codeforces 风格)
1#include <bits/stdc++.h>
2using namespace std;
3using ll = long long;
4
5// 2x2 矩阵,a[i][j] 表示第 i 行第 j 列
6// 表示 | a[0][0] a[0][1] |
7// | a[1][0] a[1][1] |
8
9// 矩阵乘法:C = A * B
10// C[i][j] = sum_k A[i][k] * B[k][j]
11void mul(ll A[2][2], ll B[2][2], ll C[2][2]) {
12 ll tmp[2][2] = {{0, 0}, {0, 0}};
13 for (int i = 0; i < 2; i++) {
14 for (int j = 0; j < 2; j++) {
15 for (int k = 0; k < 2; k++) {
16 tmp[i][j] += A[i][k] * B[k][j];
17 }
18 }
19 }
20 // 结果写回 C
21 for (int i = 0; i < 2; i++) {
22 for (int j = 0; j < 2; j++) {
23 C[i][j] = tmp[i][j];
24 }
25 }
26}
27
28// 矩阵快速幂:把 A 变成 A^n,结果写在 res 里
29// 原理:二进制拆分 n
30// 若 n 的某一位是 1,就把当前的 A^{2^i} 乘进答案
31// 然后 A 自身平方:A -> A^2 -> A^4 -> ...
32void qpow(ll A[2][2], ll n, ll res[2][2]) {
33 // 单位矩阵 I:乘任何矩阵都不变
34 res[0][0] = 1; res[0][1] = 0;
35 res[1][0] = 0; res[1][1] = 1;
36
37 ll base[2][2];
38 for (int i = 0; i < 2; i++) {
39 for (int j = 0; j < 2; j++) {
40 base[i][j] = A[i][j];
41 }
42 }
43
44 while (n > 0) {
45 if (n & 1) {
46 // n 当前最低位是 1:res = res * base
47 ll tmp[2][2];
48 mul(res, base, tmp);
49 for (int i = 0; i < 2; i++) {
50 for (int j = 0; j < 2; j++) {
51 res[i][j] = tmp[i][j];
52 }
53 }
54 }
55 // base = base * base
56 ll tmp[2][2];
57 mul(base, base, tmp);
58 for (int i = 0; i < 2; i++) {
59 for (int j = 0; j < 2; j++) {
60 base[i][j] = tmp[i][j];
61 }
62 }
63 n >>= 1; // 处理下一位
64 }
65}
66
67// 用矩阵快速幂求第 n 个斐波那契数
68// 转移:
69// | F(n+1) | | 1 1 |^{n-1} | F(1) | | F(1)=1 |
70// | F(n) | = | 1 0 | * | F(0) | , | F(0)=0 |
71//
72// 所以 F(n) = (转移矩阵)^{n-1} 的左上角 a[0][0]
73ll fib(ll n) {
74 if (n < 1) return 0; // 题目要求正索引;按需改
75
76 // 转移矩阵 Q
77 ll Q[2][2] = {{1, 1}, {1, 0}};
78 ll P[2][2];
79 qpow(Q, n - 1, P);
80
81 // P * (1, 0)^T = (P[0][0], P[1][0])
82 // 其中 P[0][0] = F(n), P[1][0] = F(n-1)
83 return P[0][0];
84}
85
86int main() {
87 ios::sync_with_stdio(false);
88 cin.tie(nullptr);
89
90 ll n;
91 cin >> n;
92 cout << fib(n) << '\n';
93 return 0;
94}
| Haskell | C++ |
|---|---|
Matrix a = ((a,a),(a,a)) | ll A[2][2] 普通二维数组 |
multiply | mul |
power | qpow(二进制快速幂) |
apply 后取 fst | P[0][0] 即 F(n) |
说明:上面没取模。CF 上通常要 mod 的话,在 mul 里每次加法后 % MOD 即可。
方法二:使用线性递推的 fast doubling(推荐)
1fibonacci' :: Integral a => a -> a
2fibonacci' n
3 | n < 1 = error "fibonacci': index must be positive"
4 | otherwise = fst (fibPair n)
5 where
6 fibPair 0 = (0, 1)
7 fibPair k =
8 let (a, b) = fibPair (k `div` 2)
9 c = a * (b * 2 - a)
10 d = a * a + b * b
11 in if even k then (c, d) else (d, c + d)
Fast doubling 公式:
同样 O(log n),但无需矩阵,代码更紧凑,实际常数也更小。
方法三:调用 P29 的尾递归(朴素 O(n))
1fibonacci' :: Integral a => a -> a
2fibonacci' n = fibonacci n -- 复用 P29 的尾递归实现
对于 n 不大(< 10^6)的场景,O(n) 尾递归已经足够,代码最简单。
方法对比
| 方法 | 复杂度 | 代码量 | 适用场景 |
|---|---|---|---|
| 矩阵快速幂 | O(log n) | ~20 行 | 教学展示矩阵运算 |
| Fast doubling | O(log n) | ~8 行 | 生产推荐,简洁高效 |
| P29 尾递归 | O(n) | 1 行 | n 较小时最直接 |
测试
1>>> fibonacci' 100
2354224848179261915075
3>>> fibonacci' 20 == fibonacci 20
4True
5>>> fibonacci' 1
61