JZOJ 6693. 2020.06.05省选模拟紫色彼岸樱推迟绽放 （自然数幂拆上升幂+组合意义矩阵乘法+NTT）

Posted 2020-12-07 coldchair

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)*(inom{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 inom{n}{L} sum_{i=0}^{k/2} A^{k-i}*inom{k-i}{i}*n^i)

尝试把(n^i)拆成上升幂以搞掉前面的(inom{n}{L})。

(n^i)直接搞成上升幂(n^{i↑})不太好弄，因为提不出(n!)，拆成((n+1)^{i↑})比较好。

(=sum_{i=0}^{k/2} A^{k-i}*inom{k-i}{i} sum_{n=0}^R inom{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}*inom{k-i}{i} sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}*sum_{n=0}^R inom{n}{L}*(n+1)^{j↑})

(=sum_{i=0}^{k/2} A^{k-i}*inom{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}*inom{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}*inom{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}*inom{k-i}{i} sum_{j=0}^{k/2}S(i+1,j+1)*(-1)^{i-j}* inom{R+j+1}{L+j+1}frac{(L+j)!}{L!})

然后我们终于把R搞掉了，为了接下的方便，换换元：

(=sum_{i=0}^{k/2} inom{R+j+1}{L+j+1}frac{(L+j)!}{L!} sum_{m=0}^{k/2} S(m+1,i+1)*A^{k-m}*inom{k-m}{m})

干掉斯特林数：

(=sum_{i=0}^{k/2} inom{R+j+1}{L+j+1}frac{(L+j)!}{L!} sum_{m=0}^{k/2} A^{k-m}*inom{k-m}{m} frac{1}{(i+1)!}sum_{j=0}^{k/2+1} inom{i+1}{j}*(-1)^{i+1-j}*j^{m+1})

(=sum_{i=0}^{k/2} inom{R+j+1}{L+j+1}frac{(L+j)!}{L!} frac{1}{(i+1)!}sum_{j=0}^{k/2+1} inom{i+1}{j}*(-1)^{i+1-j}*sum_{m=0}^{k/2} A^{k-m}*inom{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("
")
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
", ans);
	}
}