[BZOJ4913][SDOI2017]遗忘的集合

luogu
bzoj

sol

我们设\(a_i\in\{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_i\ln(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_i\sum_{j=1}^{\infty}\frac{x^{ij}}{j} \]

令\(T=ij\)，代入并交换求和顺序得

\[\ln F(x)=\sum_{T=1}^{\infty}x^T\sum_{i|T}a_i\times \frac iT \]

对\(F(x)\)求对数后就可以知道\(\sum_{i|T}a_i\times\frac iT\)，再反演一下（可以不用莫比乌斯反演，直接\(O(n\ln 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\n",Ans);
	for (int i=1;i<=n;++i) if (f[i]) printf("%d ",i);
	return puts(""),0;
}