感觉这种东西 WA 了还是比较离奇的。

#include<bits/stdc++.h>
#define ri register int
#define ll long long
#define ld long double
using namespace std;
const int maxn=1e6+7;
inline ll qpow(ll a,ll b,ll p){
	ll ans=1;
	while(b){
		if(b&1)ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
int mod=1e9+7;
namespace PolyTimes{
	const ld pi=acos(-1);
	struct cpn{
		ld a,b;
		cpn(const ld &a=0,const ld &b=0):a(a),b(b){};
		friend cpn operator + (cpn &a,cpn &b){
			return (cpn){a.a+b.a,a.b+b.b};
		}
		friend cpn operator - (cpn &a,cpn &b){
			return (cpn){a.a-b.a,a.b-b.b};
		}
		friend cpn operator * (const cpn &a,const cpn &b){
			return (cpn){a.a*b.a-a.b*b.b,a.b*b.a+a.a*b.b};
		}
//		friend cpn operator / (const cpn &a,const ld x){
//			return (cpn){a.a/x,a.b/x};
//		}
//		friend cpn operator / (const cpn &a,const int x){
//			return (cpn){a.a/x,a.b/x};
//		}
	};
	cpn A[maxn],B[maxn],T[maxn];
	int rev[maxn];
	void calc_rev(int limit){
		for(ri i=1;i<limit;i++)rev[i]=(rev[i>>1]>>1)|(i&1?limit>>1:0);
	}
	void FFT(cpn f[],bool tp,int limit){
		for(ri i=0;i<limit;i++){
			if(i<rev[i])swap(f[i],f[rev[i]]);
		}
		for(ri len=2;len<=limit;len<<=1){
			int pa=len>>1;
			cpn dwg(cos(2*pi/len),sin(2*pi/len)*(tp?1:-1));
//			cout<<pi<<' '<<len<<' '<<dwg.a<<' '<<dwg.b<<endl;
			for(ri i=0;i<limit;i+=len){
				cpn bf(1,0);
				for(ri j=i;j<i+pa;j++){
					cpn t=bf*f[j+pa];
					f[j+pa]=f[j]-t;
					f[j]=f[j]+t;
					bf=bf*dwg;
				}
			}
		}
		if(tp==0){
			for(ri i=0;i<limit;i++){
				f[i].a/=limit;f[i].b/=limit;
			}
		}
	}
	void times(int f[],int g[],int limit){
		fill(A,A+limit,cpn(0,0));
		fill(B,B+limit,cpn(0,0));
		fill(T,T+limit,cpn(0,0));
		for(ri i=0;i<limit>>1;i++){
			A[i]=(cpn){f[i]&((1<<16)-1),0};
			B[i]=(cpn){f[i]>>16,0};
			T[i]=(cpn){g[i]&((1<<16)-1),g[i]>>16};
		}
		calc_rev(limit);
//		for(ri i=0;i<limit;i++)cout<<rev[i]<<' ';cout<<endl;
		FFT(A,1,limit);FFT(B,1,limit);FFT(T,1,limit);
//		for(ri i=0;i<limit>>1;i++)cout<<A[i].a<<'+'<<A[i].b<<'i'<<' ';cout<<endl;
//		for(ri i=0;i<limit>>1;i++)cout<<B[i].a<<'+'<<B[i].b<<'i'<<' ';cout<<endl;
//		for(ri i=0;i<limit>>1;i++)cout<<T[i].a<<'+'<<T[i].b<<'i'<<' ';cout<<endl;
		for(ri i=0;i<limit;i++){
			A[i]=A[i]*T[i];
			B[i]=B[i]*T[i];
		}
//		for(ri i=0;i<limit;i++)cout<<A[i].a<<'+'<<A[i].b<<'i'<<' ';cout<<endl;
//		for(ri i=0;i<limit;i++)cout<<B[i].a<<'+'<<B[i].b<<'i'<<' ';cout<<endl;
		FFT(A,0,limit);FFT(B,0,limit);
		for(ri i=0;i<limit;i++){
			ll a=(ll)(A[i].a+0.49),b=(ll)(A[i].b+0.49),c=(ll){B[i].a+0.49},d=(ll){B[i].b+0.49};
			a%=mod;b%=mod;c%=mod;d%=mod;
			f[i]=(a+((b+c)%mod*(1<<16))%mod+d*(1<<16)%mod*(1<<16)%mod)%mod;
		}
	}
	int f[maxn],g[maxn];
	void test(){
		int n,m;
		cin>>n>>m;
		cin>>mod;
		for(ri i=0;i<=n;i++){
			cin>>f[i];
		}
		for(ri i=0;i<=m;i++){
			cin>>g[i];
		}
		int limit=1;
		while(limit<=n+m)limit<<=1;
		times(f,g,limit<<1);
		for(ri i=0;i<=n+m;i++)cout<<f[i]<<' ';
		cout<<endl;
		exit(0);
	}
}
namespace INV{
	int g[2][maxn];
	inline void get_inv(int h[],int limit){
		int cur=0;
		memset(g[cur],0,sizeof(g[cur]));
		g[cur][0]=qpow(h[0],mod-2,mod);
		for(ri len=2;len<=limit;len<<=1){
			cur^=1;
			fill(g[cur],g[cur]+(len<<1),0);
			for(ri i=0;i<len;i++){
				g[cur][i]=g[cur^1][i];
			}
			PolyTimes::times(g[cur],h,len<<1);
			for(ri i=0;i<len;i++){
				g[cur][i]=mod-g[cur][i];
			}
			g[cur][0]=(g[cur][0]+2)%mod;
			PolyTimes::times(g[cur],g[cur^1],len<<1);
		}
		for(ri i=0;i<limit;i++){
			h[i]=g[cur][i];
		}
	}
	int f[maxn];
	void test(){
		int n;
		cin>>n;
		for(ri i=0;i<n;i++){
			cin>>f[i];
		}
		int limit=1;
		while(limit<n)limit<<=1;
		get_inv(f,limit);
		for(ri i=0;i<n;i++){
			cout<<f[i]<<' ';
		}
		cout<<endl;
		exit(0);
	}
}
namespace Ln{
	int inv[maxn];
	void init(int limit){
		for(ri i=1;i<=limit;i++)inv[i]=qpow(i,mod-2,mod);
	}
	inline void differential(int f[],int limit){
		for(ri i=0;i<limit-1;i++){
			f[i]=1ll*f[i+1]*(i+1)%mod;
		}
		f[limit-1]=0;
	}
	inline void integral(int f[],int limit){
		for(ri i=limit-1;i;i--){
			f[i]=1ll*f[i-1]*inv[i]%mod;
		}
		f[0]=0;
	}
	int g[maxn];
	inline void get_ln(int f[],int limit){
		for(ri i=0;i<limit;i++)g[i]=f[i];
		differential(g,limit);
		INV::get_inv(f,limit);
		PolyTimes::times(f,g,limit<<1);
		integral(f,limit);
	}
}
namespace Exp{
	int g[2][maxn];
	inline void get_exp(int f[],int limit){
		g[0][0]=1;
		g[0][1]=0;
		int cur=0;
		for(ri len=2;len<=limit;len<<=1){
			cur^=1;
			fill(g[cur],g[cur]+(len<<1),0);
			for(ri i=0;i<len;i++)g[cur][i]=g[cur^1][i];
			Ln::get_ln(g[cur^1],len);
			for(ri i=0;i<len;i++){
				g[cur^1][i]=f[i]-g[cur^1][i]+mod;
				if(g[cur^1][i]>=mod)g[cur^1][i]-=mod;
			}
			g[cur^1][0]=(g[cur^1][0]+1);
			if(g[cur^1][0]>=mod)g[cur^1][0]-=mod;
			PolyTimes::times(g[cur],g[cur^1],len<<1);
		}
		for(ri i=0;i<limit;i++)f[i]=g[cur][i];
	}
}
int n;
int f[maxn];
int mu[maxn];
bool np[maxn];
void seive(int limit){
	vector<int>pri;
	mu[1]=1;
	for(ri i=2;i<=limit;i++){
		if(!np[i]){
			pri.emplace_back(i);
			mu[i]=-1;
		}
		for(auto p:pri){
			if(i*p>limit)break;
			np[i*p]=1;
			if(i%p==0){
				mu[i*p]=0;
				break;
			}
			mu[i*p]=-mu[i];
		}
	}
}
int h[maxn];
signed main(){
	ios::sync_with_stdio(0);
//	INV::test();
	cin>>n>>mod;
	f[0]=1;
	for(ri i=1;i<=n;i++){
		cin>>f[i];
	}
	int limit=1;
	while(limit<=n)limit<<=1;
	Ln::init(limit);
	Ln::get_ln(f,limit);
//	for(ri i=0;i<limit;i++){
//		cout<<f[i]<<' ';
//	}
//	cout<<endl;
	for(ri i=0;i<limit;i++){
		f[i]=1ll*f[i]*i%mod;
	}
//	for(ri i=1;i<=n;i++)cout<<f[i]<<' ';cout<<endl;
	seive(n);
//	for(ri i=1;i<=n;i++)cout<<mu[i]<<' ';cout<<endl;
	for(ri i=1;i<=n;i++){
		for(ri j=1;i*j<=n;j++){
			h[i*j]=(h[i*j]+1ll*f[i]*mu[j]%mod)%mod;
		}
	}
	int cnt=0;
	for(ri i=1;i<=n;i++){
		if(h[i])cnt++;
	}
	cout<<cnt<<endl;
	for(ri i=1;i<=n;i++){
		if(h[i])cout<<i<<' ';
	}
	cout<<endl;
	return 0;
}
/*
9 998244353
1 2 2 3 4 5 6 7 8
*/