【UOJ574】多线程计算【二元二项式反演】【定积分】【矩阵】【NTT 卷积】

题意：有 $n×mn\times m$ 的网格，每个结点在 $[0, 1)$ 内的一个随机时刻被点亮。有 $h$ 个数对 $x_i,y_i$ ，对于一个瞬间状态，如果存在一个 $x_i,y_i$ 使得恰好有 $x_i$ 行 $y_i$ 列被点亮，则每单位时间产生值为斐波拉契第 $k$ 项的贡献， $k$ 为点亮的结点个数。求期望贡献之和模 $998244353$ 。

$1≤h≤n×m≤5×1051\leq h\leq n\times m\leq 5\times 10^5$

恰好不好搞，先容斥。显然一个时刻最多一个条件满足，所以求所有条件的贡献和就可以了。（其实多个条件是给部分分用的）

设 $f (x, y)$ 表示钦定 $x$ 行 $y$ 列全亮的贡献之和， $g (x, y)$ 表示恰好有 $x$ 行 $y$ 列全亮的贡献。

有

$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)$

由二项式反演的推广可知

$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)$

文末有这个东西的证明。

现在的问题是怎么求 $f (i, j)$ 。因为钦定必须亮的行列数给定后，可以计算出必须亮的结点的个数为 $n m - (n - i) (m - j)$ ，于是转换为下面这个问题：

有 $x + y$ 个结点，其中前 $x$ 个必须亮，后 $y$ 个随意，求总贡献。

（这个 $x, y$ 和前面的 $x, y$ 没有关系）

注意要乘个组合数。

这个贡献不是离散的，我们先要找点结论出来。考虑一个时刻 $t$ ，所有结点亮的概率都是独立的 $t$ 。那个斐波拉契实际上就是矩阵的幂，设转移矩阵为 $M$ 。

先用 $y = 0$ 的找点灵感，产生的贡献就是 $M^xt^x$ 。

求个积分

$∫01Mxtxdt=1x+1Mx\int_0^1 M^xt^x dt=\frac{1}{x+1}M^x$

然后就可以用类似的思路推广出所有多项式的情况。

开始想的再套一个容斥上去，但式子十分精神污染。然后想的积分的时候强制后面的不选，然后推出了一模一样的式子……

考虑这些式子本质肯定都相同，你推不出来是数学不好。所以考虑运用一些成熟的结论多跳几步。

注意到后面任意选里面带的是组合数，可以考虑直接用二项式定理，然后直接把矩阵积分。即

$∫01(tM)x(tM+1−t)ydt\int_0^1(tM)^x(tM+1-t)^y dt$

右边直接拆不好搞，但注意到可以把 $t$ 合并

$∫01(tM)x(t(M−1)+1)ydt\int_0^1(tM)^x(t(M-1)+1)^y dt$

直接考察这个式子的 $x + i$ 项系数，它等于 $Mx(yi)(M−1)iM^x\binom yi (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)^i$

暴力做这个东西，然后前面反演式也暴力算，可以做到 $O((nm)2)\Omicron((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)!}$

找不到更好的姿势，就直接糊了个矩阵 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)$

$(−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)!}$

两维是独立的，所以横着对每一行卷一次，再竖着对每一列卷一次就可以了。

于是又写成了代数大杂烩……

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <vector>
#define MAXN (1<<20)+5
using 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=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=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=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)=∑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)$