欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

20200802 T3 我永远喜欢【生成函数容斥,拉格朗日反演】

程序员文章站 2022-03-23 15:33:34
题目描述有 nnn 种颜色的石子,每种 cic_ici​ 个,记一个石子序列首尾相接后极长连续段的长度为 lil_ili​,求所有石子序列的 1∏li!\frac 1{\prod l_i!}∏li​!1​ 的和。n≤105,∑ci≤2∗105n\le10^5,\sum c_i\le2*10^5n≤105,∑ci​≤2∗105题目分析先考虑去掉首尾相接的情况怎么做。因为限制了极长,容易想到把每种颜色分成几段,然后合并,但是不好保证相同颜色不被并在一起。先不谈容斥做法,题解给出了一种用生成函数解决的...

题目描述

nn 种颜色的石子,每种 cic_i 个,记一个石子序列首尾相接后极长连续段的长度为 lil_i,求所有石子序列的 1li!\frac 1{\prod l_i!} 的和。

n105,ci2105n\le10^5,\sum c_i\le2*10^5

题目分析

先考虑去掉首尾相接的情况怎么做。

因为限制了极长,容易想到把每种颜色分成几段,然后合并,但是不好保证相同颜色不被并在一起。

先不谈容斥做法,题解给出了一种用生成函数解决的方法:
一个极长连续段的权值是 1len!\frac 1{len!}
我们准备将每种颜色划分成了许多段然后并起来。考虑给每种长度的段给一个权值 aia_i,使得对于一个长度为 lenlen 的极长颜色段的所有划分中ai\prod a_i 的和等于 1len!\frac 1{len!},即:

F=aixi,G=ex1F=\sum a_ix^i,G=e^x-1
那么有 F+F2+F3...=F1F=GF+F^2+F^3...=\frac F{1-F}=G

从而 F=GG+1=ex1exF=\frac G{G+1}=\frac {e^x-1}{e^x}

于是问题变成了:把每种颜色划分成若干段,然后把不同颜色的所有段一起排列,排列后相邻段颜色可以相同,长度为 ii 的段的权值是 aia_i,整个序列的权值是每段的乘积,求所有方案的序列权值之和。

那么我们只需要对每种颜色求出划分为 ii 段的权值之和,然后不同的段用 EGF 乘起来即可。
mm 个石子划分为 ii 段的权值之和即 gi=[xm]Fig_i=[x^m]F^i


然后考虑加上首尾相接的条件,在环上把 1 号点标出来,相当于现在会有颜色段跨越 NN11,如果按序列的方法做会漏掉这些方案的贡献,考虑循环移位,把跨越 NN11 的颜色段顺次往后移,直到它和上一段的分界点为 N ][ 1N~][~1,那么一个这样的方案唯一对应了一个序列上的方案。

所以只需要把第一段的权值额外乘上一个长度,就可以算到跨环的贡献。

那么对于一个有 mm 个石子的颜色,当它的第一段为整个序列的第一段时,它分 ii 段的贡献之和应当是:
gi=[xm]Fi1xF=[xm1](Fi)1i=[xm]Fimig_i=[x^m]F^{i-1}xF'\\=[x^{m-1}](F^i)'*\frac 1i\\ =[x^m]F^i*\frac mi

求出 gig_i 后得到对应的两个 EGF,然后分治NTT合并即可。

i[1,n],[xn]Fii\in[1,n],[x^n]F^i 求出 FF 的复合逆就可以用扩展拉格朗日反演计算,具体做法写在了这篇博客

FF 的复合逆 HH 满足 eH1eH=x{e^H-1\over e^H}=x,即 eH=11xe^H=\frac 1{1-x},所以 H=ln11x=i11ixiH=\ln {\frac 1{1-x}}=\sum_{i\ge 1}\frac 1ix^i

Code:

#include<bits/stdc++.h>
#define maxn 550005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
using namespace std;
void output(vector<int>&a){
	rep(i,0,a.size()-1) printf("%d ",a[i]); puts("");
}
const int mod = 998244353;
int w[maxn],r[maxn],WL,lg[maxn],fac[maxn],inv[maxn],invf[maxn];
int Pow(int a,int b){int s=1;for(;b;b>>=1,a=1ll*a*a%mod) b&1&&(s=1ll*s*a%mod);return s;}
void init(int n){
	WL=1; while(WL<=n) WL<<=1;
	w[WL>>1]=1; int wn=Pow(3,(mod-1)/WL);
	for(int i=WL/2+1;i<WL;i++) w[i]=1ll*w[i-1]*wn%mod;
	for(int i=WL/2-1;i>=1;i--) w[i]=w[i<<1];
	fac[0]=fac[1]=inv[0]=inv[1]=invf[0]=invf[1]=1;
	rep(i,2,WL-1) fac[i]=1ll*fac[i-1]*i%mod,inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod,invf[i]=1ll*invf[i-1]*inv[i]%mod,lg[i]=lg[i>>1]+1;
}
int upd(int x){return x+(x>>31&mod);}
void NTT(int *a,int len,int flg){
	if(flg^1) reverse(a+1,a+len);
	for(int i=0;i<len;i++) if(i<(r[i]=r[i>>1]>>1|(i&1?len>>1:0))) swap(a[i],a[r[i]]);
	for(int i=2,l=1,v;i<=len;l=i,i<<=1)
		for(int j=0;j<len;j+=i)
			for(int k=j,o=l;k<j+l;k++,o++)
				a[k+l]=upd(a[k]-(v=1ll*w[o]*a[k+l]%mod)),a[k]=upd(a[k]+v-mod);
	if(flg^1) for(int i=0,Inv=mod-(mod-1)/len;i<len;i++) a[i]=1ll*a[i]*Inv%mod;
}
typedef vector<int> Poly;
Poly operator * (const Poly &A,const Poly &B){
	static int a[maxn],b[maxn]; int n=A.size(),m=B.size(),L=1<<(lg[n+m-2]+1);
	if(n<=20||m<=20||n+m<=200) {rep(i,0,n+m-2) a[i]=0; rep(i,0,n-1) rep(j,0,m-1) a[i+j]=(a[i+j]+1ll*A[i]*B[j])%mod;}
	else{
		rep(i,0,L-1) a[i]=i<n?A[i]:0,b[i]=i<m?B[i]:0;
		NTT(a,L,1),NTT(b,L,1);
		rep(i,0,L-1) a[i]=1ll*a[i]*b[i]%mod;
		NTT(a,L,-1);
	}
	return Poly(a,a+n+m-1);
}
Poly operator + (Poly A,const Poly &B){
	if(B.size()>A.size()) A.resize(B.size());
	rep(i,0,A.size()-1) A[i]=upd(A[i]+B[i]-mod);
	return A;
}
Poly INV(const Poly &A,int n){
	static int B[maxn],a[maxn]; B[B[1]=0]=Pow(A[0],mod-2);
	for(int k=2,L=4;k<n<<1;k=L,L<<=1){
		memcpy(a,&A[0],min(n,k)<<2); rep(i,k,L-1) a[i]=B[i]=0;
		NTT(a,L,1),NTT(B,L,1); rep(i,0,L-1) B[i]=1ll*B[i]*upd(2-1ll*a[i]*B[i]%mod)%mod; NTT(B,L,-1);
		memset(B+k,0,k<<2);
	}
	return Poly(B,B+n);
}
Poly LN(const Poly &A,int n){
	Poly a(n-1),b(INV(A,n));
	rep(i,0,n-2) a[i]=1ll*A[i+1]*(i+1)%mod;
	a=a*b,b[0]=0;
	rep(i,1,n-1) b[i]=1ll*a[i-1]*inv[i]%mod;
	return b;
}
Poly EXP(const Poly &A,int n){
	Poly B{1},b;
	for(int k=2,L=4;k<n<<1;k=L,L<<=1){
		B.resize(k),b=LN(B,k); rep(i,0,min(n,k)-1) b[i]=upd((i==0)-b[i]+A[i]);
		B=B*b,B.resize(min(n,k));
	}
	return B;
}
Poly L;
void calc(Poly &g1,Poly &g2,int n){
	g1.resize(n+1),g2.resize(n+1);
	Poly A(L.begin(),L.begin()+n);
	rep(i,0,n-1) A[i]=1ll*A[i]*(mod-n)%mod;
	A=EXP(A,n);
	rep(i,1,n) g2[i]=1ll*A[n-i]*invf[i-1]%mod,g1[i]=1ll*A[n-i]*invf[i-1]%mod*inv[n]%mod;
}
Poly g1[maxn],g2[maxn];
void solve(int l,int r){
	if(l==r) return;
	int mid=(l+r)>>1;
	solve(l,mid),solve(mid+1,r);
	Poly A(g1[l]*g1[mid+1]),B(g1[l]*g2[mid+1]+g2[l]*g1[mid+1]);
	g1[l].swap(A),g2[l].swap(B);
}
int n,c[maxn],Mx,Sum;
int main()
{
	freopen("always.in","r",stdin);
	freopen("always.out","w",stdout);
	scanf("%d",&n);
	for(int i=1;i<=n;i++) scanf("%d",&c[i]),Mx=max(Mx,c[i]),Sum+=c[i];
	init(max(2*Mx,Sum));
	L=LN(Poly(inv+1,inv+Mx+1),Mx);
	for(int i=1;i<=n;i++) calc(g1[i],g2[i],c[i]);
	solve(1,n);
	int ans=0;
	rep(i,n,g2[1].size()-1) ans=(ans+1ll*g2[1][i]*fac[i-1])%mod;
	printf("%d\n",ans);
}

这里是容斥的做法,大致思路是求容斥系数的时候通过每种段数被算了多少次来计算,不过在环上的时候还要保证最后一段和第二段与第一段颜色不同,而上面这种做法则不需要考虑这个问题。

本文地址:https://blog.csdn.net/C20181220_xiang_m_y/article/details/107898725