SDOI2017 遗忘的集合

亦凉 2021-11-01 09:38 305阅读 0赞

题目描述：

[luogu][]

bzoj

题解：

生成函数+多项式ln（+反演？）

首先如果我们已知$S$，那$S$中$i$的生成函数就是$1+x^i+x^\{2i\}+……$。

整理一下就是$ \\frac\{1\} \{1-x^i\}$。

所以设$S$的生成函数为$F(x)$，则$F(x)= \\prod \_\{i \\in S\} \\frac\{1\}\{1-x^i\} $

由于我们并不知道有哪些是答案，我们可以引入一个系数$a\_i \\in \{0,1\}$，表示数$i$不取/取。

所以$F(x) = \\prod\_\{i \\ge 1\}(\\frac\{1\}\{1-x^i\})^\{a\_i\}$

其实都明白$i$是有上界的，因为整个式子都是模$x^\{n+1\}$，所以上界为$n$。

然后对上式取**负对数**，得$-ln(F(x)) = \\sum\_\{i \\ge 1\} a\_i \* ln(1-x^i)$。

同时$ln(1-x^i)=\\sum\_\{j \\ge 1\}- \\frac\{x^\{ij\}\}\{j\}$

然后左右去掉负号，得$ln(F(x)) = \\sum\_\{i \\ge 1\} a\_i \\sum\_\{j \\ge 1\} \\frac\{x^\{ij\}\}\{j\}$

把带$x$的提前，得$ln(F(x)) = \\sum\_\{T \\ge 1\}x^T \\sum\_\{i | T\} \\frac\{i\*a\_i\}\{T\}$

所以把$F(x)$取个$ln$，然后系数就可以$O(nlnn)$筛出了。

代码：

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

#include<cmath>
    #include<cstdio>
    #include<cstring>
    #include<algorithm>
    using namespace std;
    typedef long long ll;
    const int N = (1<<21)+50;
    const int M = 32768;
    const long double Pi = acos(-1.0);
    template<typename T>
    inline void read(T&x)
    {
        T f = 1,c = 0;char ch=getchar();
        while(ch<'0'||ch>'9'){
         if(ch=='-')f=-1;ch=getchar();}
        while(ch>='0'&&ch<='9'){c=c*10+ch-'0';ch=getchar();}
        x = f*c;
    }
    int n,MOD;
    void Mod(int&x){
         if(x>=MOD)x-=MOD;}
    int fastpow(int x,int y)
    {
        int ret = 1;
        while(y)
        {
            if(y&1)ret=1ll*ret*x%MOD;
            x=1ll*x*x%MOD;y>>=1;
        }
        return ret;
    }
    struct cp
    {
        long double x,y;
        cp(){}
        cp(long double x,long double y):x(x),y(y){}
        cp operator + (const cp&a)const{
         return cp(x+a.x,y+a.y);}
        cp operator - (const cp&a)const{
         return cp(x-a.x,y-a.y);}
        cp operator * (const cp&a)const{
         return cp(x*a.x-y*a.y,x*a.y+y*a.x);}
    };
    int to[N],lim,LL[N],L;
    void fft(cp*a,int len,int k)
    {
        for(int i=0;i<len;++i)
            if(i<to[i])swap(a[i],a[to[i]]);
        for(int i=1;i<len;i<<=1)
        {
            cp w0(cos(Pi/i),k*sin(Pi/i));
            for(int j=0;j<len;j+=(i<<1))
            {
                cp w(1,0);
                for(int o=0;o<i;o++,w=w*w0)
                {
                    cp w1 = a[j+o],w2 = w*a[j+o+i];
                    a[j+o] = w1+w2;
                    a[j+o+i] = w1-w2;
                }
            }
        }
        if(k==-1)
            for(int i=0;i<len;++i)a[i].x/=len,a[i].y/=len;
    }
    cp a[N],b[N],c[N],d[N];
    void get_lim(int len)
    {
    //    lim = 1,L = 0;
    //    while(lim<=len)lim<<=1,L++;
        lim = len;L=LL[len];
        for(int i=1;i<lim;++i)to[i]=((to[i>>1]>>1)|((i&1)<<(L-1)));
    }
    void mtt(int*F,int*G,int*H,int len)
    {
        get_lim(len<<1);
        for(register int i=0;i<lim;++i)a[i]=b[i]=cp(0,0);
        for(register int i=0;i<len;++i)
            a[i]=cp(F[i]&(M-1),F[i]>>15),b[i]=cp(G[i]&(M-1),G[i]>>15);
        fft(a,lim,1),fft(b,lim,1);
        for(register int i=0;i<lim;++i)
        {
            int j = (lim-i)&(lim-1);
            c[j]=cp(0.5*(a[i].x+a[j].x),0.5*(a[i].y-a[j].y))*b[i];
            d[j]=cp(0.5*(a[i].y+a[j].y),0.5*(a[j].x-a[i].x))*b[i];
        }
        fft(c,lim,1),fft(d,lim,1);
        for(register int i=0;i<lim;++i)
        {
            int kaa = (ll)(c[i].x/lim+0.5)%MOD;
            int kab = (ll)(c[i].y/lim+0.5)%MOD;
            int kba = (ll)(d[i].x/lim+0.5)%MOD;
            int kbb = (ll)(d[i].y/lim+0.5)%MOD;
            Mod(H[i] = ((ll)kaa+((ll)(kab+kba)<<15)%MOD+((ll)kbb<<30)%MOD)%MOD+MOD);
        }
    }
    int H[N],ny[N],T[N];
    void get_inv(int len,int*F,int*G)
    {
        if(len==1){G[0]=fastpow(F[0],MOD-2);return ;}
        get_inv(len>>1,F,G),mtt(G,G,T,len>>1),mtt(T,F,H,len);
        for(register int i=0;i<len;++i)Mod(G[i]=G[i]*2%MOD-H[i]+MOD);
    }
    void up(int*a,int len)
    {
        for(register int i=1;i<len;++i)a[i-1]=1ll*a[i]*i%MOD;
        a[len-1] = 0;
    }
    void down(int*a,int len)
    {
        for(register int i=len-1;i;--i)a[i]=1ll*a[i-1]*ny[i]%MOD;
        a[0]=0;
    }
    void get_ln(int len,int*F,int*G)
    {
        get_inv(len,F,G);up(F,len);
        mtt(F,G,G,len);
        down(G,len);
    }
    int F[N],G[N];
    int main()
    {
        read(n),read(MOD);F[0]=1;
        lim=LL[2]=1;while(lim<=n)lim<<=1,LL[lim<<1]=LL[lim]+1;
        ny[1]=1;
        for(register int i=2;i<(lim<<1);++i)ny[i]=1ll*(MOD-MOD/i)*ny[MOD%i]%MOD;
        for(register int i=1;i<=n;++i)
            read(F[i]),F[i]%=MOD;
        get_ln(lim,F,G);
        int cnt = 0;
        for(register int i=1;i<=n;++i)G[i]=1ll*G[i]*i%MOD;
        for(register int i=1;i<=n;++i)
            for(register int j=i+i;j<=n;j+=i)Mod(G[j]+=MOD-G[i]);
        for(register int i=1;i<=n;++i)if(G[i])cnt++;
        printf("%d\n",cnt);
        for(register int i=1;i<=n;++i)
            if(G[i])printf("%d ",i);
        puts("");
        return 0;
    }

（我还重新学了一遍MTT）

转载于:https://www.cnblogs.com/LiGuanlin1124/p/11050943.html

[luogu]: https://www.luogu.org/problemnew/show/P3784
[ContractedBlock.gif]: https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif
[ExpandedBlockStart.gif]: /images/20211101/fd4adc7105a740b1949676bd3abb051d.png