Statement
给出 \(k,p,L\)，数序列 \(a\)，满足如下条件：
- \(1\le a_i\le k\)
- \(\sum_i a_i=L\)
- \(\nexists i,a_i\ge p\land a_{i+1}\ge p\)
答案对 \(20201114\) 取模，\(p\le k\le 1000, L\le 10^{18}\).
Solution
30 pts
注意到可以 dp，记 \(f(i,0/1)\) 为凑出 \(i\) 的方案数，其中最后一个数是否 \(<p\).
\[ \begin{aligned} f(i,1)&=\sum_{1\le j<p} f(i-j,1)+f(i-j,0)\\ f(i,0)&=\sum_{p\le j\le k} f(i-j,1)\\ f(0,1)&=1 \end{aligned} \]答案为 \(f(L,0)+f(L,1)\)，时间 \(O(kl)\)，预计得分 30 pts.
我做了个没什么用的前缀和优化：记 \(g(i,0/1)\) 为 \(f(i,0/1)\) 的前缀和，推式子可得
答案为 \(g(L,0)-g(L-1,0)+g(L,1)-g(L-1,1)\)，时间 \(O(L)\)，预计得分 30 pts.
60 pts
用矩阵快速幂优化上述 dp，\(f\) 或 \(g\) 都可以，时间 \(O(k^3\log L)\)，带有 8 倍常数，预计得分 60 pts.
100 pts
考虑把 \(f(i,0)\) 拆成若干 \(f(i,1)\)，这样就可以线性递推了！
不用写任意模数乘法，直接暴力乘就行，预处理系数也可以直接暴力，这是 \(O(k^2\log L)\) 的，期望 100 pts.
也可以用 \(g\) 算，把 \(g(i,0)\) 拆成若干 \(g(i,1)\) 就行了，预处理系数时会更优一点，需要稍微推一下式子。
你当然可以再写一个任意模数乘法做到 \(O(k\log k\log L)\)。
Code
#include <bits/stdc++.h>
using namespace std;
#define rep(i, j, k) for (int i = j; i <= k; ++i)
#define reo(i, j, k) for (int i = j; i >= k; --i)
typedef long long ll;
const int N = 2e5 + 10, P = 20201114;
int k, p, n, xishu[N], f[N][2];
ll L;
ll Mod(ll x) {
return (x % P + P) % P;
}
struct Poly {
vector<ll> f;
ll & operator[] (int x) {
return f[x];
}
Poly() {
f.clear();
}
Poly operator* (Poly rhs) {
Poly Res;
int Slhs = f.size(), Srhs = rhs.f.size();
Res.f.assign(Slhs + Srhs - 1, 0ll);
rep(i, 0, Slhs - 1)
rep(j, 0, Srhs - 1)
(Res[i + j] += f[i] * rhs[j] % P) %= P;
return Res;
}
void Modulo() {
int sz = f.size();
if (sz < p + k) {
return;
}
reo(i, sz - 1, p + k - 1) {
ll x = f[i];
if (x) {
rep(j, 1, p + k - 1) {
(f[i - j] += x * xishu[j] % P) %= P;
}
f[i] = 0;
}
}
f.resize(p + k);
}
void print() {
for (auto v : f) cout << v << ' ';
cout << '\n';
}
} base, Res;
int main() {
// freopen("jump.in", "r", stdin), freopen("jump.out", "w", stdout);
ios::sync_with_stdio(false), cin.tie(nullptr);
cin >> k >> p >> L;
f[0][1] = 1;
rep(i, 1, k + p) {
rep(j, 1, p - 1) {
if (i - j >= 0) (f[i][1] += f[i - j][1] + f[i - j][0]) %= P;
else break;
}
rep(j, p, k) {
if (i - j >= 0) (f[i][0] += f[i - j][1]) %= P;
else break;
}
}
rep(i, 1, p - 1) {
++xishu[i];
rep(j, p, k) {
++xishu[i + j];
}
}
if (L <= k) {
ll res = f[L][1];
rep(i, L - k, L - p) {
if (i >= 0) (res += f[i][1]) %= P;
else break;
}
cout << res << '\n';
return 0;
}
ll res = 0, Target = L - k - 1;
base.f.assign(2, 0ll), base[1] = 1;
Res.f.assign(1, 0ll), Res[0] = 1;
for (; Target; Target >>= 1, base = base * base, base.Modulo()) {
if (Target & 1) Res = Res * base, Res.Modulo();
}
rep(i, 0, k + p - 1) {
if (i >= (int)Res.f.size()) break;
(res += Res[i] * f[i + 1][1] % P) %= P;
}
rep(i, 1, k - p) {
int sz = Res.f.size() + 1;
Res.f.resize(sz);
reo(i, sz, 1) Res[i] = Res[i - 1];
Res[0] = 0;
Res.Modulo();
rep(i, 0, k + p - 1) {
if (i >= (int)Res.f.size()) break;
(res += Res[i] * f[i + 1][1] % P) %= P;
}
}
rep(i, 1, p) {
int sz = Res.f.size() + 1;
Res.f.resize(sz);
reo(i, sz, 1) Res[i] = Res[i - 1];
Res[0] = 0;
Res.Modulo();
}
rep(i, 0, k + p - 1) {
if (i >= (int)Res.f.size()) break;
(res += Res[i] * f[i + 1][1] % P) %= P;
}
cout << Mod(res) << '\n';
return 0;
}