题意：有 n×mn\times mn×m 的网格，每个结点在 [0,1)[0,1)[0,1) 内的一个随机时刻被点亮。有 hhh 个数对 xi,yix_i,y_ixi,yi，对于一个瞬间状态，如果存在一个 xi,yix_i,y_ixi,yi 使得恰好有 xix_ixi 行 yiy_iyi 列被点亮，则每单位时间产生值为斐波拉契第 kkk 项的贡献，kkk 为点亮的结点个数。求期望贡献之和 模 998244353998244353998244353。
1≤h≤n×m≤5×1051\leq h\leq n\times m\leq 5\times 10^51≤h≤n×m≤5×105
恰好不好搞，先容斥。显然一个时刻最多一个条件满足，所以求所有条件的贡献和就可以了。（其实多个条件是给部分分用的）
设 f(x,y)f(x,y)f(x,y) 表示钦定 xxx 行 yyy 列全亮的贡献之和，g(x,y)g(x,y)g(x,y) 表示恰好有 xxx 行 yyy 列全亮的贡献。
有
f(x,y)=∑i=xn∑j=ym(ix)(jy)g(i,j)f(x,y)=\sum_{i=x}^n\sum_{j=y}^{m}\binom ix\binom jyg(i,j)f(x,y)=i=x∑nj=y∑m(xi)(yj)g(i,j)
由二项式反演的推广可知
g(x,y)=∑i=xn∑j=ym(−1)i+j−x−y(ix)(jy)f(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m(-1)^{i+j-x-y}\binom ix\binom jyf(i,j)g(x,y)=i=x∑nj=y∑m(−1)i+j−x−y(xi)(yj)f(i,j)
文末有这个东西的证明。
现在的问题是怎么求 f(i,j)f(i,j)f(i,j)。因为钦定必须亮的行列数给定后，可以计算出必须亮的结点的个数为 ，于是转换为下面这个问题：
有 x+yx+yx+y 个结点，其中前 xxx 个必须亮，后 yyy 个随意，求总贡献。（这个 x,yx,yx,y 和前面的 x,yx,yx,y 没有关系）
注意要乘个组合数。
这个贡献不是离散的，我们先要找点结论出来。考虑一个时刻 ttt，所有结点亮的概率都是独立的 ttt。那个斐波拉契实际上就是矩阵的幂，设转移矩阵为 MMM。
先用 y=0y=0y=0 的找点灵感，产生的贡献就是 MxtxM^xt^xMxtx。
求个积分
∫01Mxtxdt=1x+1Mx\int_0^1 M^xt^x dt=\frac{1}{x+1}M^x∫01Mxtxdt=x+11Mx
然后就可以用类似的思路推广出所有多项式的情况。
开始想的再套一个容斥上去，但式子十分精神污染。然后想的积分的时候强制后面的不选，然后推出了一模一样的式子……
考虑这些式子本质肯定都相同，你推不出来是数学不好。所以考虑运用一些成熟的结论多跳几步。
注意到后面任意选里面带的是组合数，可以考虑直接用二项式定理，然后直接把矩阵积分。即
∫01(tM)x(tM+1−t)ydt\int_0^1(tM)^x(tM+1-t)^y dt∫01(tM)x(tM+1−t)ydt
右边直接拆不好搞，但注意到可以把 ttt 合并
∫01(tM)x(t(M−1)+1)ydt\int_0^1(tM)^x(t(M-1)+1)^y dt∫01(tM)x(t(M−1)+1)ydt
直接考察这个式子的 x+ix+ix+i 项系数，它等于 Mx(yi)(M−1)iM^x\binom yi (M-1)^iMx(iy)(M−1)i
所以积分后的总贡献就是
∑i=0y1x+i+1Mx(yi)(M−1)i\sum_{i=0}^y \frac 1{x+i+1}M^x\binom yi (M-1)^ii=0∑yx+i+11Mx(iy)(M−1)i
暴力做这个东西，然后前面反演式也暴力算，可以做到 O((nm)2)\Omicron((nm)^2)O((nm)2)
注意到这个式子是卷积的形式
f(k)=∑i=0k1n−k+i+1Mn−kk!i!(k−i)!(M−1)i=Mn−kk!∑i=0k(M−1)ii!1(n−k+i+1)(k−i)!f(k)=\sum_{i=0}^k\frac 1{n-k+i+1}M^{n-k}\frac{k!}{i!(k-i)!}(M-1)^i\\=M^{n-k}k!\sum_{i=0}^k\frac{(M-1)^i}{i!}\frac 1{(n-k+i+1)(k-i)!}f(k)=i=0∑kn−k+i+11Mn−ki!(k−i)!k!(M−1)i=Mn−kk!i=0∑ki!(M−1)i
找不到更好的姿势，就直接糊了个矩阵 NTT 上去，荣获赛时最慢代码。
后面的反演式
g(x,y)=∑i=xn∑j=ym(−1)i+j−x−y(ix)(jy)f(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m(-1)^{i+j-x-y}\binom ix\binom jyf(i,j)g(x,y)=i=x∑nj=y∑m(−1)i+j−x−y(xi)(yj)f(i,j)
(−1)x+yx!y!g(x,y)=∑i=xn∑j=ym(−1)i+ji!j!f(i,j)(i−x)!(j−y)!(-1)^{x+y}x!y!g(x,y)=\sum_{i=x}^n\sum_{j=y}^m\frac{(-1)^{i+j}i!j!f(i,j)}{(i-x)!(j-y)!}(−1)x+yx!y!g(x,y)=i=x∑nj=y∑m(i−x)!(j−y)!(−1)i+ji!j!f(i,j)
两维是独立的，所以横着对每一行卷一次，再竖着对每一列卷一次就可以了。
于是又写成了代数大杂烩……
#include <iostream>#include <cstdio>#include <cstring>#include <cctype>#include <vector>#define MAXN (1<<20)+5using namespace std;const int MOD=998244353;typedef long long ll;inline int add(const int& x,const int& y){return x+y>=MOD? x+y-MOD:x+y;}inline int dec(const int& x,const int& y){return x<y? x-y+MOD:x-y;}inline int qpow(int a,int p){int ans=1;while (p){if (p&1) ans=(ll)ans*a%MOD;a=(ll)a*a%MOD,p>>=1;}return ans;}#define inv(x) qpow(x,MOD-2)inline int read(){int ans=0;char c=getchar();while (!isdigit(c)) c=getchar();while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();return ans;}int fac[MAXN],finv[MAXN];struct mat{int e[2][2];inline mat(const int& v=0){e[0][0]=e[1][1]=v;e[0][1]=e[1][0]=0;}inline int* operator [](const int& i){return e[i];}inline const int* operator [](const int& i)const{return e[i];}};inline mat operator +(const mat& a,const mat& b){mat c;for (int i=0;i<2;i++)for (int j=0;j<2;j++)c[i][j]=add(a[i][j],b[i][j]);return c;}inline mat operator -(const mat& a,const mat& b){mat c;for (int i=0;i<2;i++)for (int j=0;j<2;j++)c[i][j]=dec(a[i][j],b[i][j]);return c;}inline mat operator *(mat a,const int& b){for (int i=0;i<2;i++)for (int j=0;j<2;j++)a[i][j]=(ll)a[i][j]*b%MOD;return a;}inline mat operator *(const mat& a,const mat& b){mat c;for (int i=0;i<2;i++)for (int k=0;k<2;k++)for (int j=0;j<2;j++)c[i][j]=(c[i][j]+(ll)a[i][k]*b[k][j])%MOD;return c;}mat P[MAXN],Q[MAXN];vector<int> f[MAXN],g[MAXN];inline int C(const int& n,const int& m){return (ll)fac[n]*finv[m]%MOD*finv[n-m]%MOD;}int l,lim,r[MAXN];inline void init(){lim=1<<l;for (int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));}int rt[2][24];void ntt(mat* a,int type){for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1;int Wn=rt[type][L+1];for (int s=0;s<lim;s+=len)for (int k=0,w=1;k<mid;k++,w=(ll)w*Wn%MOD){mat x=a[s+k],y=a[s+mid+k]*w;a[s+k]=x+y,a[s+mid+k]=x-y;}}if (type){int t=inv(lim);for (int i=0;i<lim;i++) a[i]=a[i]*t;}}void ntt(int* a,int type){for (int i=0;i<lim;i++) if (i<r[i]) swap(a[i],a[r[i]]);for (int L=0;L<l;L++){int mid=1<<L,len=mid<<1;int Wn=rt[type][L+1];for (int s=0;s<lim;s+=len)for (int k=0,w=1;k<mid;k++,w=(ll)w*Wn%MOD){int x=a[s+k],y=(ll)a[s+mid+k]*w%MOD;a[s+k]=add(x,y),a[s+mid+k]=dec(x,y);}}if (type){int t=inv(lim);for (int i=0;i<lim;i++) a[i]=(ll)a[i]*t%MOD;}}mat F[MAXN],res[MAXN];int A[MAXN],B[MAXN];inline int calc(int x,int y){mat ans=P[x]*res[y]*fac[y];return add(ans[1][0],ans[1][1]);}int main(){rt[0][23]=qpow(3,119),rt[1][23]=inv(rt[0][23]);for (int i=22;i>=0;i--){rt[0][i]=(ll)rt[0][i+1]*rt[0][i+1]%MOD;rt[1][i]=(ll)rt[1][i+1]*rt[1][i+1]%MOD;}int n,m;n=read(),m=read();fac[0]=1;for (int i=1;i<=n*m;i++) fac[i]=(ll)fac[i-1]*i%MOD;finv[n*m]=qpow(fac[n*m],MOD-2);for (int i=n*m-1;i>=0;i--) finv[i]=(ll)finv[i+1]*(i+1)%MOD;P[0]=Q[0]=mat(1);mat p,q;p[0][0]=p[0][1]=p[1][0]=1;q[0][1]=q[1][0]=1,q[1][1]=MOD-1;for (int i=1;i<=n*m;i++) P[i]=P[i-1]*p,Q[i]=Q[i-1]*q;for (int i=0;i<=n;i++) f[i].resize(m+1),g[i].resize(m+1);for (l=0;(1<<l)<=2*n*m;l++);init();for (int i=0;i<=n*m;i++) F[i]=Q[i]*finv[i],A[i]=(ll)finv[i]*inv(n*m+1-i)%MOD;ntt(F,0),ntt(A,0);for (int i=0;i<lim;i++) res[i]=F[i]*A[i];ntt(res,1);for (int i=0;i<=n;i++)for (int j=0;j<=m;j++)f[i][j]=(ll)C(n,i)*C(m,j)%MOD*calc(i*m+j*n-i*j,(n-i)*(m-j))%MOD;//for (int i=0;i<=n;i++,puts(""))//for (int j=0;j<=m;j++)//printf("%d ",f[i][j]);//for (int i=n;i>=0;i--)//for (int j=m;j>=0;j--)//{//g[i][j]=f[i][j];//if (i<n) g[i][j]=add(g[i][j],g[i+1][j]);//if (j<m) g[i][j]=add(g[i][j],g[i][j+1]);//if (i<n&&j<m) g[i][j]=dec(g[i][j],g[i+1][j+1]);//}for (int i=0;i<=n;i++)for (int j=0;j<=m;j++){f[i][j]=(ll)f[i][j]*fac[i]%MOD*fac[j]%MOD;if ((i+j)&1) f[i][j]=dec(0,f[i][j]);}for (l=0;(1<<l)<=(m<<1);l++);init();for (int T=0;T<=n;T++){for (int i=0;i<lim;i++) A[i]=B[i]=0;for (int i=0;i<=m;i++) A[i]=f[T][i],B[i]=finv[m-i];ntt(A,0),ntt(B,0);for (int i=0;i<lim;i++) A[i]=(ll)A[i]*B[i]%MOD;ntt(A,1);for (int i=m;i<=2*m;i++) f[T][i-m]=A[i];}for (l=0;(1<<l)<=(n<<1);l++);init();for (int T=0;T<=m;T++){for (int i=0;i<lim;i++) A[i]=B[i]=0;for (int i=0;i<=n;i++) A[i]=f[i][T],B[i]=finv[n-i];ntt(A,0),ntt(B,0);for (int i=0;i<lim;i++) A[i]=(ll)A[i]*B[i]%MOD;ntt(A,1);for (int i=n;i<=2*n;i++) f[i-n][T]=A[i];}for (int i=0;i<=n;i++)for (int j=0;j<=m;j++){f[i][j]=(ll)f[i][j]*finv[i]%MOD*finv[j]%MOD;if ((i+j)&1) f[i][j]=dec(0,f[i][j]);}int ans=0;for (int T=read();T;T--){int x,y;x=read(),y=read();ans=add(ans,f[x][y]);}cout<<ans;return 0;} 附：二元二项式反演证明g(x,y)=∑i=xn∑j=ym(ix)(jy)[i−x=0][j−y=0]g(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m\binom ix\binom jy[i-x=0][j-y=0]g(i,j)g(x,y)=i=x∑nj=y∑m(xi)(yj)[i−x=0][j−y=0]g(i,j)
g(x,y)=∑i=xn∑j=ym(ix)(jy)∑a=0i−x∑b=0j−y(−1)a+b(i−xa)(j−yb)g(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m\binom ix\binom jy\sum_{a=0}^{i-x}\sum_{b=0}^{j-y}(-1)^{a+b}\binom{i-x}a\binom{j-y}bg(i,j)g(x,y)=i=x∑nj=y∑m(xi)(yj)a=0∑i−xb=0∑j−y(−1)a+b(ai−x)(bj−y)g(i,j)
g(x,y)=∑i=xn∑j=ym∑a=0i−x∑b=0j−y(−1)a+b(ix+a)(x+ax)(jy+b)(y+by)g(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m\sum_{a=0}^{i-x}\sum_{b=0}^{j-y}(-1)^{a+b}\binom i{x+a}\binom{x+a}{x}\binom{j}{y+b}\binom{y+b}{y}g(i,j)g(x,y)=i=x∑nj=y∑ma=0∑i−xb=0∑j−y(−1)a+b(x+ai)(xx+a)(y+bj)(yy+b)g(i,j)
g(x,y)=∑i=xn∑j=ym∑a=xi∑b=yj(−1)a+b−x−y(ia)(ax)(jb)(by)g(i,j)g(x,y)=\sum_{i=x}^n\sum_{j=y}^m\sum_{a=x}^{i}\sum_{b=y}^{j}(-1)^{a+b-x-y}\binom i{a}\binom{a}{x}\binom{j}{b}\binom{b}{y}g(i,j)g(x,y)=i=x∑nj=y∑ma=x∑ib=y∑j(−1)a+b−x−y(ai)(xa)(bj)(yb)g(i,j)
g(x,y)=∑a=xn∑b=ym(−1)a+b−x−y(ax)(by)∑i=xn∑j=ym(ix)(jy)g(i,j)g(x,y)=\sum_{a=x}^{n}\sum_{b=y}^{m}(-1)^{a+b-x-y}\binom{a}{x}\binom{b}{y}\sum_{i=x}^n\sum_{j=y}^m\binom ix \binom jyg(i,j)g(x,y)=a=x∑nb=y∑m(−1)a+b−x−y(xa)(yb)i=x∑nj=y∑m(xi)(yj)g(i,j)
g(x,y)=∑i=xn∑j=ym(−1)i+j−x−y(ix)(jy)f(i,j)g(x,y)=\sum_{i=x}^{n}\sum_{j=y}^{m}(-1)^{i+j-x-y}\binom{i}{x}\binom{j}{y}f(i,j)g(x,y)=i=x∑nj=y∑m(−1)i+j−x−y(xi)(yj)f(i,j)
