一个 (R imes C) 的棋盘，你有 (Q) 组询问，每次询问国王走 (R-1) 步从 ((1,x)) 到达 ((R,y)) 有多少种方案。你只需要输出答案对 (998244353) 取模的结果。
首先不难有一个每组询问(O(C^3 log R))的做法，即设转移矩阵 (T) 满足 (T[i][j] = [ |i-j| leq 1 ],)然后矩阵快速幂。
因为转移是向量乘矩阵的形式，所以在(C)固定的情况下，对于任意的(x,y,)答案都是一个关于(R)的线性递推，并且有相同的特征多项式。
(T)的特征多项式(f(lambda))为(det(lambda I - T),)即我们需要求如下矩阵的行列式:
(egin{equation*} egin{bmatrix} lambda - 1 & -1 & & \ -1 & lambda - 1 & -1 & \ & -1 & lambda - 1 & ddots \ & & -1 & ddots \ end{bmatrix} end{equation*})
不难发现行列式/特征多项式满足递推式 : (F_C = (lambda - 1)F_{C-1} - F_{C-2})
可以通过矩阵乘法在(Theta(log C))的时间内求出(F_C)的一个点值，那么求出所有单位根的点值然后再(IDFT,)就可以在(Theta(Clog C))的时间内获得(F_C)的系数了。
记(x^R)对特征多项式(F_C)取模之后的每一项系数为(r_0,r_1,..r_{C-1},)那么(ans_R = sumlimits_{i=0}^{C-1} ans_i imes r_i.)
考虑翻折容斥，因为(i<C,)所以上下界不会同时越界，因此可以把有边界的(ans)转化成(3)个无边界的查询。
对于无边界的情况，我们要求(sumlimits_{i=0}^{C-1} r_i (x^{-1} + x^0 + x^1)^i,)分治(FFT)即可。
然后就可以(Theta(1))查询了。
(Theta(Clog C log R) - Theta(1).)
code :
#include <bits/stdc++.h>
#define LL unsigned long long
#define uint unsigned int
using namespace std;
template <typename T> void read(T &x){
static char ch; x = 0,ch = getchar();
while (!isdigit(ch)) ch = getchar();
while (isdigit(ch)) x = x * 10 + ch - '0',ch = getchar();
}
inline void write(int x){ if (x > 9) write(x/10); putchar(x%10+'0'); }
const uint P = 998244353,N = 262144;
inline uint power(uint x,int y){
static uint r; r = 1; while (y){ if (y&1) r = (LL)r * x % P; x = (LL)x * x % P,y >>= 1; } return r;
}
uint rt[21],irt[21],R[N<<1],inv[N<<1|1],wn[N<<2],iwn[N<<2];
inline int getR(int n){
static int l,L; l = 0,L = 1; while (L <= n) L<<=1,++l;
for (register int i = 1; i < L; ++i) R[i] = (R[i>>1]>>1) | ((i&1)<<l-1);
return L;
}
inline void NTT(uint *A,int n){
register int i,j,k; uint v;
for (i = 1; i < n; ++i) if (i < R[i]) swap(A[i],A[R[i]]);
for (i = 1; i < n; i<<=1) for (j = 0; j < n; j += i<<1) for (k = j; k < i+j; ++k)
v = (LL)wn[(i<<1)+k-j] * A[k+i] % P,A[k+i] = (A[k]<v)?(A[k]+P-v):(A[k]-v),A[k] = (A[k]+v>=P)?(A[k]+v-P):(A[k]+v);
}
inline void iNTT(uint *A,int n){
register int i,j,k; uint v;
for (i = 1; i < n; ++i) if (i < R[i]) swap(A[i],A[R[i]]);
for (i = 1; i < n; i<<=1) for (j = 0; j < n; j += i<<1) for (k = j; k < i+j; ++k)
v = (LL)iwn[(i<<1)+k-j] * A[k+i] % P,A[k+i] = (A[k]<v)?(A[k]+P-v):(A[k]-v),A[k] = (A[k]+v>=P)?(A[k]+v-P):(A[k]+v);
for (i = 0,v = inv[n]; i < n; ++i) A[i] = (LL)A[i] * v % P;
}
inline void Inv(uint *A,uint *B,int n){
static int i,l,L; static uint _A[N];
L = 1; while (L <= n) L <<= 1; n = L; B[0] = power(A[0],P-2);
for (L = l = 1; L <= n; L <<= 1,++l){
for (i = 0; i < L<<1; ++i) R[i] = (R[i>>1]>>1) | ((i&1)<<l-1);
memcpy(_A,A,L<<2),memset(_A+L,0,L<<2),memset(B+L,0,L<<2);
NTT(_A,L<<1),NTT(B,L<<1);
for (i = 0; i < L<<1; ++i) B[i] = (LL)B[i] * (2 - (LL)_A[i] * B[i] % P + P) % P;
iNTT(B,L<<1),memset(B+L,0,L<<2);
}
}
int K,C;
namespace Get_Characteristic_polynomial{
struct Mat{
LL a00,a01,a10,a11;
Mat(int v00=0,int v01=0,int v10=0,int v11=0){ a00=v00,a01=v01,a10=v10,a11=v11; }
}I(1,0,0,1);
inline Mat operator * (Mat A,Mat B){
return Mat((A.a00 * B.a00 + A.a01 * B.a10) % P,(A.a00 * B.a01 + A.a01 * B.a11) % P,
(A.a10 * B.a00 + A.a11 * B.a10) % P,(A.a10 * B.a01 + A.a11 * B.a11) % P);
}
inline Mat power(Mat A,int n){
static Mat T; T = I; while (n){ if (n&1) T = T * A; A = A * A,n >>= 1; } return T;
}
inline uint power(uint x,int y){
static uint r; r = 1; while (y){ if (y&1) r = (LL)r * x % P; x = (LL)x * x % P,y >>= 1; } return r;
}
inline uint calc(uint w){
static Mat T; --w,T = power(Mat(w,P-1,1,0),C);
return (T.a10 * w + T.a11) % P;
}
inline void solve(uint *G){
int L = getR(C);
for (register int i = 0; i < L; ++i) G[i] = calc(wn[L+i]);
iNTT(G,L);
}
}
uint r[N];
namespace Linear_sequence{
int L; uint F[N],G[N],Gr[N];
inline void poly_mod(uint *A){
register int i; int n = C<<1,m = C; while (n && !A[n]) --n; if (n < m) return;
reverse(A,A+n+1); memcpy(F,A,C<<2),memset(F+C,0,L-C<<2);
NTT(F,L); for (i = 0; i < L; ++i) F[i] = (LL)F[i] * Gr[i] % P; iNTT(F,L);
memset(F+n-m+1,0,L-n+m-1<<2); reverse(F,F+n-m+1);
NTT(F,L); for (i = 0; i < L; ++i) F[i] = (LL)F[i] * G[i] % P; iNTT(F,L);
reverse(A,A+n+1); for (i = 0; i <= n; ++i) A[i] = (A[i]<F[i])?(A[i]+P-F[i]):(A[i]-F[i]);
}
uint A[N],B[N];
inline void solve(){
register int i; int k = K-1; L = getR(C<<1);
reverse(G,G+C+1),Inv(G,Gr,C),NTT(Gr,L),reverse(G,G+C+1),NTT(G,L);
A[0] = 1,B[1] = 1,NTT(B,L);
while (k){
if (k&1){
NTT(A,L);
for (i = 0; i < L; ++i) A[i] = (LL)A[i] * B[i] % P;
iNTT(A,L),poly_mod(A);
}
if (k>>=1){
for (i = 0; i < L; ++i) B[i] = (LL)B[i] * B[i] % P;
iNTT(B,L),poly_mod(B),NTT(B,L);
}
}
memcpy(r,A,C<<2);
}
}
struct poly{
vector<uint>F; int z;
inline void free(){ vector<uint>().swap(F); }
inline void clear(){ for (int i = 0; i < F.size(); ++i) F[i] = 0; }
inline void add(int id,uint v){ F[z+id] = (F[z+id]+v >= P) ? (F[z+id]+v-P) : (F[z+id]+v); }
inline uint ask(int id){ return (z+id < 0 || z+id >= F.size()) ? 0 : F[z+id]; }
};
namespace polyopt{
uint F[N<<1],G[N<<1];
int n,m,Lim;
inline void Mul(poly &A,poly &B,poly &C){ // C = A * B
register int i;
C.z = A.z + B.z; n = A.F.size()-1; m = B.F.size()-1; Lim = getR(n+m);
memset(F+n+1,0,Lim-n-1<<2); for (i = 0; i <= n; ++i) F[i] = A.F[i];
memset(G+m+1,0,Lim-m-1<<2); for (i = 0; i <= m; ++i) G[i] = B.F[i];
NTT(F,Lim),NTT(G,Lim); for (i = 0; i < Lim; ++i) F[i] = (LL)F[i] * G[i] % P; iNTT(F,Lim);
C.F.resize(n+m+1); for (i = 0; i <= n+m; ++i) C.F[i] = F[i];
}
inline void Mul2(poly &A,poly &B){ // B = A * A
register int i;
B.z = A.z<<1; n = A.F.size()-1; Lim = getR(n<<1);
memset(F+n+1,0,Lim-n-1<<2); for (i = 0; i <= n; ++i) F[i] = A.F[i];
NTT(F,Lim); for (i = 0; i < Lim; ++i) F[i] = (LL)F[i] * F[i] % P; iNTT(F,Lim);
B.F.resize(n<<1|1); for (i = 0; i <= n<<1; ++i) B.F[i] = F[i];
}
inline void Add(poly &A,poly &B,poly &C){ //C = A + B
register int i;
C.free(),C.z = max(A.z,B.z),C.F.resize(C.z+max(A.F.size()-A.z,B.F.size()-B.z)),C.clear();
for (i = 0; i < A.F.size(); ++i) C.add(i-A.z,A.F[i]);
for (i = 0; i < B.F.size(); ++i) C.add(i-B.z,B.F[i]);
}
}
using polyopt::Mul;
using polyopt::Mul2;
using polyopt::Add;
namespace DivideFFT{
poly tmp[18],F; int L,id[N];
inline void divide(poly &T,int L,int R){
if (L == R){ T.F.push_back(r[L]),T.z = 0; return; }
int mid = L+R>>1; poly Le,Ri; divide(Le,L,mid); divide(Ri,mid+1,R);
Mul(Ri,tmp[id[mid-L+1]],Ri),Add(Le,Ri,T),Ri.free(),Le.free();
// Ri.F.clear(),Le.F.clear(),Ri.F.shrink_to_fit(),Le.F.shrink_to_fit();
// vector<int>().swap(Ri.F),vector<int>().swap(Le.F);
}
inline void solve(){
int i,j; L = 1; while (L < C) L <<= 1;
tmp[0].F.resize(3),tmp[0].F[0] = tmp[0].F[1] = tmp[0].F[2] = tmp[0].z = 1;
for (i = 1,j = 2; j <= L; ++i,j<<=1) id[j] = i,Mul2(tmp[i-1],tmp[i]);
divide(F,1,L),Mul(F,tmp[0],F),F.add(0,r[0]);
}
inline int query(int x,int y){
return ((F.ask(x-y) + P - F.ask(x+y)) % P + P - F.ask(C*2+2-x-y)) % P;
}
}
int main(){
register int i,j,t; uint v;
for (i = 1,j = 2; i <= 20; ++i,j <<= 1) rt[i] = power(3,(P-1)/j),irt[i] = power(rt[i],P-2);
for (i = 1,j = 2; i <= 20 && j <= N<<1; ++i,j <<= 1){
for (wn[j] = 1,t=j+1,v=rt[i]; t < j<<1; ++t) wn[t] = (LL)wn[t-1] * v % P;
for (iwn[j] = 1,t=j+1,v=irt[i]; t < j<<1; ++t) iwn[t] = (LL)iwn[t-1] * v % P;
}
for (inv[0] = inv[1] = 1,i = 2; i <= N<<1; ++i) inv[i] = (LL)(P-P/i) * inv[P%i] % P;
read(K),read(C); // K = 1000000000,C = 100000;
Get_Characteristic_polynomial::solve(Linear_sequence::G); // cerr << clock() << '
';
Linear_sequence::solve(); // cerr << clock() << '
';
DivideFFT::solve(); // cerr << clock() << '
';
int q,x,y; read(q);
while (q--) read(x),read(y),write(DivideFFT::query(x,y)),putchar('
');
return 0;
}