P3784 [SDOI2017] 遗忘的集合

题目描述

https://www.luogu.com.cn/problem/P3784

简要题意：假设现在有一个集合 $S$，$S$ 中只含若干个 $[1,n]$ 的正整数，我们令 $f(k)$ 为把 $k$ 表示成 $S$ 中的元素的和的方案数，每个数都能用无限次，每种方案都是无序的，现在给定 $f(1)$ 到 $f(n)$ 这 $n$ 个模 $p$ 的值，求构造一个字典序最小的集合 $S$，保证 $p$ 是素数，但不一定是 $998244353$

$n\le 2^{18}$

Solution

我们不妨令 $f(x)$ 的生成函数为 $F(x)$，根据一些经典结论容易得到 $F(x)=\prod_{a_i\in S}\frac{1}{1-x^{a_i}}=\exp(\sum_{a_i\in S}\ln\frac{1}{1-x^{a_i}})=\exp(\sum_{a_i\in S}\sum_{n=1}\frac{x^{na_i}}{n})$

两边取一下 $\ln$，能够得到 $\ln F(x)=\sum_{a_i\in S}\sum_{n=1}\frac{x^{na_i}}{n}$，我们令 $G(x)$ 表示 $S$ 的生成函数，含义为如果 $S$ 有 $a_i$ 这个元素，那么 $x^{a_i}$ 的系数为 $1$，否则为 $0$，那么能够得到 $f(k)=\frac{1}{k}\sum_{d|k}\frac{d}{k}g(d)$，容易发现这是一个莫比乌斯反演的形式，那么答案一定是唯一的

唯一的难点在于任意模数多项式求逆，用三模NTT写确实慢，时间复杂度 $O(n\log n)$

#include <iostream>
#include <vector>
#include <algorithm>
#define maxn 270010
#define ll long long
using namespace std;

int p;
ll pow_mod(ll x, ll n, int p) {
    ll s = 1;
    for (; n; n >>= 1, x = x * x % p)
        if (n & 1) s = s * x % p;
    return s;
}

const int mod1 = 998244353, mod2 = 1004535809, mod3 = 469762049, G = 3;
const ll mod12 = 1ll * mod1 * mod2;
const int inv1 = pow_mod(mod1, mod2 - 2, mod2), inv2 = pow_mod(mod12 % mod3, mod3 - 2, mod3);
struct Int {
    int x, y, z;

    Int() { x = y = z = 0; }
    Int(int x, int y, int z) : x(x), y(y), z(z) {}
    Int(int v) : x(v), y(v), z(v) {}

    static inline int add(int x, int y, int p) { return (x += y) >= p ? x - p : x; }
    
    inline friend Int operator + (const Int &u, const Int &v) {
        return Int(add(u.x, v.x, mod1), add(u.y, v.y, mod2), add(u.z, v.z, mod3));
    }
    inline friend Int operator - (const Int &u, const Int &v) {
        return Int(add(u.x, mod1 - v.x, mod1), add(u.y, mod2 - v.y, mod2), add(u.z, mod3 - v.z, mod3));
    }
    inline friend Int operator * (const Int &u, const Int &v) {
        return Int(1ll * u.x * v.x % mod1, 1ll * u.y * v.y % mod2, 1ll * u.z * v.z % mod3);
    }
    
    inline int get() const {
        ll v = 1ll * add(y, mod2 - x, mod2) * inv1 % mod2 * mod1 + x;
        return (1ll * add(z, mod3 - v % mod3, mod3) * inv2 % mod3 * (mod12 % p) % p + v) % p;
    }
};

typedef vector<Int> Poly;
#define len(a) ((int) a.size())
namespace Pol {
    Poly operator - (const Int &v, const Poly &a) {
        Poly res(a); 
        for (int i = 0; i < len(res); ++i) res[i] = p - res[i]; 
        res[0] = res[0] + v; return res;
    }
    Poly operator - (const Poly &a, const Int &v) {
        Poly res(a); res[0] = res[0] - v; return res;
    }
    Poly operator * (const Poly &a, const Int &v) {
        Poly res(a);
        for (int i = 0; i < len(res) ; ++i) res[i] = res[i] * v; 
        return res;
    }
    
    const int N = 4200000;
    const int Gi1 = pow_mod(G, mod1 - 2, mod1);
    const int Gi2 = pow_mod(G, mod2 - 2, mod2);
    const int Gi3 = pow_mod(G, mod3 - 2, mod3);

    int P[N];
    void init_P(int n) {
        int l = 0; while ((1 << l) < n) ++l;
        for (int i = 0; i < n; ++i) P[i] = (P[i >> 1] >> 1) | ((i & 1) << l - 1);
    }
    void NTT(Poly &a, int type) {
        static Int w[N]; int n = len(a); 
        for (int i = 0; i < n; ++i) if (i < P[i]) swap(a[i], a[P[i]]);
        for (int i = 2, m = 1; i <= n; m = i, i *= 2) {
            Int wn = Int(pow_mod(type > 0 ? G : Gi1, (mod1 - 1) / i, mod1),
                        pow_mod(type > 0 ? G : Gi2, (mod2 - 1) / i, mod2),
                        pow_mod(type > 0 ? G : Gi3, (mod3 - 1) / i, mod3));
            w[0] = Int(1); for (int j = 1; j < m; ++j) w[j] = wn * w[j - 1];
            for (int j = 0; j < n; j += i) 
                for (int k = 0; k < m; ++k) {
                    Int t1 = a[j + k], t2 = a[j + k + m] * w[k];
                    a[j + k] = t1 + t2;
                    a[j + k + m] = t1 - t2;
                }
        }
        if (type < 0) {
            Int inv = Int(pow_mod(n, mod1 - 2, mod1), pow_mod(n, mod2 - 2, mod2), pow_mod(n, mod3 - 2, mod3));
            for (int i = 0; i < n; ++i) a[i] = a[i] * inv;
        }
    }
    Poly restore(Poly &a) {
        for (int i = 0; i < len(a); ++i) a[i] = a[i].get();
        return a;
    }
    Poly operator * (const Poly &A, const Poly &B) {
        int n = 1, n1 = len(A), n2 = len(B); while (n < n1 + n2 - 1) n <<= 1; init_P(n);
        Poly a(n), b(n); copy_n(A.begin(), n1, a.begin()); copy_n(B.begin(), n2, b.begin());
        NTT(a, 1); NTT(b, 1);
        for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
        NTT(a, -1); Poly ans(n1 + n2 - 1);
        for (int i = 0; i < n1 + n2 - 1; ++i) ans[i] = a[i].get();
        return ans;
    }
    Poly Der(const Poly &a) {
        Poly res(a);
        for (int i = 0; i < len(a) - 1; ++i) res[i] = Int(i + 1) * res[i + 1];
        res[len(a) - 1] = Int(0); restore(res); return res;
    }
    Poly Inte(const Poly &a) {
        static int inv[N];
        if (!inv[1]) {
            inv[1] = 1; 
            for (int i = 2; i < N; ++i) inv[i] = 1ll * (p - p / i) * inv[p % i] % p;
        }
        Poly res(a); res.resize(len(a) + 1);
        for (int i = len(a); i; --i) res[i] = res[i - 1] * Int(inv[i]);
        res[0] = Int(0); restore(res); return res;
    }
    Poly Inv(const Poly &a) {
        Poly res(1, Int(pow_mod(a[0].get(), p - 2, p)));
        int n = 1; while (n < len(a)) n <<= 1; 
        for (int k = 2; k <= n; k <<= 1) {
            Poly t(a); t.resize(k); 
            t = t * res; t.resize(k); 
            res = res * (2 - t); res.resize(k);
        } res.resize(len(a)); return res;
    }
    Poly Ln(const Poly &a) {
        Poly res = Inte(Der(a) * Inv(a));
        res.resize(len(a)); return res;
    }
}  // namespace Pol

int mu[maxn], pri[maxn], cnt; bool isp[maxn];
void init_isp(int n) {
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!isp[i]) pri[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <= cnt && i * pri[j] <= n; ++j) {
            isp[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            mu[i * pri[j]] = -mu[i];
        }
    }
}

int n;

ll f[maxn], g[maxn];
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr); cout.tie(nullptr);

    cin >> n >> p; Poly A(n + 1); A[0] = Int(1); init_isp(n); 
    for (int i = 1, x; i <= n; ++i) cin >> x, A[i] = Int(x);
    A = Pol::Ln(A); for (int i = 1; i <= n; ++i) f[i] = A[i].get(), f[i] = f[i] * i % p; 
    for (int i = 1; i <= n; ++i)
        for (int j = i; j <= n; j += i) g[j] = (g[j] + mu[j / i] * f[i]) % p;
    vector<int> ans;
    for (int i = 1; i <= n; ++i) if ((g[i] + p) % p) ans.push_back(i);
    cout << ans.size() << "\n";
    for (auto t : ans) cout << t << " "; 
    return 0;
}