P3784-[SDOI2017]遗忘的集合【多项式ln,MTT,莫比乌斯反演】

1|0正题

题目链接:https://www.luogu.com.cn/problem/P3784

1|1题目大意

你若干个在 $[1,n]$ 的不同数字组成序列 $a$ 。

记录 $f(x)$ 表示将 $x$ 无序拆分成 $a$ 中数字的和的方案数（一个数字可以使用多次）。

现在给出所有的 $f(x)\%p\ (x\in[1,n])$ ，要求构造一组字典序最小的 $a$ 满足这个条件。

$n\leq 2^{18},p\in[10^6,2^{30})\cap Pri$

1|2解题思路

先考虑给出 $a$ 怎么求 $f$ ，考虑使用生成函数，对于一个数字 $a$ 可以表示成 $\frac{1}{1-x^{a}}$

F (x) = \sum_{i = 0}^{n} f (i) x^{i} = \prod_{i = 1}^{n} \frac{1}{1 - x^{a_{i}}}

$F(x)=\sum_{i=0}^nf(i)x^i=\prod_{i=1}^n\frac{1}{1-x^{a_i}}$

\ln F (x) = \sum_{i = 1}^{n} \ln \frac{1}{1 - x^{a_{i}}}

$\ln F(x)=\sum_{i=1}^n\ln\frac{1}{1-x^{a_i}}$

x (\ln F (x))^{'} = \sum_{i = 1}^{n} \frac{a_{i} x^{a_{i}}}{1 - x^{a_{i}}}

$x(\ \ln F(x)\ )'=\sum_{i=1}^n\frac{a_ix^{a_i}}{1-x^{a_i}}$

x (\ln F (x))^{'} = \sum_{i = 1}^{n} \sum_{j = 1}^{\infty} a_{i} x^{i j}

$x(\ \ln F(x)\ )'=\sum_{i=1}^n\sum_{j=1}^{\infty}a_ix^{ij}$

如果我们求出 $x(\ \ln F(x)\ )'$ ，就可以莫反得到每一个 $a_i$ 是否有了。

因为模数很丑，所以要用任意模数。

时间复杂度： $O(n\log n)$

1|3code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const ll N=1<<20,T=1<<15;
const double Pi=acos(-1);
struct complex{
	double x,y;
	complex(double xx=0,double yy=0)
	{x=xx;y=yy;return;}
}A[N],B[N],C[N],D[N],E[N],J[N],I[N],w[N];
struct Poly{
	ll a[N],n;
}F,G,t1;
ll n,P,cnt,pri[N/10],r[N],g[N],f[N],mu[N];
bool v[N];
complex operator+(const complex &x,const complex &y)
{return complex(x.x+y.x,x.y+y.y);}
complex operator-(const complex &x,const complex &y)
{return complex(x.x-y.x,x.y-y.y);}
complex operator*(const complex &x,const complex &y)
{return complex(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
ll power(ll x,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*x%P;
		x=x*x%P;b>>=1;
	}
	return ans;
}
void FFT(complex *f,ll n,ll op){
	for(ll i=0;i<n;i++)
		if(i<r[i])swap(f[i],f[r[i]]);
	for(ll p=2;p<=n;p<<=1){
		ll len=p>>1;
		for(ll k=0;k<n;k+=p)
			for(ll i=k;i<k+len;i++){
				complex tmp=w[n/len*(i-k)];
				if(op==-1)tmp.y=-tmp.y;
				complex tt=f[i+len]*tmp;
				f[i+len]=f[i]-tt;
				f[i]=f[i]+tt;
			}
	}
	if(op==-1){
		for(ll i=0;i<n;i++)
			f[i]=complex(f[i].x/n+0.5);
	}
	return;
}
void Mul(Poly &F,Poly &G,Poly &H){
	ll l=1;
	while(l<F.n+G.n-1)l<<=1;
	for(ll i=0;i<l;i++)
		r[i]=r[i>>1]>>1|((i&1)?(l>>1):0);
	for(ll k=1;k<l;k<<=1)
		for(ll i=0;i<k;i++)
			w[l/k*i]=complex(cos(Pi/k*i),sin(Pi/k*i));
	for(ll i=0;i<F.n;i++)A[i]=complex(F.a[i]/T,0),B[i]=complex(F.a[i]%T,0);
	for(ll i=0;i<G.n;i++)C[i]=complex(G.a[i]/T,0),D[i]=complex(G.a[i]%T,0);
	for(ll i=F.n;i<l;i++)A[i]=B[i]=complex(0,0);
	for(ll i=G.n;i<l;i++)C[i]=D[i]=complex(0,0);
	for(ll i=0;i<l;i++)E[i]=J[i]=I[i]=complex(0,0);
	FFT(A,l,1);FFT(B,l,1);FFT(C,l,1);FFT(D,l,1);
	for(ll i=0;i<l;i++){
		E[i]=A[i]*C[i];
		J[i]=A[i]*D[i]+C[i]*B[i];
		I[i]=B[i]*D[i];
	}
	FFT(E,l,-1);FFT(J,l,-1);FFT(I,l,-1);
	for(ll i=0;i<l;i++){
		H.a[i]=(ll)E[i].x*T%P*T%P;
		(H.a[i]+=(ll)J[i].x*T%P)%=P;
		(H.a[i]+=(ll)I[i].x)%=P;
	}
	H.n=F.n+G.n-1;
	return;
}
void CalcInv(ll n,Poly &F,Poly &G){
	if(n==1){G.a[0]=power(F.a[0],P-2);G.n=1;return;}
	CalcInv(n>>1,F,G);
	Mul(G,G,t1);t1.n=n;
	ll pn=F.n;F.n=n;Mul(t1,F,t1);
	t1.n=n;F.n=pn;
	for(ll i=0;i<G.n;i++)
		G.a[i]=(2ll*G.a[i]-t1.a[i]+P)%P;
	for(ll i=G.n;i<n;i++)G.a[i]=P-t1.a[i];
	G.n=n;return;
}
void GetInv(Poly &F,Poly &G){
	ll l=1;
	while(l<F.n)l<<=1;
	CalcInv(l,F,G);G.n=F.n;
	return;
}
void GetD(Poly &F){
	for(ll i=0;i<F.n-1;i++)
		F.a[i]=F.a[i+1]*(i+1)%P;
	F.n--;return;
}
void GetJ(Poly &F){
	for(ll i=F.n;i>0;i--)
		F.a[i]=F.a[i-1]*power(i,P-2)%P;
	F.a[0]=0;F.n++;return;
}
void GetLn(Poly &F){
	GetInv(F,G);
	GetD(F);
	Mul(F,G,F);GetJ(F);
	return;
}
void Prime(ll n){
	mu[1]=1;
	for(ll i=2;i<=n;i++){
		if(!v[i])pri[++cnt]=i,mu[i]=-1;
		for(ll j=1;j<=cnt&&i*pri[j]<=n;j++){
			v[i*pri[j]]=1;
			if(i%pri[j]==0)break;
			mu[i*pri[j]]=-mu[i];
		}
	}
	return;
}
signed main()
{
	scanf("%lld%lld",&n,&P);
	for(ll i=1;i<=n;i++)scanf("%lld",&F.a[i]);
	F.n=n+1;F.a[0]=1;
	GetLn(F);GetD(F);
	Prime(n);
	for(ll i=1;i<=n;i++)f[i]=F.a[i-1];
	for(ll i=1;i<=n;i++)
		for(ll j=i;j<=n;j+=i)
			(g[j]+=f[i]*mu[j/i]+P)%=P;
	ll ans=0;
	for(ll i=1;i<=n;i++)
		if(g[i])ans++;
	printf("%lld\n",ans);
	for(ll i=1;i<=n;i++)
		if(g[i])printf("%lld ",i);
	return 0;
}

__EOF__

本文作者：QuantAsk
本文链接：https://www.cnblogs.com/QuantAsk/p/16204147.html
关于博主：退役OIer，GD划水选手
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！