线性递推式

GitHub跳转原题关系图 Canvas 关系图返回列表

用生成函数求出连续前缀，再以特征多项式取模计算 x^l，最后卷积得到区间内的线性递推值。

OJ: shumeng

题目 ID: CSP201809E

难度:省选/NOI-

标签:线性递推生成函数多项式NTT

日期: 2026-07-31 16:21

形式化题目

给定 $m$ 阶线性递推

a_n=\sum_{i=1}^{\min(n,m)}k_i a_{n-i}\pmod {998244353},\qquad a_0=1,

以及下标区间 $[l,r]$ ，输出 $a_l, a_{l+1}, \dots, a_r$ 。其中 $m \le 10^5$ 、 $l \le 10^{12}$ ，区间长度 $r-l+1 \le 3\times 10^5$ 。

思路

朴素做法

先看一个直接按递推式逐项计算的朴素程序：

cpp

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/**
 * Author by Rainboy blog: https://rainboylv.com github: https://github.com/rainboylvx
 * rbook: -> https://rbook.roj.ac.cn  https://rbook2.roj.ac.cn
 * rainboy的学习导航网站: https://idx.roj.ac.cn
 * create_at: 2026-07-31 16:21
 * update_at: 2026-08-17 22:41
 */
// brute.cpp：直接按递推式逐项计算，只用于小下标对拍。
#include <bits/stdc++.h>
using namespace std;

const int MOD = 998244353;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int m;
    long long left, right;
    cin >> m >> left >> right;
    vector<int> coefficient(m + 1, 0); // k_1..k_m
    for (int i = 1; i <= m; i++) cin >> coefficient[i];

    // 按定义逐项计算 a[i] = sum k_j * a[i-j]，复杂度 O(mr)，只能跑小数据。
    vector<int> sequence(right + 1, 0);
    sequence[0] = 1;
    for (int i = 1; i <= right; i++) {
        long long value = 0;
        for (int j = 1; j <= m && j <= i; j++) {
            value += coefficient[j] * 1LL * sequence[i - j] % MOD;
            if (value >= MOD) value -= MOD;
        }
        sequence[i] = (int)value;
    }
    for (long long i = left; i <= right; i++) cout << sequence[i] << '\n';

    return 0;
}

它计算每个 $a_i$ 时枚举最近的 $m$ 项，复杂度 $O(mr)$ ，只能用于小下标和对拍。完整数据需要同时处理很大的 $l$ 与较长的输出区间。

生成函数得到前缀

令

C(x)=1-k_1x-k_2x^2-\cdots-k_mx^m.

递推式等价于生成函数

A(x)=\sum_{n\geq 0}a_nx^n=\frac{1}{C(x)}.

因此只要用多项式求逆求出 $A(x)$ 的前 $m+(r-l)-1$ 项，就能得到后面卷积所需的 $a_0$ 到 $a_{m+r-l-1}$ 。代码使用牛顿迭代

R_{new}=R(2-CR)\pmod {x^L}

每次把已知长度翻倍；乘法由 NTT 完成。实现中把一次迭代放在同一个频域里计算，减少一次正变换。

计算 $x^l$ 对特征多项式取模

定义递推的特征多项式

P(x)=x^m-k_1x^{m-1}-k_2x^{m-2}-\cdots-k_m。

在模 $P(x)$ 的意义下，所有多项式都可以压成次数小于 $m$ 的形式。二进制快速幂计算

R_l(x)=x^l\bmod P(x)=\sum_{j=0}^{m-1}c_jx^j。

每次乘法后需要做多项式取模：把 $V$ 和 $P$ 反转后，用预先求出的 $\operatorname{rev}(P)^{-1}$ 求商的反转形式，再乘回 $P$ 得到余式。 $P$ 与这个逆多项式在所有取模操作中固定不变，所以代码预先缓存它们的 NTT 结果。

一次卷积得到整个区间

线性递推的特征多项式关系说明，对任意 $t \ge 0$ 都有

a_{l+t}=\sum_{j=0}^{m-1}c_j a_{j+t}。

把 $c_j$ 逆序后与已求出的 $a_i$ 做一次卷积，卷积下标 $m-1+t$ 正好是上式的右侧。因此只需读取 $t=0,1,\ldots,r-l$ 的这些系数。

代码

cpp

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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239/**
 * Author by Rainboy blog: https://rainboylv.com github: https://github.com/rainboylvx
 * rbook: -> https://rbook.roj.ac.cn  https://rbook2.roj.ac.cn
 * rainboy的学习导航网站: https://idx.roj.ac.cn
 * create_at: 2026-07-31 16:21
 * update_at: 2026-08-17 22:41
 */
#include <bits/stdc++.h>
using namespace std;

const int MOD = 998244353;
const int ROOT = 3;

int ntt_limit;              // 一次 NTT 的长度
int ntt_log;                // ntt_limit 对应的指数（2^ntt_log = ntt_limit）
vector<int> ntt_root;       // 各次单位根旋转因子
vector<int> ntt_reverse;    // 蝶形变换所需的位逆序表
int reduction_ntt_size;     // 取模运算中固定乘法的 NTT 长度
vector<int> characteristic_transform;      // 特征多项式的 NTT 结果（预缓存）
vector<int> inverse_reversed_transform;    // rev(P)^{-1} 的 NTT 结果（预缓存）

int power_mod(long long base, int exponent) {
    long long result = 1;
    while (exponent > 0) {
        if (exponent & 1) result = result * base % MOD;
        base = base * base % MOD;
        exponent >>= 1;
    }
    return (int)result;
}

void prepare_ntt(int minimum_limit) {
    ntt_limit = 1;
    ntt_log = 0;
    while (ntt_limit < minimum_limit) {
        ntt_limit <<= 1;
        ntt_log++;
    }

    ntt_root.assign(ntt_limit, 1);
    int primitive = power_mod(ROOT, (MOD - 1) / ntt_limit);
    for (int i = 1; i < ntt_limit; i++) {
        ntt_root[i] = (int)(ntt_root[i - 1] * 1LL * primitive % MOD);
    }

    ntt_reverse.assign(ntt_limit, 0);
    for (int i = 1; i < ntt_limit; i++) {
        ntt_reverse[i] = (ntt_reverse[i >> 1] >> 1) | ((i & 1) << (ntt_log - 1));
    }
}

void ntt(vector<int> &value, bool inverse) {
    int n = (int)value.size();
    int current_log = __builtin_ctz((unsigned int)n);
    int reverse_shift = ntt_log - current_log;
    for (int i = 1; i < n; i++) {
        int j = ntt_reverse[i] >> reverse_shift;
        if (i < j) swap(value[i], value[j]);
    }
    for (int half = 1; half < n; half <<= 1) {
        int step = ntt_limit / (half << 1);
        for (int start = 0; start < n; start += half << 1) {
            for (int i = 0; i < half; i++) {
                int left = value[start + i];
                int root_index = i * step;
                if (inverse && root_index != 0) root_index = ntt_limit - root_index;
                int right = (int)(value[start + i + half] * 1LL * ntt_root[root_index] % MOD);
                value[start + i] = left + right < MOD ? left + right : left + right - MOD;
                value[start + i + half] = left - right < 0 ? left - right + MOD : left - right;
            }
        }
    }
    if (inverse) {
        int inv_n = power_mod(n, MOD - 2);
        for (int i = 0; i < n; i++) value[i] = (int)(value[i] * 1LL * inv_n % MOD);
    }
}

vector<int> multiply(const vector<int> &left, const vector<int> &right) {
    if (left.empty() || right.empty()) return vector<int>();
    if ((int)left.size() * (int)right.size() <= 4096) {
        vector<int> result(left.size() + right.size() - 1, 0);
        for (int i = 0; i < (int)left.size(); i++) {
            for (int j = 0; j < (int)right.size(); j++) {
                result[i + j] = (result[i + j] + left[i] * 1LL * right[j]) % MOD;
            }
        }
        return result;
    }
    int size = 1;
    int need = (int)left.size() + (int)right.size() - 1;
    while (size < need) size <<= 1;
    vector<int> first(left), second(right);
    first.resize(size);
    second.resize(size);
    ntt(first, false);
    ntt(second, false);
    for (int i = 0; i < size; i++) first[i] = (int)(first[i] * 1LL * second[i] % MOD);
    ntt(first, true);
    first.resize(need);
    return first;
}

vector<int> multiply_with_fixed_right(const vector<int> &left, int right_size,
                                      const vector<int> &right_transform) {
    int need = (int)left.size() + right_size - 1;
    // 右侧多项式已经完成 NTT，只需变换左侧并做一次逆变换。
    vector<int> value(left);
    value.resize(reduction_ntt_size);
    ntt(value, false);
    for (int i = 0; i < reduction_ntt_size; i++) {
        value[i] = (int)(value[i] * 1LL * right_transform[i] % MOD);
    }
    ntt(value, true);
    value.resize(need);
    return value;
}

vector<int> polynomial_inverse(const vector<int> &polynomial, int need) {
    vector<int> result(1, power_mod(polynomial[0], MOD - 2));
    while ((int)result.size() < need) {
        int length = min(need, (int)result.size() << 1);
        vector<int> prefix(polynomial.begin(), polynomial.begin() + min((int)polynomial.size(), length));
        int transform_size = 1;
        // 同一频域中还要完成 R * (2 - A * R)，长度需覆盖完整的中间乘积，避免循环卷积串入低位。
        while (transform_size < (int)prefix.size() + 2 * (int)result.size() - 1) transform_size <<= 1;
        vector<int> prefix_transform(prefix);
        vector<int> result_transform(result);
        prefix_transform.resize(transform_size);
        result_transform.resize(transform_size);
        ntt(prefix_transform, false);
        ntt(result_transform, false);
        for (int i = 0; i < transform_size; i++) {
            int product = (int)(prefix_transform[i] * 1LL * result_transform[i] % MOD);
            int correction = 2 - product;
            if (correction < 0) correction += MOD;
            result_transform[i] = (int)(result_transform[i] * 1LL * correction % MOD);
        }
        ntt(result_transform, true);
        result_transform.resize(length);
        result = result_transform;
    }
    return result;
}

vector<int> polynomial_remainder(const vector<int> &value, const vector<int> &characteristic,
                                 const vector<int> &inverse_reversed) {
    int degree = (int)characteristic.size() - 1;
    if ((int)value.size() <= degree) return value;
    int quotient_size = (int)value.size() - degree;
    vector<int> reversed_value(value.rbegin(), value.rend());
    reversed_value.resize(quotient_size);
    vector<int> reversed_quotient;
    if ((int)reversed_value.size() * (int)inverse_reversed.size() <= 4096) {
        reversed_quotient = multiply(reversed_value, inverse_reversed);
    } else {
        reversed_quotient = multiply_with_fixed_right(reversed_value, inverse_reversed.size(),
                                                      inverse_reversed_transform);
    }
    reversed_quotient.resize(quotient_size);
    vector<int> quotient(reversed_quotient.rbegin(), reversed_quotient.rend());
    vector<int> product;
    if ((int)quotient.size() * (int)characteristic.size() <= 4096) {
        product = multiply(quotient, characteristic);
    } else {
        product = multiply_with_fixed_right(quotient, characteristic.size(), characteristic_transform);
    }
    vector<int> result(degree, 0);
    for (int i = 0; i < degree; i++) {
        int subtract = i < (int)product.size() ? product[i] : 0;
        result[i] = value[i] - subtract < 0 ? value[i] - subtract + MOD : value[i] - subtract;
    }
    return result;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int m;
    long long left, right;
    cin >> m >> left >> right;
    vector<int> coefficient(m + 1, 0); // k_1..k_m
    for (int i = 1; i <= m; i++) cin >> coefficient[i];
    int count = (int)(right - left + 1);
    prepare_ntt(2 * (m + count + 5));

    // 生成函数分母 C(x)=1-k1*x-...-km*x^m，序列前缀由 1/C(x) 的前若干项给出。
    vector<int> denominator(m + 1, 0);
    denominator[0] = 1;
    for (int i = 1; i <= m; i++) denominator[i] = coefficient[i] == 0 ? 0 : MOD - coefficient[i];

    // 递推的特征多项式 P(x)=x^m-k1*x^(m-1)-...-km，以及 rev(P) 的逆。
    vector<int> characteristic(m + 1, 0);
    characteristic[m] = 1;
    for (int i = 0; i < m; i++) {
        characteristic[i] = coefficient[m - i] == 0 ? 0 : MOD - coefficient[m - i];
    }
    vector<int> reversed_characteristic(characteristic.rbegin(), characteristic.rend());
    vector<int> inverse_reversed = polynomial_inverse(reversed_characteristic, m);

    // 取模操作反复用到 P 与 rev(P)^{-1}，它们各自只做一次 NTT。
    reduction_ntt_size = 1;
    while (reduction_ntt_size < 2 * m - 1) reduction_ntt_size <<= 1;
    characteristic_transform = characteristic;
    characteristic_transform.resize(reduction_ntt_size);
    ntt(characteristic_transform, false);
    inverse_reversed_transform = inverse_reversed;
    inverse_reversed_transform.resize(reduction_ntt_size);
    ntt(inverse_reversed_transform, false);

    // 二进制快速幂求 x^left mod P(x)，乘法后都要取模压回次数小于 m。
    vector<int> result(m, 0), base(m, 0);
    result[0] = 1;
    if (m == 1) base[0] = coefficient[1];
    else base[1] = 1;
    long long exponent = left;
    while (exponent > 0) {
        if (exponent & 1) {
            result = polynomial_remainder(multiply(result, base), characteristic, inverse_reversed);
        }
        exponent >>= 1;
        if (exponent > 0) {
            base = polynomial_remainder(multiply(base, base), characteristic, inverse_reversed);
        }
    }

    // 求出序列前缀 a[0..m+count-2]。
    int prefix_length = m + count - 1;
    vector<int> sequence = polynomial_inverse(denominator, prefix_length);
    // 特征多项式关系 a[l+t] = sum c_j * a[j+t]，把 c 逆序后与序列卷积即可整段读出。
    vector<int> reversed_result(result.rbegin(), result.rend());
    vector<int> convolution = multiply(reversed_result, sequence);
    for (int i = 0; i < count; i++) {
        cout << convolution[m - 1 + i] << '\n';
    }

    return 0;
}

复杂度

设 $T=r-l+1$ ，记长度为 $N$ 的 NTT 乘法复杂度为 $M(N)=O(N\log N)$ 。生成函数求逆和最后的卷积共需 $O(M(m+T))$ ，二进制计算 $x^l\bmod P(x)$ 需 $O(M(m)\log l)$ 。总时间复杂度为

O\bigl(M(m+T)+M(m)\log l\bigr)，

空间复杂度为 $O(m+T)$ 。

总结

面对“递推阶数很高、下标很大、还要输出一段连续值”的问题，可以先把递推写成生成函数，再把大下标转化为特征多项式上的幂。前缀由多项式求逆得到，大下标由 $x^l$ 取模处理，最后用一次卷积同时恢复整段答案。

由于原题 1 秒的时限相对于本题复杂度上限极其不合理，本评测时限改为 5 秒，同时数据范围进行了一些微调。但是我们的 std 之一可以保证在 1 秒内同时通过官方数据和民间数据。

时间限制： ~~1.0 秒~~ 5.0 秒

空间限制： 256 MB

题目描述

给出取模意义下的线性递推式：

a_n \equiv \sum_{i=1}^{\min(n,m)}{k_ia_{n-i} } \pmod Q

上式对 $n \in \mathbb{N}^+$ 成立，有 $a_0=1$ ，并且满足 $0\le a_n < Q$ 。其中 $m,k_1,k_2,\dots,k_m,Q$ 为给定的非负整数， $m,Q > 0$ ， $a \equiv b \pmod Q$ 表示 $a$ 和 $b$ 除以 $Q$ 的余数相等。

已知 $Q=998244353$ 。给出非负整数 $l,r\ (l\le r)$ ，求 $a_l, a_{l+1}, \dots, a_r$ 对 $Q$ 取模的结果。

输入格式

从标准输入读入数据。

输入的第一行包含 $3$ 个非负整数 $m,l,r$ 。其中 $l \le r$ 。

第二行包含 $m$ 个非负整数 $k_1,k_2,\dots,k_m$ ，保证 $0 \le k_1,k_2,\dots,k_m < Q$ 。

输出格式

输出到标准输出。

输出 $r-l+1$ 行，每行一个正整数，分别表示 $a_l, a_{l+1},\dots,a_r$ 对 $Q$ 取模的结果。

input1

1
23 3 6
2 0 4

output1

样例 1 解释

$k_1 = 2,k_2=0,k_3=4$ 。需要求出 $a_3,a_4,a_5,a_6$ 。

$a_1 = k_1 a_0 = 2 \times 1 = 2$ ；

$a_2 = k_1 a_1 + k_2 a_0 = 2 \times 2 + 0 \times 1 = 4$ ；

$a_3 = k_1 a_2 + k_2 a_1 + k_3 a_0 = 2 \times 4 + 0 \times 2 + 4 \times 1 = 12$ ；

$a_4 = k_1 a_3 + k_2 a_2 + k_3 a_1 = 2 \times 12 + 0 \times 4 + 4 \times 2 = 32$ ；

$a_5 = k_1 a_4 + k_2 a_3 + k_3 a_2 = 2 \times 32 + 0 \times 12 + 4 \times 4 = 80$ ；

$a_6 = k_1 a_5 + k_2 a_4 + k_3 a_3 = 2 \times 80 + 0 \times 32 + 4 \times 12 = 208$ ；

input2

1
22 1 11
1 1

output2

样例 2 解释

因为 $k_1=k_2=1$ ，因此这组样例就是斐波那契数列 $a_n = a_{n-1} + a_{n-2}$ 。

input3

1
210 10 20
532737790 634932889 335818534 101179174 977780682 695192541 779962395 295668292 157661238 325351676

output3

1
2
3
4
5
6
7
8
9
10
11119744921
651421717
601080475
163399777
291546699
108479226
406175654
344671679
459752012
489415425
349454810

子任务

对于所有数据，保证 $1\le m\le 10^5,1\le l\le 10^{12},0\le r-l\le 3\times 10^5$ 。

子任务 1（97 分）：传统计分，每组测试点满足以下限制：

测试点编号	$m=r-l=$	$l=$
$1$	$10$	$10$
$2$	$10^3$	$10^3$
$3\sim 4$	$10^5$	$10^5$
$5\sim 6$	$10^2$	$10^{12}$
$7\sim 8$	$10^3$	$10^{12}$
$9\sim 10$	$10^5$	$10^{12}$

子任务 2（3 分）：捆绑测试，无特殊限制。