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

正题

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

---

题目大意

你若干个在\([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\)

---

解题思路

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

\[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}} \]

\[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_ix^{ij} \]

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

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

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

---

code

#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;
}