[ÑÖÏ Simulation Competition] fading (matrix acceleration, cyclic convolution, Gauss elimination)

Topic

4s ,512mb

Answer key

set up $dp[x][y]$ Indicates arrival $x,y$ The expected time ,$dp[0][0]=0$, Easy to get transfer
$$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$$

But because of $x$ Too big , and $y<400$ , So this is not the formula we want , We need fewer principal components to facilitate Gaussian elimination . We have to put $dp[x][y-1]$ an , Expand to infinity , No, $x$ Item ：
$$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}$$

It's a question from the beginning $\{\{dp[x-1][y]\},1\}$ To $\{\{dp[x][y]\},1\}$ Linear transformation of , So we can use matrix multiplication . Calculate the transfer matrix $n-1$ After power , set $dp[0][0\sim m-1]$ The main element , We can get it $m$ The coefficients of the equations , Then the principal component is obtained by Gauss elimination .

In the process, the transfer matrix is preprocessed $1,2,4,...2^{\log n}$ Power , When asking, double the number $dp[X_i][Y_i]$, Matrix multiplication can be obtained from $O(m^3)$ Turn into $O(m^2)$ , The old routine .

thus , The time complexity is $O(m^3\log n+t{m}^2\log n)$ , Just can't pass .

Look back at the formula , We found that $\frac{p{q}^i}{1-q^m}$ Follow $y$ It doesn't matter. , in other words , This is a cyclic convolution .

Make $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$ , that
$$f_k=f_{k-1}\ast g+h$$

This is a linear transformation of polynomials , We pack 2×2 Matrix multiplication with polynomial matrix of , use NTT Convolution simulates cyclic convolution , The complexity of a transformation starts from $O(m^3)$ Turned into $O(m\log m)$ .

The final complexity is $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;
}