題意：有 $n\times m$ 的網格，每個結點在 $[0,1)$ 內的一個隨機時刻被點亮。有 $h$ 個數對 $x_i,y_i$，對于一個瞬間狀態，如果存在一個 $x_i,y_i$ 使得恰好有 $x_i$ 行 $y_i$ 列被點亮，則每單位時間產生值為斐波拉契第 $k$ 項的貢獻，$k$ 為點亮的結點個數。求期望貢獻之和 模 $998244353$。

$1\le h\le n\times m\le 5\times 10^5$

恰好不好搞，先容斥。顯然一個時刻最多一個條件滿足，所以求所有條件的貢獻和就可以了。（其實多個條件是給部分分用的）

設 $f(x,y)$ 表示欽定 $x$ 行 $y$ 列全亮的貢獻之和，$g(x,y)$ 表示恰好有 $x$ 行 $y$ 列全亮的貢獻。

有

$$f(x,y)=\sum _{i=x}^n\sum _{j=y}^m\left(\genfrac{}{}{0ex}{}{i}{x}\right)\left(\genfrac{}{}{0ex}{}{j}{y}\right)g(i,j)$$

由二項式反演的推廣可知

$$g(x,y)=\sum _{i=x}^n\sum _{j=y}^m(?1{)}^{i+j?x?y}(\genfrac{}{}{0ex}{}{i}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{y})f(i,j)$$

文末有這個東西的證明。

現在的問題是怎么求 $f(i,j)$。因為欽定必須亮的行列數給定后，可以計算出必須亮的結點的個數為 $nm?(n?i)(m?j)$，于是轉換為下面這個問題：

• 有 $x+y$ 個結點，其中前 $x$ 個必須亮，后 $y$ 個隨意，求總貢獻。

（這個 $x,y$ 和前面的 $x,y$ 沒有關系）

注意要乘個組合數。

這個貢獻不是離散的，我們先要找點結論出來。考慮一個時刻 $t$，所有結點亮的概率都是獨立的 $t$。那個斐波拉契實際上就是矩陣的冪，設轉移矩陣為 $M$。

先用 $y=0$ 的找點靈感，產生的貢獻就是 $M^x{t}^x$。

求個積分

$${\int }_0^1{M}^x{t}^{x}dt=\frac{1}{x+1}{M}^x$$

然后就可以用類似的思路推廣出所有多項式的情況。

開始想的再套一個容斥上去，但式子十分精神污染。然后想的積分的時候強制后面的不選，然后推出了一模一樣的式子……

考慮這些式子本質肯定都相同，你推不出來是數學不好。所以考慮運用一些成熟的結論多跳幾步。

注意到后面任意選里面帶的是組合數，可以考慮直接用二項式定理，然后直接把矩陣積分。即

$${\int }_0^1(tM{)}^x(tM+1?t{)}^{y}dt$$

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

$${\int }_0^1(tM{)}^x(t(M?1)+1{)}^{y}dt$$

直接考察這個式子的 $x+i$ 項系數，它等于 $M^x\left(\genfrac{}{}{0ex}{}{y}{i}\right)(M?1{)}^i$

所以積分后的總貢獻就是

$$\sum _{i=0}^y\frac{1}{x+i+1}{M}^x\left(\genfrac{}{}{0ex}{}{y}{i}\right)(M?1{)}^i$$

暴力做這個東西，然后前面反演式也暴力算，可以做到 $\mathrm{O}((nm{)}^2)$

注意到這個式子是卷積的形式

$$\begin{gathered} 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)!} \end{gathered}$$

找不到更好的姿勢，就直接糊了個矩陣 NTT 上去，榮獲賽時最慢代碼。

后面的反演式

$$g(x,y)=\sum _{i=x}^n\sum _{j=y}^m(?1{)}^{i+j?x?y}(\genfrac{}{}{0ex}{}{i}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{y})f(i,j)$$

$$(?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)=\sum _{i=x}^n\sum _{j=y}^m\left(\genfrac{}{}{0ex}{}{i}{x}\right)\left(\genfrac{}{}{0ex}{}{j}{y}\right)[i?x=0][j?y=0]g(i,j)$$

$$g(x,y)=\sum _{i=x}^n\sum _{j=y}^m\left(\genfrac{}{}{0ex}{}{i}{x}\right)\left(\genfrac{}{}{0ex}{}{j}{y}\right)\sum _{a=0}^{i?x}\sum _{b=0}^{j?y}(?1{)}^{a+b}(\genfrac{}{}{0ex}{}{i?x}{a}\left)\right(\genfrac{}{}{0ex}{}{j?y}{b})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}(\genfrac{}{}{0ex}{}{i}{x+a}\left)\right(\genfrac{}{}{0ex}{}{x+a}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{y+b}\left)\right(\genfrac{}{}{0ex}{}{y+b}{y})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}(\genfrac{}{}{0ex}{}{i}{a}\left)\right(\genfrac{}{}{0ex}{}{a}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{b}\left)\right(\genfrac{}{}{0ex}{}{b}{y})g(i,j)$$

$$g(x,y)=\sum _{a=x}^n\sum _{b=y}^m(?1{)}^{a+b?x?y}(\genfrac{}{}{0ex}{}{a}{x}\left)\right(\genfrac{}{}{0ex}{}{b}{y}\left)\sum _{i=x}^n\sum _{j=y}^m\right(\genfrac{}{}{0ex}{}{i}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{y})g(i,j)$$

$$g(x,y)=\sum _{i=x}^n\sum _{j=y}^m(?1{)}^{i+j?x?y}(\genfrac{}{}{0ex}{}{i}{x}\left)\right(\genfrac{}{}{0ex}{}{j}{y})f(i,j)$$

總結

以上是生活随笔為你收集整理的【UOJ574】多线程计算【二元二项式反演】【定积分】【矩阵】【NTT 卷积】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。