[SDOI2017]遗忘的集合

题面

DP学得好好的为什么突然想起写多项式呢？

这要从几天前说起

一开始是看到一道FWT的DP([SNOI2017]遗失的答案)

然后在洛谷上搜

没有搜到

但是搜到了这个题

仔细一看

MTT+莫比乌斯反演!

好像很有意思的样子

顺便写一写(4)次DFT的MTT

嘿嘿嘿……

于是LYC沉迷多项式无法自拔

说这个题

这个题是真的毒瘤

明明就有唯一解

偏要说什么输出字典序最小的

现在考虑计算方案数

一般来说会想到DP

然而总不可能DP回去吧

这时LYC耳边响起了语文老师的谆谆教诲

生命不能让阅(数)读(学)缺席

阅(数)读(学)让我们找到回家的路

所以只要找到方案数与集合的函数关系

反演将“让我们找到回家的路”

---

现在我们的问题是

已知集合求方案数的表达式

似乎除了生成函数也没什么好用的

设(a_i)表示(i)是否在集合中

设序列(a)的生成函数为(A)

方案数序列的生成函数为(F)

[F(x)=prod_{i=1}^{n}{(frac{1}{1-x^i})^{a_i}} \ ln(F(x))=sum_{i=1}^{n}{(a_i*ln(frac{1}{1-x^i}))} \-ln(F(x))=sum_{i=1}^{n}{(a_i*ln({1-x^i}))} ]

考虑计算函数(G_n(x)=ln(1-x^n))

根据链式法则

[G_n'(x)=frac{-nx^{n-1}}{1-x^n} \ecause frac{1}{1-x^n}=sum_{i=1}^{infty}{x^{i*n}} \ herefore G_n'(x)=sum_{i=0}^{infty}{-nx^{i*n+(n-1)}} \ herefore G_n(x)=sum_{i=0}^{infty}frac{-nx^{i*n+(n-1)+1}}{i*n+(n-1)+1} \=-sum_{i=0}^{infty}{frac{nx^{i*n+n}}{i*n+n}} \=-sum_{i=0}^{infty}{frac{x^{(i+1)*n}}{i+1}} \=-sum_{i=1}^{infty}{frac{x^{i*n}}{i}} ]

带入原式

[\-ln(F(x))=-sum_{i=1}^{n}{(a_i*sum_{j=0}^{infty}{frac{x^{ij}}{j}})} \ln(F(x))=sum_{t=1}^{infty}{x^t*sum_{i|t}{a_i*frac{i}{t}}} ]

令(ln(F(x)))对应的数列为(f)

[n*f_n=sum_{d|n}{a_d*d} ]

莫比乌斯反演即可求出(a_i)

注意只保证(p)为质数

要写MTT

---

话说洛谷上有个讨论问什么输出(a_i)是对的(本来要输出(i))

因为程序里反演求出的序列是({a_i*i})

而本来的(a_i)是零一串

当然求出的(a_i)要么是(0)要么是(i)啦

#include<bits/stdc++.h>

using namespace std;

#define gc c=getchar()
#define ll long long
#define db double 

template<typename T>
inline void read(T&x){
    x=0;T k=1;char gc;
    while(!isdigit(c)){if(c=='-')k=-1;gc;}
    while(isdigit(c)){x=x*10+c-'0';gc;}x*=k;
}

struct comp{
    db r,i;
    comp(){}
    comp(const db &_r,const db &_i):r(_r),i(_i){}
};

inline comp operator + (const comp &a,const comp &b){
    return comp(a.r+b.r,a.i+b.i);
}

inline comp operator - (const comp &a,const comp &b){
    return comp(a.r-b.r,a.i-b.i);
}

inline comp operator * (const comp &a,const comp &b){
    return comp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
}

inline comp operator / (const comp &a,const db &b){
    return comp(a.r/b,a.i/b);
}

inline comp conj(const comp &a){
    return comp(a.r,-a.i);
}

const db PI=acos(-1.0);
const int N=1<<21|7;

int rev[N];
comp Wn[N];
int Now_Len;
inline void revinit(int n){
    for(int i=0;i<n;++i)rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
    for(int i=0;i<n;++i)Wn[i]=comp(cos(PI*i/n),sin(PI*i/n));
    Now_Len=n;
}

inline void DFT(comp *A,int n){
    if(Now_Len!=n)revinit(n);
    for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int i=1,l=0;i<n;i<<=1,l++){
        for(int j=0;j<n;j+=i<<1){
            for(int k=0;k<i;++k){
                comp u=A[j+k],v=A[i+j+k]*Wn[(n>>l)*k];
                A[j+k]=u+v;
                A[i+j+k]=u-v;
            }
        }
    }
}

inline void IDFT(comp *A,int n){
    if(Now_Len!=n)revinit(n);
    for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
    for(int i=1,l=0;i<n;i<<=1,l++){
        for(int j=0;j<n;j+=i<<1){
            for(int k=0;k<i;++k){
                comp u=A[j+k],v=A[i+j+k]*conj(Wn[(n>>l)*k]);
                A[j+k]=u+v;
                A[i+j+k]=u-v;
            }
        }
    }
    for(int i=0;i<n;++i){
        A[i]=A[i]/n;
    }
}

#define DFT(A) DFT(A,n)
#define IDFT(A) IDFT(A,n)

int p;
comp a[N],b[N],c[N],d[N];
inline void MTT(int *A,int *B,int *C,int lenA,int lenB){
    int lenC=lenA+lenB-1,n;
    for(n=1;n<lenC;n<<=1);
    for(int i=0;i<lenA;++i)a[i]=comp(A[i]&32767,A[i]>>15);
    for(int i=0;i<lenB;++i)b[i]=comp(B[i]&32767,B[i]>>15);
    for(int i=lenA;i<n;++i)a[i].r=a[i].i=0;
    for(int i=lenB;i<n;++i)b[i].r=b[i].i=0;
    DFT(a);DFT(b);
    for(int i=0;i<n;++i){
        int j=(n-i)&(n-1);
        comp A0=(a[i]+conj(a[j]))*comp(0.5,0);
        comp A1=(a[i]-conj(a[j]))*comp(0,-0.5);
        comp B0=(b[i]+conj(b[j]))*comp(0.5,0);
        comp B1=(b[i]-conj(b[j]))*comp(0,-0.5);
        c[i]=A0*B0+A0*B1*comp(0,1);
        d[i]=A1*B0+A1*B1*comp(0,1);
    }
    IDFT(c);IDFT(d);
    for(int i=0;i<lenC;++i){
        ll s1=(ll)(c[i].r+0.5)%p;
        ll s2=(ll)(c[i].i+0.5)%p;
        ll s3=(ll)(d[i].r+0.5)%p;
        ll s4=(ll)(d[i].i+0.5)%p;
        C[i]=(s1+((s2+s3)<<15)+(s4<<30))%p;
    }
}

inline ll qpow(ll a,ll b){
    ll ans=1;
    while(b){
        if(b&1)(ans*=a)%=p;
        (a*=a)%=p;
        b>>=1;
    }
    return ans;
}

int Tmp_Inv[N];
inline void inverse(int *A,int *Inv,int lenA){
    Inv[0]=qpow(A[0],p-2);
    for(int i=2;i<=lenA;i<<=1){
        MTT(A,Inv,Tmp_Inv,i,i);
        MTT(Tmp_Inv,Inv,Tmp_Inv,i,i);
        for(int j=0;j<i;++j)Inv[j]=(2ll*Inv[j]-Tmp_Inv[j]+p)%p;
    }
}

inline void differentiate(int *A,int *Dif,int len){
    for(int i=1;i<len;++i)Dif[i-1]=(ll)A[i]*i%p;
    Dif[len-1]=0;
}

inline void integrate(int *A,int *Int,int len){
    for(int i=len-1;i>=1;--i)Int[i]=(ll)A[i-1]*qpow(i,p-2)%p;
    Int[0]=0;
}

int Inv_ln[N];
inline void ln(int *A,int *Ln,int n){
    inverse(A,Inv_ln,n);
    differentiate(A,Ln,n); 
    MTT(Ln,Inv_ln,Ln,n,n);
    integrate(Ln,Ln,n<<1);
}

int tot;
int mu[N],pri[N];
bool mark[N];

inline void init(int n){
    mu[1]=1;
    for(int i=2;i<=n;++i){
        if(!mark[i])pri[++tot]=i,mu[i]=-1;
        for(int j=1,tmp;j<=tot&&(tmp=i*pri[j])<=n;++j){
            mark[tmp]=1;
            if(i%pri[j]==0){
                mu[tmp]=0;
                break;
            }
            mu[tmp]=-mu[i];
        }
    }
}

int F[N],G[N];
int main(){
    int n,m;read(n);read(p);
    for(int i=1;i<=n;++i)read(F[i]);F[0]=1;
    for(m=1;m<n;m<<=1);
    ln(F,F,m);
    for(int i=0;i<=n;++i)F[i]=(ll)F[i]*i%p;
    init(n);
    for(int i=1;i<=n;++i){
        for(int j=1;i*j<=n;++j){
            G[i*j]+=F[i]*mu[j];
        }
    }
    vector<int>ans;
    for(int i=1;i<=n;++i)if(G[i])ans.push_back(i);
    printf("%d
",ans.size());
    for(int i=0;i<ans.size();++i){
        if(i)putchar(' ');
        printf("%d",ans[i]);
    }
}