萌新刚学多项式,25pts求助
查看原帖
萌新刚学多项式,25pts求助
406941
Register_int-std=c++14楼主2022/8/15 15:10

rt.不知道是什么玄学问题……

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int MAXN = 1 << 20;
const ll moda = 998244353, modb = 1004535809, modc = 469762049;

int p[MAXN], mu[MAXN], tot;
bool vis[MAXN];

inline 
void init(int n) {
	mu[1] = 1;
	for (int i = 2; i <= n; i++) {
		if (!vis[i]) p[++tot] = i, mu[i] = -1;
		for (int j = 1; j <= tot; j++) {
			if (i * p[j] > n) break;
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) break;
			mu[i * p[j]] = -mu[i];
		}
	}
}

inline 
ll qpow(ll b, ll p, ll m) {
    ll res = 1;
    while (p) {
        if (p & 1) res = res * b % m;
        b = b * b % m, p >>= 1;
    }
    return res;
}

const ll modab = moda * modb, inva = qpow(moda, modb - 2, modb), invab = qpow(modab % modc, modc - 2, modc);

struct num {
	ll a, b, c;
	num() {}
	num(ll t) : a(t), b(t), c(t) {}
	num(ll a, ll b, ll c) : a(a), b(b), c(c) {}
	num reduce() { return num((a % moda + moda) % moda, (b % modb + modb) % modb, (c % modc + modc) % modc); }
	num operator + (const num &rhs) const { return num(a + rhs.a, b + rhs.b, c + rhs.c).reduce(); }
	num operator - (const num &rhs) const { return num(a - rhs.a, b - rhs.b, c - rhs.c).reduce(); }
	num operator * (const num &rhs) const { return num(a * rhs.a, b * rhs.b, c * rhs.c).reduce(); }
	ll get(ll mod) {
		ll t = (b - a + modb) % modb * inva % modb * moda + a;
		return ((c - t % modc + modc) % modc * invab % modc * (modab % mod) % mod + t) % mod;
	}
};

inline 
num inv(num x) {
	return num(qpow(x.a, moda - 2, moda), qpow(x.b, modb - 2, modb), qpow(x.c, modc - 2, modc));
} 

int rev[MAXN];

num w[MAXN];

inline 
int getrev(int n) {
	int l = 1;
	while (l < n << 1) l <<= 1;
	for (int i = 1; i < l; i++) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? l >> 1 : 0);
	w[0] = num(1);
	num t = num(qpow(3, (moda - 1) / l, moda), qpow(3, (modb - 1) / l, modb), qpow(3, (modc - 1) / l, modc));
	for (int i = 1; i <= l; i++) w[i] = w[i - 1] * t;
	return l;
}

inline 
void ntt(num *f, int n, int t) {
	for (int i = 0; i < n; i++) {
		if (i < rev[i]) swap(f[i], f[rev[i]]);
	}
	for (int i = 1; i < n; i <<= 1) {
		int p = n / i >> 1;
        for (int j = 0; j < n; j += i << 1) {
            for (int k = j; k < i + j; k++) {
                num x = (t ? w[n - p * (k - j)] : w[p * (k - j)]) * f[i + k];
                f[i + k] = f[k] - x, f[k] = f[k] + x;
            }
        }
    }
    if (t) {
    	num q = num(qpow(n, moda - 2, moda), qpow(n, modb - 2, modb), qpow(n, modc - 2, modc));
        for (int i = 0; i < n; i++) f[i] = f[i] * q;
	}
}

num f1[MAXN];

void finv(num *f, num *g, int n) {
	if (n == 1) return g[0] = inv(f[0]), void();
	finv(f, g, n + 1 >> 1);
	int l = getrev(n);
	for (int i = 0; i < n; i++) f1[i] = f[i];
	for (int i = n; i < l; i++) f1[i] = num(0);
	ntt(f1, l, 0), ntt(g, l, 0);
	for (int i = 0; i < l; i++) g[i] = (num(2) - f1[i] * g[i]) * g[i];
	ntt(g, l, 1);
	for (int i = n; i < l; i++) g[i] = 0;
}

inline 
void drt(num *f, num *g, int n) {
	for (int i = 1; i < n; i++) g[i - 1] = num(i) * f[i];
	g[n - 1] = num(0);
}

inline 
void igt(num *f, num *g, int n) {
	g[0] = num(0);
	for (int i = 1; i < n; i++) g[i] = f[i - 1] * inv(i);
}

num f2[MAXN], g2[MAXN];

inline 
void fln(num *f, num *g, int n) {
	drt(f, f2, n), finv(f, g2, n);
	int l = getrev(n);
	ntt(f2, l, 0), ntt(g2, l, 0);
	for (int i = 0; i < l; i++) f2[i] = f2[i] * g2[i];
	ntt(f2, l, 1), igt(f2, g, n);
}

num f[MAXN], g[MAXN];

int n, m, x;

int a[MAXN], ans;

int main() {
	scanf("%d%d", &n, &m);
	init(n);
	f[0] = num(1);
	for (int i = 1; i <= n; i++) scanf("%d", &x), f[i] = num(x % m);
	fln(f, g, n + 1);
	for (int i = 1; i <= n; i++) g[i] = g[i] * num(i);
	for (int i = 1; i <= n; i++) {
		for (int j = i; j <= n; j += i) a[j] += mu[i] * g[j / i].get(m);
	}
	for (int i = 1; i <= n; i++) {
		if (a[i]) ans++;
	}
	printf("%d\n", ans);
	for (int i = 1; i <= n; i++) {
		if (a[i]) printf("%d ", i);
	}
}
2022/8/15 15:10
加载中...