运用快速幂算法求解递推式的第N项

将递推式转化为矩阵式,并使用矩阵快速幂算法求解矩阵幂运算,从而解得原递推式的第N项值。

算法题目

这是一道来自蓝桥杯的题目。
问题描述
已知递推式:

初始值为:\( F(1,1) = 2 \) , \( F(1,2) = 3 \) , \( F(2,1) = 1 \) , \( F(2,2) = 4 \) , \( F(3,1) = 6 \) , \( F(3,2) = 5 \)。
输入 \( N \),输出 \( F(n,1) \) 和 \( F(n,2) \) ,由于答案可能很大,你只需要输出答案除以 99999999 的余数。
输入格式
输入第一行包括一个整数 \( N \)。
输出格式
输出两行,第一行为 \( F(n,1) \) 除以 99999999 的余数,第二行为 \( F(n,2) \) 除以 99999999 的余数。
样例输入
4
样例输出
14
21
数据规模和约定
\( 1 \leqslant N \leqslant 10^{18} \)

问题分析

  这个题目数据量很大,如果是简单地递归,肯定会栈溢出或者超时。即使把递归改成循环,也只能避免栈溢出,难逃超时的命运。
  递推式主要涉及加法运算,必须从头到尾把每一项都计算出来才可以得到答案,很难再做出优化。将递推式改为矩阵式,则涉及到矩阵的幂运算,而幂运算可以用 快速幂算法 优化,只需做 log(N) 次乘法运算就可以得到答案。

斐波那契数列的矩阵式

  我没有找到这个题目矩阵式推导的详细过程,所以只能自己琢磨。我知道一个很经典的递推式子,即斐波那契数列的递推式。研究怎么推导斐波那契数列的矩阵式,也许就能知道怎么推导这个蓝桥杯题目中的递推式的矩阵式。
  以下是我整理的斐波那契数列矩阵式的推导过程。

已知斐波那契数列的递推式

初始值为

设有矩阵 \( M \),使得下式成立

不妨设矩阵 \( M \)为

则由公式 \( \eqref{2} \) 可得

显然,\( C=1,D=0 \)。由公式 \( \eqref{1} \) 可知 \( A = B = 1 \)。
所以矩阵 \( M\) 为

于是公式 \( \eqref{2} \) 可以写为

将 \( F{n-1},F{n-2} \) 部分也替换为矩阵式,则有

进一步可以得到斐波那契数列的矩阵式

快速幂算法

  使用矩阵式来计算斐波那契数列的第 \( n \) 项,关键在于计算矩阵 \( M \) 的 \( n \) 次方。
  最简单地方法是连续做 \( n-1 \) 次乘法。考虑到这里其实是对矩阵进行幂运算,所以可以使用一种叫 快速幂 的算法来加速运算,只需要做 log(N) 次乘法运算就能得到结果。

  快速幂算法的思想是,把指数部分用 二进制 表示,分割出更小的任务。
  举个例子(这里直接使用了 OJ-Wiki.org 中的例子)

  因为 \( n \) 有 \( \lfloor \log_2 n \rfloor + 1 \) 个二进制位,因此当知道了 \( a^1,a^2,a^4,a^8,\cdots,a^{2 \lfloor \log_2 n \rfloor} \)后,就可以只计算 \( \Theta(\log n) \) 次乘法得到 \( a^n \)。

  由等式 \( \eqref{11} \) 可知,要计算 \( 3^{13} \),只需要将指数对应二进制位为 \( 1 \) 的整系数幂乘起来就可以了。

  循环求整数快速幂的代码实现

1
2
3
4
5
6
7
8
9
10
11
12
long long binpow(long long a, long long b) {
long long res = 1;
while (b > 0) {
// 如果当前二进制位为1,把系数幂乘进答案
if (b & 1) res = res * a;
// 求a^n
a = a * a;
// 整体右移,则最低位为下一个二进制位
b >>= 1;
}
return res;
}

  对于矩阵的快速幂算法,只需要把 long long 替换为自己的矩阵类,把整数乘法替换为矩阵乘法,把 res 的初始值替换为单位矩阵就可以了。
  未完全实现的矩阵快速幂:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// 矩阵类为简单实现,都为 maxn X maxn 的方阵
const int maxn = 8;
// 矩阵快速幂
Matrix quickPow(Matrix M, int e)
{
// 初始化ans为单位矩阵,eye()返回一个单位矩阵
Matrix ans = eye(),base = M;
// 快速幂运算
while (e > 0)
{
// MatrixMul 是对两个矩阵做乘法运算的函数
if (e & 1) ans = MatrixMul(ans, base);
base = MatrixMul(base, base);
e >>= 1;
}
return ans;
}

  有了矩阵快速幂算法,那么使用矩阵式来计算斐波那契数列的第 N 项已经不是问题了。

题目的矩阵式推导

  回到最开始的题目。题目中的递推式比斐波那契数列的递推式复杂的多,但是推导它的矩阵式并没有难多少。
同样地,假设存在矩阵 \( M \),使得下式成立

(注意:因为原 递推式 的右边存在 \( F(n-1,1),F(n-1,2),F(n-3,1),F(n-3,2),5,3 \) ,所以矩阵式的右边也会出现它们,将矩阵式右边的下标统统加 1,就得到了矩阵式的左边。)

设 \( M \) 为

由等式 \( \eqref{21} \) 可得一系列关系式,再结合递推式,可得矩阵 \( M \)为

……
(懒得写推导过程了,和斐波那契数列真的几乎一样)

附上这一题的解题代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#include <bits/stdc++.h>
using namespace std;
#define int int64_t
const int mod = 99999999ll;
const int maxn = 8;

struct Matrix
{
int M[maxn][maxn] = { 0 };
};

// 普通矩阵乘法
Matrix MatrixMul(Matrix A, Matrix B)
{
Matrix tM;
for (int i = 0; i < maxn; i++)
for (int j = 0; j < maxn; j++)
{
for (int k = 0; k < maxn; k++)
{
tM.M[i][j] += ( (A.M[i][k])%mod * (B.M[k][j])%mod )%mod;

}
tM.M[i][j] %= mod;
}
return tM;
}
// 矩阵快速幂
Matrix quickPow(Matrix M, int e)
{
Matrix ans,base = M;
// 初始化ans为单位矩阵
for (int i = 0; i < maxn; i++)
{
for (int j = 0; j < maxn; j++)
{
if (i == j) ans.M[i][j] = 1;
}
}
while (e)
{
if (e & 1) ans = MatrixMul(ans, base);
base = MatrixMul(base, base);
e >>= 1;
}
return ans;
}

signed main()
{
int T[maxn][maxn] =
{
{0,1,0,0,2,0,1,0},
{1,0,0,0,3,2,0,1},
{1,0,0,0,0,0,0,0},
{0,1,0,0,0,0,0,0},
{0,0,1,0,0,0,0,0},
{0,0,0,1,0,0,0,0},
{0,0,0,0,0,0,1,0},
{0,0,0,0,0,0,0,1}
};
int P[maxn] = { 6,5,1,4,2,3,5,3 };

int S1 = 2, S2 = 1, S3 = 6; // F(x,1) 用
int V1 = 3, V2 = 4, V3 = 5; // F(x,2) 用

Matrix M, ans;
memcpy(&M.M, T, sizeof(M.M));

int N;
cin >> N;


if (N == 1)
{
cout << S1 << endl << V1 << endl;
return 0;
}

if (N == 2)
{
cout << S2 << endl << V2 << endl;
return 0;
}

if (N == 3)
{
cout << S3 << endl << V3 << endl;
return 0;
}

ans = quickPow(M, N - 3); // 矩阵式从n>=4开始

int Fn1 = 0, Fn2 = 0;
for (int i = 0; i < maxn; i++)
{
Fn1 += (ans.M[0][i]%mod * P[i]%mod)%mod;
Fn2 += (ans.M[1][i]%mod * P[i]%mod)%mod;
}

cout << Fn1%mod << endl;
cout << Fn2%mod << endl;

return 0;
}

参考文献

[1] 快速幂 - OI Wiki. https://oi-wiki.org/math/quick-pow/

作者

uint128.com

发布于

2020-01-30

更新于

2022-08-22

许可协议

评论