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

题意：有 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)。因为钦定必须亮的行列数给定后，可以计算出必须亮的结点的个数为 nm−(n−i)(m−j)nm-(n-i)(m-j)nm−(n−i)(m−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(n−k+i+1)(k−i)!1

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

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

【HAOI2018】染色【反向二项式反演】【NTT卷积】
传送门题意:NNN个位置染MMM种颜色,恰好出现SSS次的颜色数量恰好为kkk时的愉悦度为wkw_kwk,求所有方案的愉悦度之和.对100453580910045358091004535809取模 ...
【LOJ6072】苹果树【折半搜索】【矩阵树定理】【二项式反演】
题意:有好坏两种点共 nnn 个,每个好点有权值,把这 nnn 个点连成一棵树,一个好点为有用的当且仅当它至少与一个好点相邻,求所有有用的点的权值和不超过 limlimlim 的方案数. n≤40n\ ...
matlab判断矩阵不可约,用Matlab计算二元域GF(2)上的不可约多项式
1 二元域 GF(2) 上的不可约多项式二元域 GF(2)={0,1} 上的运算规则如下: 加法:+ 0 1 0 0 1 1 1 0 乘法:⋅ 0 1 0 0 0 1 0 1 二元域 GF(2) 上 ...
ACM数论之旅17---反演定理第一回二项式反演（神说要有光于是就有了光(´・ω・`)）...
终于讲到反演定理了,反演定理这种东西记一下公式就好了,反正我是证明不出来的~(-o￣▽￣)-o 首先,著名的反演公式我先简单的写一下o(￣ヘ￣*o) 比如下面这个公式 f(n) = g(1) + g ...
二项式反演[bzoj3622]已经没有什么好害怕的了
前言继续学习容斥的技巧! 题意简介题面链接题目大意给出两个数组a,ba,ba,b 求有多少种对应方式使得有恰好kkk对匹配(i,j)(i,j)(i,j)满足ai>bja_i>b_j ...
集合计数二项式反演_对计数数据使用负二项式
集合计数二项式反演 The Negative Binomial distribution is a discrete probability distribution that you should ...
LOJ3119 CTS2019 随机立方体概率、容斥、二项式反演
传送门为了方便我们设\(N\)是\(N,M,L\)中的最小值,某一个位置\((x,y,z)\)所控制的位置为集合\(\{(a,b,c) \mid a = x \text{或} b = y \text ...
二项式反演（学习笔记）
q w q qwq qwq机房最后一个学二项式反演的人众所周知二项式反演可以表示成 f n = ∑ i = 0 n ( − 1 ) i × C n i × g i ⟺ g n = ∑ i = 0 ...
洛谷·幼儿园篮球题【including范德蒙德卷积，二项式反演
初见安~时隔良久我又回来写多项式了[靠还是放在题目前面吧,简单讲一下这两个东西. 一.范德蒙德卷积可以理解为:在两个有n个石子和m个石子的堆里面共选k个石子的方案数.这样这个等式的成立就很显然了. ...

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

附：二元二项式反演证明

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

最新文章

热门文章