【BZOJ 3684】大朋友和多叉树

题目

拉格朗日反演的板子题;

拉格朗日反演用来求多项式的复合逆,大概就是有两个多项式(G(x),F(x)),满足(g_0=0,f_0=0),且(F(G(x))=x),或是(G(F(x))=x);就有

[[x^n]F(x)=frac{1}{n}[x^{-1}]frac{1}{G^{n}(x)} ]

证明看这位老鸽的吧。

有点吓人的是(x^{-1})是个啥,还有(G^n(x))没有逆咋求啊;

先让右边整体平移一个(x^n),有

[[x^n]F(x)=frac{1}{n}[x^{n-1}]frac{x^n}{G^{n}(x)} ]

我们发现(x^n)(G^n(x))都没有常数项,于是我们可以都往左平移一位,就有

[[x^n]F(x)=frac{1}{n}[x^{n-1}]frac{1}{(frac{G(x)}{x})^n} ]

于是我们利用多项式(exp)就能在(O(nlog n))的时间内求出(F(x))某一项的值了;

对于这道题,我们设答案的生成函数为(F(x))(G(x)=sum_{iin S}x^i),不难发现几个儿子就是搞几个(F(x))卷起来;

于是(G(F(x))+x=F(x)),移项可得(F(x)-G(F(x))=x),构造一个多项式(D(x)=x-sum_{iin S}x^i),则(D(F(x))=x),套上面的反演公式就可以了;

代码

#include<bits/stdc++.h>
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read() {
	char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
	while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();return x;
}
const int maxn=262144+5;
const int mod=950009857;const int G[]={7,(mod+1)/7};
int rev[maxn],C[maxn],D[maxn],H[maxn],K[maxn],d[maxn];
int a[maxn],b[maxn],_o[2][105],inv[maxn],n,len,m,c[maxn];
inline int dqm(int x) {return x<0?x+mod:x;}
inline int qm(int x) {return x>=mod?x-mod:x;}
inline int ksm(int a,int b) {
	int S=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)S=1ll*S*a%mod;return S;
}
inline void NTT(int *f,int o) {
	for(re int i=0;i<len;i++)if(i<rev[i])std::swap(f[i],f[rev[i]]);
	for(re int w=0,i=2;i<=len;i<<=1,++w) {
		int ln=i>>1,og1=_o[o][w];
		if(!og1) og1=_o[o][w]=ksm(G[o],(mod-1)/i);
		for(re int og=1,t,l=0;l<len;l+=i,og=1) 
			for(re int x=l;x<l+ln;++x) {
				t=1ll*f[x+ln]*og%mod;og=1ll*og*og1%mod;
				f[x+ln]=dqm(f[x]-t);f[x]=qm(f[x]+t);
			}
	}if(!o)return;int Inv=ksm(len,mod-2);
	for(re int i=0;i<len;i++)f[i]=1ll*f[i]*Inv%mod;
}
inline void init(int n) {
	len=1;while(len<n+n-1)len<<=1;
	for(re int i=0;i<len;++i)rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
}
void Inv(int n,int *A,int *B) {
	if(n==1) {B[0]=ksm(A[0],mod-2);return;}
	Inv((n+1)>>1,A,B);init(n);
	for(re int i=0;i<n;i++)C[i]=A[i];NTT(B,0),NTT(C,0);
	for(re int i=0;i<len;i++)B[i]=dqm(2ll*B[i]%mod-1ll*C[i]*B[i]%mod*B[i]%mod);
	NTT(B,1);for(re int i=0;i<len;i++)C[i]=0;
	for(re int i=n;i<len;i++)B[i]=0;
}
inline void Ln(int n,int *A,int *B) {
	for(re int i=1;i<n;i++)D[i-1]=1ll*i*A[i]%mod;
	Inv(n,A,H);len=1;init(n);NTT(H,0),NTT(D,0);
	for(re int i=0;i<len;i++)B[i]=1ll*D[i]*H[i]%mod,D[i]=H[i]=0;
	NTT(B,1);for(re int i=n-1;i;--i)B[i]=1ll*B[i-1]*inv[i]%mod;B[0]=0;
	for(re int i=n;i<len;i++)B[i]=0;
}
void Exp(int n,int *A,int *B) {
	if(n==1) {B[0]=1;return;}
	Exp((n+1)>>1,A,B);Ln(n,B,K);init(n);
	for(re int i=0;i<n;i++)K[i]=dqm(A[i]-K[i]);
	K[0]=qm(K[0]+1);NTT(K,0),NTT(B,0);
	for(re int i=0;i<len;++i)B[i]=1ll*B[i]*K[i]%mod,K[i]=0;
	NTT(B,1);for(re int i=n;i<len;i++)B[i]=0;
}
int main() {
	scanf("%d%d",&n,&m);inv[1]=1;
	for(re int i=2;i<=n;++i)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
	for(re int x,i=1;i<=m;i++) {
		scanf("%d",&x);if(x-1<n) a[x-1]=mod-1;
	}
	a[0]=1;Ln(n,a,c);for(re int i=0;i<n;i++)c[i]=1ll*c[i]*n%mod;Exp(n,c,b);
	Inv(n,b,d);printf("%d
",1ll*inv[n]*d[n-1]%mod);
	return 0;
}
原文地址:https://www.cnblogs.com/asuldb/p/12307387.html