JZOJ 6693. 【2020.06.05省选模拟】紫色彼岸樱推迟绽放 （自然数幂拆上升幂+组合意义矩阵乘法+NTT）
Description：
幽幽子饿了，妖梦需要给幽幽子准备食物。
有 T 天，每天幽幽子划分成了 k 个时段，妖梦需要安排每一天的日程。
第 i 天妖梦准备了 D+i-1 道菜，每道菜有无数个。第 1 个时段是早餐，幽幽子会选择 L 道不同的菜吃。
接下来 k-1 个时段，每个时段可以选择 D+i-1 道菜中的一道吃或者选择 A 个活动中的一个参加，但是出于健康考虑，幽幽子不能连续两个时段同时吃菜。
k 和 A 是幽幽子事先决定好的，她给出了 Q 个询问，每次询问给出 L,D,T，问这 T 天每一天的安排的方案数之和。
由于答案可能很大，输出答案对 998244353 取模之后的结果。
https://gmoj.net/senior/#main/show/6693
前置知识：
自然数幂的两种拆法：
1.\(n^m=\sum_{i=0}^m S(m,i)*(\binom{n}{i}*i!=n^{i↓})\)
这个是有组合意义的，\(n^m\)意义为\(m\)个求放\(n\)个盒子里，那么枚举放了至少一个球的盒子数\(i\)，分配即可。
2.\(n^m=\sum_{i=0}^m S(m,i)*n^{i↑}*(-1)^{m-i}\)
这个的意义不明，好像要用斯特林数的性质证明。
第二类斯特林数的容斥：
\(S(n,m)=\frac{1}{m!} \sum_{i=0}^m C(m,i)*(-1)^{m-i}*i^n\)
题解：
本题答案是（\(k--\)后）：
\(\sum_{n=0}^R \binom{n}{L} \sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i}*n^i\)
尝试把\(n^i\)拆成上升幂以搞掉前面的\(\binom{n}{L}\)。
\(n^i\)直接搞成上升幂\(n^{i↑}\)不太好弄，因为提不出\(n!\)，拆成\((n+1)^{i↑}\)比较好。
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{n=0}^R \binom{n}{L}*\sum_{j=0}^{k/2+1}S(i+1,j+1)*(-1)^{i-j}*(n+1)^{j↑}\)
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}*\sum_{n=0}^R \binom{n}{L}*(n+1)^{j↑}\)
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}*\sum_{n=0}^R \frac{n!*(n+j)!}{L!*(n-L)!*n!}\)
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}*\sum_{n=0}^R \frac{(n+j)!}{L!*(n-L)!}\)
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}*\sum_{n=0}^R \frac{(n+j)!*(L+j)!}{(L+i)!*(n-L)!*L!}\)
\(=\sum_{i=0}^{k/2} A^{k-i}*\binom{k-i}{i} \sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}* \binom{R+j+1}{L+j+1}\frac{(L+j)!}{L!}\)
然后我们终于把R搞掉了，为了接下的方便，换换元：
\(=\sum_{i=0}^{k/2} \binom{R+j+1}{L+j+1}\frac{(L+j)!}{L!} \sum_{m=0}^{k/2} S(m+1,i+1)*A^{k-m}*\binom{k-m}{m}\)
干掉斯特林数：
\(=\sum_{i=0}^{k/2} \binom{R+j+1}{L+j+1}\frac{(L+j)!}{L!} \sum_{m=0}^{k/2} A^{k-m}*\binom{k-m}{m} \frac{1}{(i+1)!}\sum_{j=0}^{k/2+1} \binom{i+1}{j}*(-1)^{i+1-j}*j^{m+1}\)
\(=\sum_{i=0}^{k/2} \binom{R+j+1}{L+j+1}\frac{(L+j)!}{L!} \frac{1}{(i+1)!}\sum_{j=0}^{k/2+1} \binom{i+1}{j}*(-1)^{i+1-j}*\sum_{m=0}^{k/2} A^{k-m}*\binom{k-m}{m}*j^{m+1}\)
对\(m\)那层循环，实际上可以对每个\(j\)快速算出，一种方法是多项式多点求值，另一种是发现它和一开始的式子长的很像，根据组合意义，可以写出一个矩阵乘法优化的dp。
然后发现再卷积一次就可以算出每个\(i\)后面的系数，那么每次枚举\(i\)即可\(O(k)\)计算一次询问。
Code
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i < _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("\n")
using namespace std;
const int mo = 998244353;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
namespace ntt {
const int nm = 1 << 18;
ll w[nm], a[nm], b[nm];
int r[nm];
void build() {
for(int i = 1; i < nm; i *= 2) {
w[i] = 1; ll v = ksm(3, (mo - 1) / 2 / i);
ff(j, 1, i) w[i + j] = w[i + j - 1] * v % mo;
}
}
void dft(ll *a, int n, int f) {
ff(i, 0, n) {
r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
if(i < r[i]) swap(a[i], a[r[i]]);
} ll v;
for(int i = 1; i < n; i *= 2) for(int j = 0; j < n; j += 2 * i) ff(k, 0, i) {
v = a[i + j + k] * w[i + k], a[i + j + k] = (a[j + k] - v) % mo, a[j + k] = (a[j + k] + v) % mo;
}
if(f == -1) {
reverse(a + 1, a + n);
v = ksm(n, mo - 2);
ff(i, 0, n) a[i] = (a[i] + mo) * v % mo;
}
}
}
using ntt :: dft;
const int M = 2e7 + 5;
ll fac[M], nf[M];
void build(int n) {
fac[0] = 1; fo(i, 1, n) fac[i] = fac[i - 1] * i % mo;
nf[n] = ksm(fac[n], mo - 2); fd(i, n, 1) nf[i - 1] = nf[i] * i % mo;
}
ll C(int n, int m) {
if(n < m) return 0;
return fac[n] * nf[m] % mo * nf[n - m] % mo;
}
int k, A, Q;
int L, D, T;
struct P {
ll a[2][2];
P() { memset(a, 0, sizeof a);}
};
P operator * (P a, P b) {
P c = P();
fo(k, 0, 1) fo(i, 0, 1) fo(j, 0, 1)
c.a[i][j] = (c.a[i][j] + a.a[i][k] * b.a[k][j]) % mo;
return c;
}
P ksm(P x, int y) {
P s = P(); s.a[0][0] = s.a[1][1] = 1;
for(; y; y /= 2, x = x * x)
if(y & 1) s = s * x;
return s;
}
const int N = 1 << 17;
int n;
ll f[N], g[N];
ll a[N], b[N];
void build() {
fo(j, 0, k / 2 + 1) {
P w = P();
w.a[0][0] = A; w.a[1][0] = A;
w.a[0][1] = -j;
w = ksm(w, k);
f[j] = (w.a[1][0] + w.a[1][1]) * j % mo;
}
int tp = 17;
fo(i, 0, k / 2 + 1) a[i] = nf[i] * (i % 2 ? -1 : 1);
fo(j, 0, k / 2 + 1) b[j] = f[j] * nf[j] % mo;
dft(a, 1 << tp, 1); dft(b, 1 << tp, 1);
ff(i, 0, 1 << tp) a[i] = a[i] * b[i] % mo;
dft(a, 1 << tp, -1);
fo(i, 0, k / 2) g[i] = a[i + 1];
}
ll solve(ll R, ll L) {
ll ans = 0;
fo(i, 0, k / 2)
ans = (ans + C(i + R + 1, i + L + 1) % mo * fac[L + i] % mo * (i % 2 ? -1 : 1) * g[i]) % mo;
ans = ans * nf[L] % mo;
return ans;
}
int main() {
ntt :: build();
freopen("bloom.in", "r", stdin);
freopen("bloom.out", "w", stdout);
build(2e7);
scanf("%d %d %d", &k, &A, &Q);
k --;
build();
fo(ii, 1, Q) {
scanf("%d %d %d", &L, &D, &T);
ll ans = solve(D + T - 1, L) - solve(D - 1, L);
ans = (ans % mo + mo) % mo;
pp("%lld\n", ans);
}
}
