bzoj 4913

好毒瘤的一道题啊...

对每个$a_{i}\in S$,设$F(x)$为用$j$个$a_{i}$构造出$ja_{i}$的生成函数，那么$F(x)=\sum_{j=1}^{∞}x^{ja_{i}}$

根据这篇博客里的内容，可以求得：$F(x)=\frac{1}{1-x^{a_{i}}}$

设$t_{i}=[i\in S]$，则用$S$集合中任意个数表示出一个数的生成函数即为$G(x)=\prod_{i=1}^{n}[\frac{1}{1-x^{i}}]^{t_{i}}$

 接下来做一些推导：

两边取负对数，得：

$-lnG(x)=\sum_{i=1}^{n}t_{i}ln(1-x^{i})$

对$ln(1-x^{i})$求导得到：

$\frac{-ix^{i-1}}{1-x^{i}}$

把下半部分恢复成等比数列求和的形式：

$-ix^{i-1}\sum_{j=0}^{∞}x^{ij}$

把外面的系数移进去：

$-\sum_{j=0}^{∞}ix^{ij+(i-1)}$

然后积分：

$\int -\sum_{j=0}^{∞}ix^{ij+(i-1)}=-\sum_{j=1}^{∞}\frac{i}{ij+i}x^{ij+i}$

约分一下，就得到了：

$-\sum_{j=0}^{∞}\frac{1}{j+1}x^{ij+i}$

令$j=j+1$，有：

$-\sum_{j=1}^{∞}\frac{1}{j}x^{ij}$

对一个函数先求导再积分得到的就是原函数，因此有：

$ln(1-x^{i})=-\sum_{j=1}^{∞}\frac{1}{j}x^{ij}$

再带入原来的表达式：

$-lnG(x)=-\sum_{i=1}^{n}t_{i}\sum_{j=1}^{∞}\frac{1}{j}x^{ij}$

约去两边的负号，也即：

$lnG(x)=\sum_{i=1}^{n}t_{i}\sum_{j=1}^{∞}\frac{1}{j}x^{ij}$

考虑令$T=ij$，则$j=\frac{T}{i}$，也即：

$lnG(x)=\sum_{i=1}^{n}t_{i}\sum_{i|T}^{∞}\frac{i}{T}x^{T}$

 优先枚举$T$,得到：

$lnG(x)=\sum_{T=1}^{∞}\sum_{i|T}t_{i}\frac{i}{T}x^{T}$

如果设$lnG(x)$的第$i$项为$g(i)x^{i}$，则$g(i)=\sum_{d|i}t_{d}\frac{i}{T}$

可以用枚举倍数的思想求解，这样的复杂度是调和级数$O(nlnn)$（即每计算出一个$d_{i}$，就枚举$i$的所有倍数$T$减去$d_{i}\frac{i}{T}$的贡献即可）

（不要忘了$G(x)$已知，但是需要多项式黑科技，因为模数不好...）

（多项式黑科技精度问题可能被卡，用另一个黑科技就好）

代码：

#include <cstdio>
#include <cmath>
#include <algorithm>
#define ll long long
#define ld long double
#define maxn 100000
using namespace std;
const ld pi=acos(-1.0);
const ll siz=32767;
int n;
ll mode;
ll Ff[(1<<20)+5],Gg[(1<<20)+5];
ll ff[(1<<20)+5];
ll tempG[(1<<20)+5];
ll ans[(1<<20)+5];
ll gg[(1<<20)+5];
ll inv[(1<<20)+5];
int l[(1<<20)+5];
int Lim=1;
void init()
{
    inv[0]=inv[1]=1;
    for(int i=2;i<=(1<<20);i++)inv[i]=(mode-mode/i)*inv[mode%i]%mode;
}
ll pow_mul(ll x,ll y)
{
    ll ret=1;
    while(y)
    {
        if(y&1)ret=ret*x%mode;
        x=x*x%mode,y>>=1;
    }
    return ret;
}

struct cp
{
    ld x,y;
    cp (){}
    cp (ld a,ld b):x(a),y(b){}
    cp operator + (const cp&a)const
    {
        return cp(a.x+x,a.y+y);
    }
    cp operator - (const cp&a)const
    {
        return cp(x-a.x,y-a.y);
    }
    cp operator * (const cp&a)const
    {
        return cp(a.x*x-a.y*y,a.x*y+a.y*x);
    }
};
cp get_w(cp a)
{
    return cp(a.x,-a.y);
}
int to[(1<<20)+5];
void FFT(cp *a,int len,int k)
{
    for(register int i=0;i<len;i++)if(i<to[i])swap(a[i],a[to[i]]);
    for(register int i=1;i<len;i<<=1)
    {
        cp w0=cp(cos(pi/i),k*sin(pi/i));
        for(register int j=0;j<len;j+=(i<<1))
        {
            cp w=cp(1,0);
            for(int o=0;o<i;o++,w=w*w0)
            {
                cp w1=a[j+o],w2=a[j+o+i]*w;
                a[j+o]=w1+w2,a[j+o+i]=w1-w2;
            }
        }
    }
}
cp A[(1<<20)+5],B[(1<<20)+5],C[(1<<20)+5],D[(1<<20)+5];
ll ret[(1<<20)+5];
inline void MTT(ll *f,ll *g,ll *re,int len)
{
    for(int i=0;i<len;i++)to[i]=((to[i>>1]>>1)|((i&1)<<(l[len]-1)));
    for(int i=0;i<len;i++)A[i]=B[i]=cp(0,0);
    for(int i=0;i<(len>>1);i++)A[i]=cp((ld)(f[i]&siz),(ld)(f[i]>>15));
    for(int i=0;i<(len>>1);i++)B[i]=cp((ld)(g[i]&siz),(ld)(g[i]>>15));
    FFT(A,len,1),FFT(B,len,1);
    for(int i=0;i<len;i++)
    {
        int j=(len-i)&(len-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,len,1),FFT(D,len,1);
    for(int i=0;i<len;i++)
    {
        ll tempA=(ll)(C[i].x/len+0.5)%mode;
        ll tempa=(ll)(C[i].y/len+0.5)%mode;
        ll tempB=(ll)(D[i].x/len+0.5)%mode;
        ll tempb=(ll)(D[i].y/len+0.5)%mode;
        re[i]=(ll)(tempA+mode+(ll)((ll)(tempa+tempB)<<15)%mode+mode+(ll)((ll)tempb<<30)%mode+mode)%mode;
    }
}
void get_inv(ll *f,ll *g,int dep)
{
    if(dep==1)
    {
        g[0]=pow_mul(f[0],mode-2);
        return;
    }
    int nxt=(dep)>>1;
    get_inv(f,g,nxt);
    MTT(g,g,tempG,dep);
    MTT(f,tempG,ret,dep<<1);
    for(register int i=0;i<dep;i++)g[i]=(2*g[i]-ret[i]+mode)%mode;
}
void get_ln(ll *f,ll *g,int dep)
{
    get_inv(f,Gg,dep);
    for(int i=1;i<dep;i++)ff[i-1]=f[i]*1ll*i%mode;
    MTT(Gg,ff,ret,dep<<1);
    for(int i=1;i<dep;i++)g[i]=ret[i-1]*pow_mul(i,mode-2)%mode;
}
template <typename T>inline void read(T &x)
{
    int f=1;
    x=0;
    char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    x*=f;
}
void print(int x)
{
    if(!x)return;
    print(x/10);
    putchar('0'+x%10);
}
int main()
{
    read(n),read(mode);
    for(register int i=1;i<=n;i++)read(Ff[i]);
    int Lim=1;
    while(Lim<=n)Lim<<=1,l[Lim]=l[Lim>>1]+1;
    l[Lim<<1]=l[Lim]+1;
    init();
    Ff[0]=1;
    get_ln(Ff,gg,Lim);
    int cnt=0;
    for(register int i=1;i<=n;i++)
    {
        ans[i]=gg[i];
        if(ans[i])cnt++;
        for(register int j=i+i;j<=n;j+=i)gg[j]=(gg[j]-ans[i]*i%mode*inv[j]%mode+mode)%mode;    
    }
    printf("%d\n",cnt);
    for(register int i=1;i<=n;i++)if(ans[i])print(i),putchar(' ');
    printf("\n");
    return 0;
}