【ÑÖÏ模拟赛】花萎（矩阵加速，循环卷积，高斯消元）

题面

4s ，512mb

题解

设 $dp[x][y]$ 表示到达 $x,y$ 的期望时间，$dp[0][0]=0$，容易得到转移式
$$dp[x][y]=p\ast dp[x-1\mathrm{m}\mathrm{o}\mathrm{d}n][y]+q\ast dp[x][y-1\mathrm{m}\mathrm{o}\mathrm{d}m]+1$$

但是由于 $x$ 过大，而 $y<400$ ，所以这不是我们想要的式子，我们需要主元少些方便高斯消元。我们得将 $dp[x][y-1]$ 展开，一直展开到无穷，没有 $x$ 项为止：
$$dp[x][y]=\frac{1}{1-q}+\sum _{i=0}^{m-1}dp[x-1][y-i\mathrm{m}\mathrm{o}\mathrm{d}m]\ast \frac{p{q}^i}{1-q^m}$$

这是个从 { { d p [ x − 1 ] [ y ] } , 1 } \{\{dp[x-1][y]\},1\} {{dp[x−1][y]},1} 到 { { d p [ x ] [ y ] } , 1 } \{\{dp[x][y]\},1\} {{dp[x][y]},1} 的线性变换，所以我们可以用矩阵乘法。算出转移矩阵的 $n-1$ 次方后，定 $dp[0][0\sim m-1]$ 为主元，我们就可以得到 $m$ 个方程的系数，然后高斯消元求出主元。

过程中预处理出转移矩阵的 $1,2,4,...2^{\log n}$ 次方，在询问的时候倍增求 $dp[X_i][Y_i]$，矩阵乘法可以从 $O(m^3)$ 变为 $O(m^2)$ ，老套路了。

这样一来，时间复杂度是 $O(m^3\log n+t{m}^2\log n)$ ，刚好过不了。

回过头来观察式子，我们发现 $\frac{p{q}^i}{1-q^m}$ 跟 $y$ 没有关系，也就是说，这就是个循环卷积。

令 $f_k={\sum }_{i=0}^{m-1}dp[k][i]\cdot x^i,g={\sum }_{i=0}^{m-1}\frac{p{q}^i}{1-q^m}\cdot x^i,h={\sum }_{i=0}^{m-1}\frac{1}{1-q}\cdot x^i$ ，那么
$$f_k=f_{k-1}\ast g+h$$

这是关于多项式的线性变换，我们打包 2×2 的多项式矩阵进行矩阵乘法，用 NTT 卷积模拟循环卷积，一次变换的复杂度就从 $O(m^3)$ 变成了 $O(m\log m)$ 。

最终复杂度为 $O(m^3+tm\log m\log n)$ 。

CODE

#include<map>
#include<set>
#include<cmath>
#include<ctime>
#include<queue>
#include<stack>
#include<random>
#include<bitset>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<unordered_map>
using namespace std;
#define MAXN 405
#define LL long long
#define ULL unsigned long long
#define ENDL putchar('\n')
#define DB double
#define lowbit(x) (-(x) & (x))
#define FI first
#define SE second
int xchar() {
	static const int maxn = 1000000;
	static char b[maxn];
	static int pos = 0,len = 0;
	if(pos == len) pos = 0,len = fread(b,1,maxn,stdin);
	if(pos == len) return -1;
	return b[pos ++];
}
// #define getchar() xchar()
LL read() {
	LL f = 1,x = 0;int s = getchar();
	while(s < '0' || s > '9') {if(s<0)return -1;if(s=='-')f=-f;s = getchar();}
	while(s >= '0' && s <= '9') {x = (x<<1) + (x<<3) + (s^48);s = getchar();}
	return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar((x%10)^48);}
void putnum(LL x) {
	if(!x) {putchar('0');return ;}
	if(x<0) putchar('-'),x = -x;
	return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}

const int MOD = 998244353;
int n,m,s,o,k;
inline int qkpow(int a,int b) {
	int res = 1;
	while(b > 0) {
		if(b & 1) res = res*1ll*a % MOD;
		a = a *1ll* a % MOD; b >>= 1;
	}return res;
}
int xm[MAXN<<2],om,rev[MAXN<<2];
void NTT(int *s,int n,int op) {
	for(int i = 1;i < n;i ++) {
		rev[i] = (rev[i>>1]>>1) | ((i&1) ? (n>>1):0);
		if(rev[i] < i) swap(s[rev[i]],s[i]);
	} om = qkpow(3,(MOD-1)/n); xm[0] = 1;
	if(op < 0) om = qkpow(om,MOD-2);
	for(int i = 1;i <= n;i ++) xm[i] = xm[i-1]*1ll*om % MOD;
	for(int k = 2,t = n>>1;k <= n;k <<= 1,t >>= 1) {
		for(int j = 0;j < n;j += k) {
			for(int i = j,l = 0;i < j+(k>>1);i ++,l += t) {
				int A = s[i],B = s[i+(k>>1)];
				s[i] = (A + xm[l]*1ll*B) % MOD;
				s[i+(k>>1)] = (A +MOD- xm[l]*1ll*B % MOD) % MOD;
			}
		}
	}
	if(op < 0) {
		LL v = qkpow(n,MOD-2);
		for(int i = 0;i < n;i ++) s[i] = s[i] * v % MOD;
	} return ;
}
int ta[MAXN<<2],tb[MAXN<<2];
void addmult(int *C,int *A,int *B,int m) {
	int le = 1; while(le < m+m) le <<= 1;
	for(int i = 0;i < le;i ++) ta[i] = tb[i] = 0;
	for(int i = 0;i < m;i ++) ta[i] = A[i],tb[i] = B[i];
	NTT(ta,le,1); NTT(tb,le,1);
	for(int i = 0;i < le;i ++) ta[i] = ta[i] *1ll* tb[i] % MOD;
	NTT(ta,le,-1);
	for(int i = 0,j = 0;i < le;i ++,j = (j+1 == m ? 0:j+1)) {
		C[j] += ta[i]; if(C[j] >= MOD) C[j] -= MOD;
	} return ;
}
int p,q;
struct mat{
	int s[2][2][401],n,m;
	mat(){memset(s,0,sizeof(s));n=m=0;}
}A,B,b[35];
mat operator * (mat a,mat b) {
	mat c; c.n = a.n;c.m = b.m;
	for(int i = 0;i < a.n;i ++) {
		for(int k = 0;k < a.m;k ++) {
			if(a.s[i][k])
			for(int j = 0;j < b.m;j ++) {
				addmult(c.s[i][j],a.s[i][k],b.s[k][j],m);
			}
		}
	} return c;
}
int a[MAXN][MAXN];
void Gauss(int n) {
	for(int i = 0;i < n;i ++) {
		for(int j = i+1;j < n;j ++) {
			while(a[j][i]) {
				int nm = a[i][i] / a[j][i];
				for(int k = i;k <= n;k ++) {
					(a[i][k] += MOD - a[j][k]*1ll*nm % MOD);
					if(a[i][k] >= MOD) a[i][k] -= MOD;
					swap(a[i][k],a[j][k]);
				}
			}
		}
	}
	for(int i = n-1;i >= 0;i --) {
		for(int j = i+1;j < n;j ++) {
			(a[i][n] += MOD - a[i][j]*1ll*a[j][n]%MOD) %= MOD;
			a[i][j] = 0;
		}
		a[i][n] = a[i][n] *1ll* qkpow(a[i][i],MOD-2) % MOD;
		a[i][i] = 1;
	} return ;
}
int main() {
	freopen("huawei.in","r",stdin);
	freopen("huawei.out","w",stdout);
	n = read();m = read();
	int P = read(),Q = read(),pq = qkpow((P+Q)%MOD,MOD-2);
	p = P*1ll*pq % MOD; q = Q*1ll*pq % MOD;
	A.n = A.m = B.n = B.m = 2;
	A.s[0][0][0] = A.s[1][1][0] = 1;
	B.s[1][1][0] = 1;
	swap(p,q);
	int pi = qkpow((MOD+1-p)%MOD,MOD-2),pw = qkpow(p,m),pwi = qkpow((MOD+1-pw)%MOD,MOD-2);
	B.s[1][0][0] = 1;
	for(int j = 0,po = pwi*1ll*q%MOD;j < m;j ++,po = po*1ll*p % MOD) {
		B.s[0][0][j] = po;
	}
	int nn = n-1;
	for(int i = 0;i <= 30 && nn > 0;i ++) {
		b[i] = B;
		if(nn & 1) {
			A = A * B;
		}
		B = B * B; nn >>= 1;
	}
	mat tl; tl.n = 1; tl.m = 2;
	for(int i = 0;i < m;i ++) tl.s[0][1][i] = pi;
	for(int i = 0;i < m;i ++) {
		a[i][i] = 1;
		if(i == 0) continue;
		mat t2 = tl * A;
		a[i][m] = (t2.s[0][0][i]*1ll*q % MOD + 1) % MOD;
		a[i][i-1] = MOD - p;
		for(int j = 0;j < m;j ++) {
			(a[i][j] += MOD - A.s[0][0][(i+m-j)%m]*1ll*q % MOD) %= MOD;
		}
	}
	Gauss(m);
	A = tl;
	for(int i = 0;i < m;i ++) A.s[0][0][i] = a[i][m];
	int T = read();
	while(T --) {
		s = read();o = read();
		mat m0 = A;
		for(int i = 0;i <= 30;i ++) {
			if(s & 1) m0 = m0 * b[i];
			s >>= 1;
		}
		AIput(m0.s[0][0][o],'\n');
	}
	return 0;
}