正题

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

---

题目大意

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

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

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

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

---

解题思路

先考虑给出$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}}$$
$$\ln F(x)=\sum _{i=1}^n\ln \frac{1}{1-x^{a_i}}$$
$$x(\ln F(x){)}^{\prime }=\sum _{i=1}^n\frac{a_i{x}^{a_i}}{1-x^{a_i}}$$
$$x(\ln F(x){)}^{\prime }=\sum _{i=1}^n\sum _{j=1}^{\infty }{a}_i{x}^{ij}$$
如果我们求出$x(\ln F(x){)}^{\prime }$，就可以莫反得到每一个$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;
}