理论

题目

K 阶斐波那契数列,f(n) = f(n-1) + f(n-2) + … + f(n-k),f(0) = f(1) = … = f(k-1) = 1, q 个询问,每次询问一个数 x,求 f(x) 的值,结果对 1000000007 取模。

针对这道题,通常有以下两个方向的解法,取决于具体的输入规模。

方案 1:数学推导优化 DP (滑动窗口法)

适用场景:$x$ 的最大值在内存和时间允许的范围内(例如 $\max(x) \le 10^7$)。

朴素 DP 分析
如果直接按公式 $f(n) = \sum_{i=1}^{k} f(n-i)$ 计算,计算一个状态需要 $O(k)$ 的时间,总时间复杂度是 $O(x \cdot k)$,这在 $x$ 和 $k$ 都比较大时必然会超时。

优化推导($O(1)$ 状态转移)
我们可以写出相邻两项的展开式:

将两式相减:

移项后得到递推式:

时间复杂度

  • 预处理阶段:我们可以用 $O(\max(x))$ 的时间,利用上述 $O(1)$ 的转移方程,把截止到最大 $x$ 的所有结果计算出来并存在数组中。
  • 查询阶段:对于 $q$ 个询问,每次只需 $O(1)$ 查表即可。
  • 总时间复杂度:$O(\max(x) + q)$。

注意:在代码实现时,减法取模需要写成 (2 * f[n-1] % MOD - f[n-k-1] % MOD + MOD) % MOD 以防止出现负数。

方案 2:矩阵快速幂

适用场景:$x$ 非常巨大(例如 $x \le 10^{18}$),无法开出那么大的数组,但 $k$ 相对较小(例如 $k \le 100$)。

理论推导
所有的常系数线性递推都可以转化为矩阵乘法。我们需要构造一个 $k \times k$ 的状态转移矩阵 $M$。

令状态列向量 $V_n = [f(n), f(n-1), \dots, f(n-k+1)]^T$。我们可以构造如下递推关系:

或者,我们也可以利用方案 1 推导出的 $f(n) = 2f(n-1) - f(n-k-1)$,构造一个 $(k+1) \times (k+1)$ 的更稀疏的矩阵,两者皆可。

根据矩阵乘法的结合律,$V_n = M^{n-(k-1)} \times V_{k-1}$。我们只需要利用快速幂算法在 $O(\log x)$ 的时间内求出 $M^{n-(k-1)}$ 即可。

矩阵快速幂,ref: https://zhuanlan.zhihu.com/p/137677246

时间复杂度

  • 单次询问:矩阵乘法耗时 $O(k^3)$,快速幂耗时 $O(\log x)$,单次为 $O(k^3 \log x)$。
  • 总时间复杂度:$O(q \cdot k^3 \log x)$。

Trade-off 分析

  1. 如果 $x \le 10^7$ 且 $q = 10^5$:选择方案 1(预处理查表)。因为方案 2 的 $O(q \cdot k^3 \log x)$ 会因为 $q$ 太大而超时。
  2. 如果 $x \le 10^{18}$ 且 $q = 100$,$k \le 50$:选择方案 2(矩阵快速幂)。因为 $x$ 太大,数组存不下,$O(N)$ 遍历也会超时。
  3. 如果 $x$ 很大,$k$ 也很大(例如 $x \le 10^{18}, k \le 10^5$):使用特征多项式取模 / Berlekamp-Massey 算法

练习

AT_abc401_c K-bonacci

https://atcoder.jp/contests/abc401/tasks/abc401_c

给定正整数 $N$ 和 $K$。我们按照以下方式定义长度为 $N+1$ 的数列 $A=(A_0,A_1,\ldots,A_N)$ 的每个元素值:

  • 当 $0 \leq i < K$ 时,$A_i = 1$
  • 当 $K \leq i$ 时,$A_i = A_{i-K} + A_{i-K+1} + \ldots + A_{i-1}$

请计算 $A_N$ 对 $10^9$ 取模后的结果。( 1 ≤ N,K ≤ 10^6 )

tle 代码:

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
#include <iostream>
using std::cin;
using std::cout;

int a[1000005];
int n, k;
const int MOD = 1e9;

int work(int x) {
if(a[x]) return a[x];
if(x < k) return 1;
for(int i = 1; i <= k ;++i) {
a[x] = (a[x] + work(x-i)) % MOD;
}
return a[x];
}

int main()
{
cin >> n >> k;
for(int i = 0; i < k; ++i) {
a[i] = 1;
}
cout << work(n);
return 0;
}

ac 代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include <iostream>
using std::cin;
using std::cout;

int n, k, a[1000005];
const int MOD = 1e9;

int main()
{
cin >> n >> k;
for(int i = 0; i < k; ++i) {
a[i] = 1;
}
a[k] = k;
for(int i = k+1; i <= n; ++i) {
a[i] = (2*a[i-1]%MOD - a[i-k-1]%MOD + MOD) % MOD;
}
cout << a[n];
return 0;
}

不过本题最自然的想法是用前缀和维护。

P1962 斐波那契数列

https://www.luogu.com.cn/problem/P1962

请你求出 $F_n \bmod 10^9 + 7$ 的值。$1\le n < 2^{63}$ 。

对于本题,有:

ac 代码:

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
#include <cstdint>
#include <iostream>
#include <cstring>

constexpr int64_t MOD = 1000000007;

class Matrix {
public:
Matrix() { memset(a, 0, sizeof(a)); }
Matrix operator * (const Matrix& rhs) const {
Matrix res;
for(int i = 1; i < 3; ++i)
for(int j = 1; j < 3; ++j) {
int64_t sum = 0;
for(int k = 1; k < 3; ++k)
// 两个 int32_t 相乘可能会溢出,所以先转成 int64_t
sum = (sum + (int64_t)a[i][k]*rhs.a[k][j]) % MOD;
res.a[i][j] = sum;
}
return res;
}

int32_t a[3][3];
};

Matrix matrix_pow(Matrix& lhs, int64_t n) {
Matrix res;
for(int i = 1; i < 3; ++i) res.a[i][i] = 1; // 单位矩阵
while(n) {
if(n & 1) res = res * lhs;
lhs = lhs * lhs;
n >>= 1;
}
return res;
}

int main()
{
int64_t n; std::cin >> n;
if(n <= 2) { std::cout << "1"; return 0;}

Matrix m;
m.a[1][1] = m.a[1][2] = m.a[2][1] = 1;
Matrix tmp = matrix_pow(m, n-2);
int64_t ans = ((tmp.a[1][1]%MOD) + (tmp.a[1][2]%MOD)) % MOD;
std::cout << ans;
return 0;
}

/*
test case:
input: 5, output: 5
input: 10, output: 55
input: 471857495, output: 495435152
*/