P30 - 矩阵快速幂计算斐波那契数

2026-09-07 00:00    #Haskell   #99题   #数论  

P30 - 矩阵快速幂计算斐波那契数

Fibonacci numbers with matrix exponentiation

官方模块:Problems.P30 核心函数:fibonacci'


← P29 斐波那契数 | P31 判断素数 →


核心公式

利用转移矩阵

(F(n+1)F(n))=(1110)n1(10) \begin{pmatrix} F(n+1) \\ F(n) \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^{n-1} \begin{pmatrix} 1 \\ 0 \end{pmatrix}

问题变成整数矩阵的快速幂。

函数签名

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}
HaskellC++
Matrix a = ((a,a),(a,a))ll A[2][2] 普通二维数组
multiplymul
powerqpow(二进制快速幂)
apply 后取 fstP[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 doublingO(log n)~8 行生产推荐,简洁高效
P29 尾递归O(n)1 行n 较小时最直接

测试

1>>> fibonacci' 100
2354224848179261915075
3>>> fibonacci' 20 == fibonacci 20
4True
5>>> fibonacci' 1
61

参考