[luogu P3784][SDOI2017]遗忘的集合

遗忘的集合

题解

生成函数万岁！！！

很明显是一道的生成函数的题，开始推吧！！！

我们知道，一个数可以被取无限次的生成函数为。

我们可以用01的序列表示是否选这个数，设每项为，则这个多项式F为

很明显，我们的目标是求出。

子命题：

证明如下：

证毕。

将原式代入并取ln

如此我们枚举，

现在就求出了，然后我们就可以用类似筛法的方式枚举倍数，大概就是有点类似欧拉筛法的方式（没了最关键的那个去重），先把每个数的值乘上它本身的编号，依次枚举它的倍数，再将它的倍数的值减去它本身的贡献，这样就能够得到它应得的值了。相信大家都明白，没什么好讲的，但T**M*E***非要让我写一下，逃...

源码

这题不能用NTT，因为模数不一定有原根，还得用FFT。

因为笔者之前打的都是NTT的模板，改FFT，害笔者调了半天。

#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<vector>
#include<map>
using namespace std;
typedef long long LL;
#define int LL
#define gc() getchar()
template<typename _T>
void read(_T &x){
	_T f=1;x=0;char s=gc();
	while(s>'9'||s<'0'){if(s=='-')f=-1;s=gc();}
	while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=gc();}
	x*=f;
}
int P;
int add(int x,int y){return x+y>=P?x+y-P:x+y;}
int dec(int x,int y){return x-y<0?x-y+P:x-y;}
int mul(int x,int y){return 1ll*x*y-1ll*x*y/P*P;}
int qkpow(int a,int s){
	int t=1;
	while(s){
		if(s&1)t=mul(t,a);
		a=mul(a,a);s>>=1;
	}
	return t;
}
const int MAXN=6e5+5;
const double PI=acos(-1.0);
int rev[MAXN],c[MAXN],d[MAXN],e[MAXN],f[MAXN],g[MAXN],h[MAXN],n,len,ans;
struct complex{
    double x,y;
    complex(double _x=0,double _y=0):x(_x),y(_y){}        
	complex operator + (complex b){return complex(x+b.x,y+b.y);}
    complex operator - (complex b){return complex(x-b.x,y-b.y);}
    complex operator * (complex b){return complex(x*b.x-y*b.y,x*b.y+y*b.x);} 
	complex operator * (double b){return complex(x*b,y*b);}   
}O[MAXN];
void FFT(int lim,complex *A,int typ){
	for(int i=0;i<lim;i++)if(i<rev[i])swap(A[i],A[rev[i]]);
	for(int i=1;i<lim;i<<=1)
		for(int j=0;j<lim;j+=(i<<1))
			for(int k=0;k<i;k++){
				complex x=A[j+k],y=O[i+k]*A[i+j+k];
				A[j+k]=x+y;A[i+j+k]=x-y;
			}
	if(typ!=-1)return ;
	reverse(A+1,A+lim);double k=1.0/lim;
	for(int i=0;i<lim;i++)A[i]=A[i]*k;
}
void Mul(int len,int *a,int *b,int *c){
	static complex A[MAXN],B[MAXN],C[MAXN],D[MAXN],F[MAXN],G[MAXN],H[MAXN];
	int lim=1,L=0;while(lim<(len<<1))lim<<=1,L++;
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<L-1);
	for(int i=1;i<lim;i<<=1)
		for(int k=0;k<i;k++)
			O[i+k]=complex(cos(PI*k/i),sin(PI*k/i));
	for(int i=0;i<len;i++){
		A[i].x=a[i]>>15;B[i].x=a[i]&32767;
		C[i].x=b[i]>>15;D[i].x=b[i]&32767;
		A[i].y=B[i].y=C[i].y=D[i].y=0;
	}
	for(int i=len;i<lim;i++)A[i]=B[i]=C[i]=D[i]=0;
	FFT(lim,A,1);FFT(lim,B,1);FFT(lim,C,1);FFT(lim,D,1);
	for(int i=0;i<lim;i++)F[i]=A[i]*C[i],G[i]=A[i]*D[i]+B[i]*C[i],H[i]=B[i]*D[i];
	FFT(lim,F,-1);FFT(lim,G,-1);FFT(lim,H,-1);
	for(int i=0;i<len;i++)c[i]=(((LL)(F[i].x+0.5)%P<<30)+((LL)(G[i].x+0.5)<<15)+((LL)(H[i].x+0.5)))%P;
	for(int i=len;i<lim;i++)c[i]=0;
}
void polyinv(int len,int *a,int *b){
	if(len==1)return (void)(b[0]=qkpow(a[0],P-2));
	polyinv(len>>1,a,b);Mul(len,a,b,c);Mul(len,c,b,d);
	for(int i=0;i<len;i++)b[i]=dec(add(b[i],b[i]),d[i]);
}
void polyln(int len,int *a,int *b){
	for(int i=1;i<len;i++)e[i-1]=mul(a[i],i);
	polyinv(len,a,f);Mul(len,e,f,b);
	for(int i=len-1;i>0;i--)b[i]=mul(b[i-1],qkpow(i,P-2));
}
signed main(){
	read(n);read(P);
	int len=1;while(len<=n)len<<=1;
	for(int i=1;i<=n;i++)read(g[i]);
	g[0]=1;polyln(len,g,h);
	for(int i=1;i<=n;i++)h[i]=mul(i,h[i]);
	for(int i=1;i<=n;i++)
		for(int j=i<<1;j<=n;j+=i)
			h[j]=dec(h[j],h[i]);
	for(int i=1;i<=n;i++)if(h[i])++ans;
	printf("%d\n",ans);
	for(int i=1;i<=n;i++)if(h[i])printf("%d ",h[i]);
    return 0;
}

谢谢！！！