sol
我们设(a_iin{0,1})表示(i)这个数有没有出现在集合中。那么(f)对应的生成函数就是：
[F(x)=prod_{i=1}^{n}(frac{1}{1-x^i})^{a_i}
]
现在给出(F(x))，要求(a_i)。
先两边取对数。
[-ln F(x)=sum_{i=1}^na_iln(1-x^i)
]
然后有一个这样子的式子：
[ln(1-x^i)=-sum_{j=1}^{infty}frac{x^{ij}}{j}
]
证明：（感谢@Cyhlnj ）
[ln F(x)=G(x)\frac{F'(x)}{F(x)}=G'(x)\frac{-ix^{i-1}}{1-x^i}=G'(x)\-sum_{j=0}^{infty} ix^{i-1+ij}=G'(x)\-sum_{j=0}^{infty}frac{ix^{i+ij}}{i+ij}=G(x)\-sum_{j=1}^{infty}frac{x^{ij}}{j}=G(x)
]
接上，把这个式子代入得
[-ln F(x)=-sum_{i=1}^na_isum_{j=1}^{infty}frac{x^{ij}}{j}
]
令(T=ij)，代入并交换求和顺序得
[ln F(x)=sum_{T=1}^{infty}x^Tsum_{i|T}a_i imes frac iT
]
对(F(x))求对数后就可以知道(sum_{i|T}a_i imesfrac iT)，再反演一下（可以不用莫比乌斯反演，直接(O(nln n))就可以了）就可以求出(a_i)了。
code
(BZOJ)上跑不过。
果然是人丑常数大。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
int gi(){
int x=0,w=1;char ch=getchar();
while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
if (ch=='-') w=0,ch=getchar();
while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
return w?x:-x;
}
const int N = 6e5+5;
const double Pi = acos(-1);
struct Complex{
double rl,im;
Complex(){rl=im=0;}
Complex(double x,double y){rl=x,im=y;}
Complex conj(){return Complex(rl,-im);}
Complex operator + (Complex b)
{return Complex(rl+b.rl,im+b.im);}
Complex operator - (Complex b)
{return Complex(rl-b.rl,im-b.im);}
Complex operator * (Complex b)
{return Complex(rl*b.rl-im*b.im,im*b.rl+rl*b.im);}
}w[N],s1[N],s2[N],s3[N],s4[N],s5[N],s6[N];
int n,mod,qmod,rev[N],f[N],Ans;
void fft(Complex *P,int len){
for (int i=0;i<len;++i) if (i<rev[i]) swap(P[i],P[rev[i]]);
for (int i=1;i<len;i<<=1)
for (int p=i<<1,j=0;j<len;j+=p)
for (int k=0;k<i;++k){
Complex x=P[j+k],y=w[len/i*k]*P[j+k+i];
P[j+k]=x+y;P[j+k+i]=x-y;
}
}
void mul(int *a,int *b,int *c,int len){
int len2=len;len<<=1;
int l=0;while ((1<<l)<len) ++l;--l;
for (int i=0;i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<l);
for (int i=0;i<len;++i) w[i]=Complex(cos(Pi*i/(len)),sin(Pi*i/(len)));
for (int i=0;i<len2;++i) s1[i]=(Complex){a[i]&32767,a[i]>>15};
for (int i=0;i<len2;++i) s2[i]=(Complex){b[i]&32767,b[i]>>15};
for (int i=len2;i<len;++i) s1[i]=s2[i]=Complex(0,0);
fft(s1,len);fft(s2,len);
for (int i=0;i<len;++i){
int j=(len-i)&(len-1);
Complex da=(s1[i]+s1[j].conj())*Complex(0.5,0);
Complex db=(s1[i]-s1[j].conj())*Complex(0,-0.5);
Complex dc=(s2[i]+s2[j].conj())*Complex(0.5,0);
Complex dd=(s2[i]-s2[j].conj())*Complex(0,-0.5);
s3[j]=da*dc;s4[j]=da*dd;s5[j]=db*dc;s6[j]=db*dd;
}
for (int i=0;i<len;++i) s1[i]=s3[i]+s4[i]*Complex(0,1);
for (int i=0;i<len;++i) s2[i]=s5[i]+s6[i]*Complex(0,1);
fft(s1,len);fft(s2,len);
for (int i=0;i<len2;++i){
int da=(long long)(s1[i].rl/len+0.5)%mod;
int db=(long long)(s1[i].im/len+0.5)%mod;
int dc=(long long)(s2[i].rl/len+0.5)%mod;
int dd=(long long)(s2[i].im/len+0.5)%mod;
c[i]=(da+((long long)(db+dc)<<15)+((long long)dd<<30))%mod;
}
for (int i=len2;i<len;++i) c[i]=0;
}
int fastpow(int a,int b){
int res=1;
while (b) {if (b&1) res=1ll*res*a%mod;a=1ll*a*a%mod;b>>=1;}
return res;
}
int tmp[N];
void GetInv(int *a,int *b,int len){
if (len==1) {b[0]=fastpow(a[0],mod-2);return;}
GetInv(a,b,len>>1);mul(a,b,tmp,len);
for (int i=0;i<len;++i) tmp[i]=mod-tmp[i];tmp[0]+=2;
mul(tmp,b,b,len);
}
int d[N],inv[N];
void Getln(int *a,int *b,int len){
for (int i=1;i<len;++i) d[i-1]=1ll*i*a[i]%mod;
GetInv(a,inv,len);mul(d,inv,b,len);
for (int i=len-1;i;--i) b[i]=1ll*b[i-1]*fastpow(i,mod-2)%mod;
}
int main(){
n=gi();mod=gi();qmod=sqrt(mod);
int len=1;while (len<=n) len<<=1;
for (int i=1;i<=n;++i) f[i]=gi();
f[0]=1;Getln(f,f,len);
for (int i=1;i<=n;++i) f[i]=1ll*f[i]*i%mod;
for (int i=1;i<=n;++i)
for (int j=i+i;j<=n;j+=i)
f[j]=(f[j]-f[i]+mod)%mod;
for (int i=1;i<=n;++i) if (f[i]) ++Ans;
printf("%d
",Ans);
for (int i=1;i<=n;++i) if (f[i]) printf("%d ",i);
return puts(""),0;
}