[LOJ2271] [SDOI2017] 遗忘的集合

题目链接

LOJ：https://loj.ac/problem/2271

洛谷：https://www.luogu.org/problemnew/show/P3784

BZOJ太伤身体死活卡不过还是算了吧...

Solution

为啥窝洛谷(rk4) ( m BZOJ)死活跑不过啊...

技不如人，肝败吓疯...

---

题目并不难

设(a_i)表示(i)有没有出现在集合中，这是我们要求的答案。

那么把背包写成生成函数就是：

[prod_{i=1}^{n}left(sum_{j=0}^{+infty}x^{ij} ight)^{a_i}=prod_{i=1}^{n}left(frac{1}{1-x^i} ight)^{a_i}=F(x) ]

两边取对数：

[sum_{i=1}^{n}a_iln frac{1}{1-x^i}=ln F(x) ]

注意到：

[ln G(x)=int frac{G'(x)}{G(x)}dx ]

把上面那个分式化一下就是：

[ln frac{1}{1-x^i}=int frac{ix^{i-1}}{1-x^i} dx=iint sum_{j=0}^{+infty}x^{ij+i-1}dx=sum_{j=1}^{+infty}frac{1}{j}x^{ij} ]

带一下：

[sum_{i=1}^{n}a_isum_{j=1}^{+infty}frac{1}{j}x^{ij} =ln F(x) ]

我们枚举(ij)：

[sum_{T=1}^{+infty}x^Tsum_{i|T}a_icdot frac{i}{T}=ln F(x)\ sum_{T=1}^{+infty}x^Tsum_{i|T}a_icdot i=Tcdot ln F(x) ]

注意到左边是一个莫反的形式，我们对(F(x))求(ln)，然后把系数搞出来，在套莫反公式就好了。

复杂度(O(nlog n))。

#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

const int maxn = 6e5+10;
const int inf = 1e9;
const lf eps = 1e-8;
const lf pi = acos(-1);
const int all = (1<<15)-1;

struct cp {
    lf r,i;
    cp () {}
    cp (lf _r,lf _i) {r=_r,i=_i;}
    cp operator + (const cp &rhs) const {return cp(r+rhs.r,i+rhs.i);}
    cp operator - (const cp &rhs) const {return cp(r-rhs.r,i-rhs.i);}
    cp operator * (const cp &rhs) const {return cp(r*rhs.r-i*rhs.i,i*rhs.r+r*rhs.i);}
    cp operator / (const int &rhs) const {return cp(r/(lf)rhs,i/(lf)rhs);}
    void print() {printf("%lf %lf
",r,i);}
}w[maxn],g[5][maxn];

cp conj(cp x) {return cp(x.r,-x.i);}

int n,m,N,bit,mxn,mod;
int pos[maxn],a[maxn],b[maxn],c[maxn],inv[maxn],f[maxn],h[maxn];

void fft(cp *r) {
    for(int i=1;i<N;i++) if(pos[i]>i) swap(r[i],r[pos[i]]);
    for(int i=1,d=mxn>>1;i<N;i<<=1,d>>=1)
        for(int j=0;j<N;j+=i<<1) 
            for(int k=0;k<i;k++) {
                cp x=r[j+k],y=w[k*d]*r[i+j+k];
                r[j+k]=x+y,r[i+j+k]=x-y;
            }
}

void ntt_get(int len) {
    for(N=1,bit=0;N<=len;N<<=1,bit++);
    for(int i=1;i<N;i++) pos[i]=pos[i>>1]>>1|((i&1)<<(bit-1));
}

void poly_mul(int *r,int *s,int *t,int len) {
    ntt_get(len);
    for(int i=0;i<len;i++) g[0][i]=cp(r[i]&all,r[i]>>15),g[1][i]=cp(s[i]&all,s[i]>>15);
    for(int i=len;i<N;i++) g[0][i]=g[1][i]=cp(0,0);
    fft(g[0]),fft(g[1]);
    for(int i=0;i<N;i++) {
        int j=(N-i)&(N-1);
        g[2][j]=(g[0][i]+conj(g[0][j]))*cp(0.5,0)*g[1][i];
        g[3][j]=(g[0][i]-conj(g[0][j]))*cp(0,-0.5)*g[1][i];
    }fft(g[2]),fft(g[3]);
    for(int i=0;i<N;i++) g[2][i]=g[2][i]/N,g[3][i]=g[3][i]/N;
    for(int i=0;i<N;i++) {
        ll pp=g[2][i].r+0.5,x=g[2][i].i+0.5,y=g[3][i].r+0.5,z=g[3][i].i+0.5;
        t[i]=((pp%mod+(((x+y)%mod)<<15)+((z%mod)<<30))%mod+mod)%mod;
    }
}

int add(int x,int y) {return x+y>mod?x+y-mod:x+y;}
int del(int x,int y) {return x-y<0?x-y+mod:x-y;}
int mul(int x,int y) {return 1ll*x*y-1ll*x*y/mod*mod;}

int tmp[10][maxn];

int qpow(int aa,int x) {
    int res=1;
    for(;x;x>>=1,aa=mul(aa,aa)) if(x&1) res=mul(res,aa);
    return res;
}

void poly_inv(int *r,int *t,int len) {
    if(len==1) return t[0]=qpow(r[0],mod-2),void();
    poly_inv(r,t,len>>1);ntt_get(len);
    poly_mul(t,r,tmp[2],len);
    poly_mul(tmp[2],t,tmp[3],len);
    for(int i=0;i<len;i++) t[i]=del(add(t[i],t[i]),tmp[3][i]);
}

void poly_der(int *r,int *t,int len) {
    ntt_get(len);
    for(int i=1;i<len;i++) t[i-1]=mul(i,r[i]);
    for(int i=len-1;i<N;i++) t[i]=0;
}

void poly_int(int *r,int *t,int len) {
    ntt_get(len);
    for(int i=0;i<len;i++) t[i+1]=mul(inv[i+1],r[i]);t[0]=0;
    for(int i=len+1;i<N;i++) t[i]=0;
}

void poly_ln(int *r,int *t,int len) {
    poly_der(r,tmp[4],len);
    poly_inv(r,tmp[5],len);
    poly_mul(tmp[4],tmp[5],tmp[6],len);
    poly_int(tmp[6],t,len);
}

int pri[maxn],tot,mu[maxn],A[maxn],vis[maxn];

void sieve() {
    mu[1]=1;
    for(int i=2;i<=mxn;i++) {
        if(!vis[i]) mu[i]=-1,pri[++tot]=i;
        for(int j=1;j<=tot&&i*pri[j]<=mxn;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    }
}

int main() {
    read(n),read(mod);ntt_get(n<<1);inv[0]=inv[1]=1;
    for(mxn=1;mxn<N;mxn<<=1);
    for(int i=0;i<=mxn;i++) w[i]=cp(cos(2.0*pi*i/mxn),sin(2.0*pi*i/mxn));
    for(int i=2;i<N;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
    for(int i=1;i<=n;i++) read(f[i]);f[0]=1;
    poly_ln(f,h,N>>1);sieve();
    for(int i=1;i<=n;i++) h[i]=mul(h[i],i);
    for(int i=1;i<=n;i++)
        for(int j=i;j<=n;j+=i)
            A[j]+=h[j/i]*mu[i];
    int ans=0;
    for(int i=1;i<=n;i++) if(A[i]) ans++;
    write(ans);
    for(int i=1;i<=n;i++) if(A[i]) printf("%d ",i);puts("");
    return 0;
}