Description：

给出一个数组 $f[]$ 模 $p$ 的值， $f[i]$ 表示使用 $S$ 内的数不考虑顺序的构成 $i$ 的方案数，求 $S$ 。

$1\leqslant n\leqslant 2^{18}$

Solution：

设 $v_i$ 表示数字 $i$ 是否在集合中，那么 $f$ 的生成函数是： $$ F(x)=\prod_{i=1}^n\Big(\frac1{1-x^i}\Big)^{v_i} $$ 看到指数有未知量很不好处理，于是两边求对数： $$ \begin{align} \ln F(x)&=\sum_{i=1}^nv_i\ln\Big(\frac1{1-x^i}\Big)\\ &=\sum_{i=1}^nv_i\sum_{k\geqslant 0}\frac1kx^{ik}\\ &=\sum_{T=0}^nx^T\sum_{i|T}v_i\frac iT \end{align} $$ 于是我们莫比乌斯反演就可以求出 $v_i$ 了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

inline int rd()

{

	register int res = 0,f = 1;register char c = getchar();

	while(!isdigit(c)){if(c == '-')f = -1;c = getchar();}

	while(isdigit(c))res = (res << 1) + (res << 3) + c - '0',c = getchar();

	return res * f;

}

#define MAXN 524300

int n,p;

int f[MAXN],g[MAXN];

int len[MAXN];

typedef long long ll;

const long double PI = acos(-1);

struct comp

{

    long double r,i;

    comp(long double r_ = 0,long double i_ = 0){r = r_;i = i_;}

};

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

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

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

comp conj(comp A) {return (comp){A.r,-A.i};}

int power(int a,int b)

{

	int res = 1;

	for(;b > 0;b = b >> 1,a = 1ll * a * a % p)if(b & 1)res = 1ll * res * a % p;

	return res;

}

int inver(int a){return power(a,p - 2);}

int rev[MAXN];

void FFT(comp f[],int l)

{

    for(int i = 0;i < l;++i)

    {

        if(i < rev[i])swap(f[i],f[rev[i]]);

    }

    for(int i = 2;i <= l;i = i << 1)

    {

        comp wn = (comp){cos(2 * PI / i),sin(2 * PI / i)};

        for(int j = 0;j < l;j += i)

        {

            comp w = (comp){1,0};

            for(int k = j;k < j + i / 2;++k)

            {

                comp u = f[k],t = w * f[k + i / 2];

                f[k] = u + t;

                f[k + i / 2] = u - t;

                w = w * wn;

            }

        }

    }

    return;

}

comp a[MAXN],b[MAXN],r1[MAXN],r2[MAXN],r3[MAXN],r4[MAXN];

void mul(int x[],int y[],int z[],int l)

{

	for(int i = 0;i < l;++i)a[i] = comp(x[i] & 32767,x[i] >> 15);

	for(int i = 0;i < l;++i)b[i] = comp(y[i] & 32767,y[i] >> 15);

	for(int i = 0;i < (l >> 1);++i)a[i + (l >> 1)] = b[i + (l >> 1)] = comp(0,0);

	FFT(a,l);FFT(b,l);

	for(int i = 0;i < l;++i)

	{

		int j = (l - i) & (l - 1);

		comp da,db,dc,dd;

		da = (a[i] + conj(a[j])) * comp(0.5,0);

		db = (a[i] - conj(a[j])) * comp(0,-0.5);

		dc = (b[i] + conj(b[j])) * comp(0.5,0);

		dd = (b[i] - conj(b[j])) * comp(0,-0.5);

		r1[j] = da * dc;

		r2[j] = da * dd;

		r3[j] = db * dc;

		r4[j] = db * dd;

	}

	for(int i = 0;i < l;++i)a[i] = r1[i] + r2[i] * comp(0,1);

	for(int i = 0;i < l;++i)b[i] = r3[i] + r4[i] * comp(0,1);

	FFT(a,l);FFT(b,l);

	for(int i = 0;i < l;++i)

	{

		int da = (ll)(a[i].r / l + 0.5) % p;

		int db = (ll)(a[i].i / l + 0.5) % p;

		int dc = (ll)(b[i].r / l + 0.5) % p;

		int dd = (ll)(b[i].i / l + 0.5) % p;

		z[i] = (da + ((ll)(db + dc) << 15) + ((ll)dd << 30)) % p;

	}

	return;

}

int k1[MAXN],k2[MAXN],k3[MAXN],k4[MAXN],k5[MAXN];

void poly_inv(int a[],int b[],int deg)

{

	if(deg == 1){b[0] = inver(a[0]);return;}

	poly_inv(a,b,deg >> 1);

	for(int i = 0;i < (deg << 1);++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len[deg << 1] - 1));

	mul(a,b,k3,deg << 1);mul(k3,b,k4,deg << 1);

	for(int i = 0;i < deg;++i)b[i] = (2 * b[i] % p - k4[i] + p) % p,b[i + deg] = 0;

	return;

}

void poly_ln(int a[],int b[],int deg)

{

	for(int i = 1;i < deg;++i)k1[i - 1] = 1ll * a[i] * i % p;k1[deg - 1] = 0;

	poly_inv(a,k2,deg);

	for(int i = 0;i < (deg << 1);++i)rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len[deg << 1] - 1));

	mul(k1,k2,k5,deg << 1);

	for(int i = 1;i < deg;++i)b[i] = 1ll * k5[i - 1] * inver(i) % p;

	b[0] = 0;

	return;

}

int main()

{

	n = rd();p = rd();

	for(int i = 1;i <= n;++i)f[i] = rd();

	int kn = 1;while(kn <= n)kn = kn << 1,len[kn] = len[kn >> 1] + 1;

	len[kn << 1] = len[kn] + 1;

	f[0] = 1;

	poly_ln(f,g,kn);

	for(int i = 1;i <= n;++i)g[i] = 1ll * g[i] * i % p;

	for(int i = 1;i <= n;++i)

		for(int j = i + i;j <= n;j += i)g[j] = (g[j] - g[i] + p) % p;

    int cnt = 0;

	for(int i = 1;i <= n;++i)if(g[i])++cnt;

	cout << cnt << endl;

	for(int i = 1;i <= n;++i)if(g[i])printf("%d ",i);

	cout << endl;

    return 0;

}