[SDOI2017]遗忘的集合(多项式ln+生成函数+莫比乌斯反演)

[SDOI2017]遗忘的集合(多项式ln+生成函数+莫比乌斯反演)

题面

略

分析

设\(a_i=[i \in S]\),那么元素\(i\)的生成函数为\((\frac{1}{1-x^i})^{a_i}\),答案的生成函数为\(f(x)=\prod_{i \geq 1}(\frac{1}{1-x^i})^{a_i}\). 现在题目已经给出了\(f(x)\)的各项系数,求\(a_i\)

为了把乘法化成加法,两边取对数,得到:

\[-\ln F(x)=\sum_{i \geq 1}a_i\ln(1-x^i) \]

根据\(\ln\)的泰勒展开\(\ln(1-x)=-\sum_{j \geq 1}\frac{x^j}{j}\)

\[\ln F(x)=\sum_{i\geq 1}a_i \sum_{j \geq 1}\frac{x^{ij}}{j} \]

交换求和顺序,令\(ij=k\),则\(j=\frac{k}{i}\)

\[\ln F(x)=\sum_{k \geq 1}x^k \sum_{i|k} a_i \frac{i}{k} \]

设\(\ln F(x)\)的第\(i\)项系数为\(g_i\),则\(ng_n=\sum_{i|n} ia_i\)

看到约数求和,想到莫比乌斯反演:\(na_n=\sum_{i|n} ig_i \mu(\frac{n}{i})\). 那么我们就可以求出\(na_n\),又因为\(a_n\)只能为0或1,当反演结果不为0的时候输出即可。显然这个解字典序最小。

直接套多项式板子,因为模数不是NTT模数,可以拆系数FFT(俗称MTT),复杂度\(O(n\log n)\)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#define maxn (1<<19)
using namespace std;
const double pi=acos(-1.0);
typedef long long ll;
ll mod; 
inline ll fast_pow(ll x,ll k){
	ll ans=1;
	while(k){
		if(k&1) ans=ans*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return ans;
}
inline ll inv(ll x){
	return fast_pow(x,mod-2);
}
struct com{
	double real;
	double imag;
	com(){
		
	} 
	com(double _real,double _imag){
		real=_real;
		imag=_imag;
	}
	friend com operator + (com p,com q){
		return com(p.real+q.real,p.imag+q.imag);
	}
	friend com operator - (com p,com q){
		return com(p.real-q.real,p.imag-q.imag);
	}
	friend com operator * (com p,com q){
		return com(p.real*q.real-p.imag*q.imag,p.real*q.imag+p.imag*q.real);
	} 
	friend com operator * (com p,double k){
		return com(p.real*k,p.imag*k);
	}
	friend com operator / (com p,double k){
		return com(p.real/k,p.imag/k);
	}
	inline com conj(){
		return com(real,-imag);
	}
};


int rev[maxn+5];
com w[maxn+5];
void fft(com *x,int n){
	for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
	for(int len=1;len<n;len*=2){
		int sz=len*2;
		for(int l=0;l<n;l+=sz){
			int r=l+len-1;
			for(int i=l;i<=r;i++){
				com tmp=x[i+len];
				x[i+len]=x[i]-tmp*w[n/sz*(i-l)];
				x[i]=x[i]+tmp*w[n/sz*(i-l)];
			}
		}
	}
}

void mul(ll *a,ll *b,ll *c,int n,int m){
	static com p[maxn+5],q[maxn+5],r[maxn+5],s[maxn+5];
	int N=1,L=0;
	while(N<n+m-1){
		N*=2;
		L++;
	}
	for(int i=0;i<n;i++){
		ll ta=(a[i]+mod)%mod;
		p[i]=com(ta>>15,ta&((1<<15)-1));
	}
	for(int i=n;i<N;i++) p[i]=com(0,0);
	for(int i=0;i<m;i++){
		ll tb=(b[i]+mod)%mod;
		q[i]=com(tb>>15,tb&((1<<15)-1));
	}
	for(int i=m;i<N;i++) q[i]=com(0,0);
	for(int i=0;i<N;i++) w[i]=com(cos(2*pi*i/N),sin(2*pi*i/N));
	for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	fft(p,N);
	fft(q,N);
	for(int i=0;i<N;i++){
		int j=(i==0?0:N-i);
		com da=(p[i]+p[j].conj())*com(0.5,0);
		com db=(p[i]-p[j].conj())*com(0,-0.5);
		com dc=(q[i]+q[j].conj())*com(0.5,0);
		com dd=(q[i]-q[j].conj())*com(0,-0.5);
		r[j]=da*dc+da*dd*com(0,1);
		s[j]=db*dc+db*dd*com(0,1);
	}
	fft(r,N);
	fft(s,N);
	for(int i=0;i<n+m-1;i++){
		ll ac=(ll)(r[i].real/N+0.5)%mod;
		ll ad=(ll)(r[i].imag/N+0.5)%mod;
		ll bc=(ll)(s[i].real/N+0.5)%mod;
		ll bd=(ll)(s[i].imag/N+0.5)%mod;
		c[i]=(((ac<<30)+((ad+bc)<<15)+bd)%mod+mod)%mod;
	}
}
void poly_inv(ll *f,ll *g,int n){
	static ll tmp[maxn+5];
	if(n==1){
		g[0]=inv(f[0]);
		return;
	}
	poly_inv(f,g,(n+1)/2);
	mul(f,g,tmp,n,n);
	mul(tmp,g,tmp,n,n);
	for(int i=0;i<n;i++) g[i]=(2*g[i]-tmp[i]+mod)%mod;//tmp[i]=f[i]*g[i]^2
} 
void poly_deriv(ll *f,ll *g,int n){
	for(int i=1;i<n;i++) g[i-1]=f[i]*i%mod;
	g[n-1]=0;
} 
void poly_inter(ll *f,ll *g,int n){
	for(int i=n-1;i>=1;i--) g[i]=f[i-1]*inv(i)%mod;
	g[0]=0;
}
void poly_ln(ll *f,ll *g,int n){
	static ll inv_ln[maxn+5];
	poly_deriv(f,g,n);
	poly_inv(f,inv_ln,n);
	mul(g,inv_ln,g,n,n);
	poly_inter(g,g,n*2);
}

int cnt;
int mu[maxn+5],prime[maxn+5];
bool vis[maxn+5];
void sieve_mu(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prime[++cnt]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}else mu[i*prime[j]]=-mu[i];
		}
	}
} 

int n;
ll f[maxn+5],lnf[maxn+5];
ll a[maxn+5];//实际上存的是a[i]*i
vector<int>ans;
int main(){
	scanf("%d",&n);
	scanf("%lld",&mod);
	sieve_mu(n);
	f[0]=1;
	for(int i=1;i<=n;i++) scanf("%lld",&f[i]);
	poly_ln(f,lnf,n+1);//对f的生成函数求ln 
	for(int i=0;i<=n;i++) lnf[i]=lnf[i]*i%mod;
	for(int i=1;i<=n;i++){
		for(int j=1;j*i<=n;j++) a[i*j]+=lnf[i]*mu[j];//莫比乌斯反演 
	}
	for(int i=1;i<=n;i++) if(a[i]) ans.push_back(i);
	printf("%d\n",(int)ans.size());
	for(int i=0;i<(int)ans.size();i++) printf("%d ",ans[i]);
}