您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> C++ >> 關於C++

BZOJ 3992 [SDOI2015] 序列統計：NTT 快速冪優化卷積

編輯：關於C++

題目大意

給定集合 S（元素 < M 的非負整數），生成長度 N 的序列（元素皆屬於 S）。求：所有可能序列中，元素乘積模 M 等於 x 的序列數量，模 1,004,535,809。

1 ≤ N ≤ 10⁹，3 ≤ M ≤ 8000（M 為質數）
1 ≤ x ≤ M-1，|S| 為集合大小

核心觀察：模 M 乘法群與原根

M 為質數 → 非零剩餘類構成循環群，存在原根 g。令每個數 a = gᵏ，則乘積對應指數相加模 (M-1)。

問題轉化為：在指數域上做卷積。

動態規劃 + NTT 優化

定義 f[i][j] = 由 i 個數組成、乘積模 M 等於 gʲ 的方案數（i 為 2 的冪）。

轉移：f[i][j] = Σₖ f[i/2][k] × f[i/2][j-k] → 循環卷積，可用 NTT 加速。

N 不一定是 2 的冪 → 採用 快速冪思想：

預計算 f[1], f[2], f[4], f[8]...（每次 NTT 卷積平方）
按 N 的二進位分解組合結果

關鍵參數

模數 P = 1,004,535,809 = 479 × 2²¹ + 1（適合 NTT 的質數）
原根 g = 3（經驗值，可驗證）
最大長度：M-1 ≤ 7999 → NTT 長度取 2¹⁴ = 16384 足夠

完整程式碼 (C++)

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const int MOD = 1004535809;
const int G = 3;           // 原根
const int MAXM = 8005;
const int MAXN = 1 << 14;  // 16384 > 2*(M-1)

int n, m, target, s;
int a[MAXM], ind[MAXM], rev[MAXN];
int A[MAXN], B[MAXN];

// 快速冪
int qpow(int a, ll e, int mod = MOD) {
    ll r = 1, x = a;
    while (e) {
        if (e & 1) r = r * x % mod;
        x = x * x % mod;
        e >>= 1;
    }
    return (int)r;
}

// 尋找模 mod 的原根
int primitive_root(int mod) {
    vector<int> fac;
    int phi = mod - 1, x = phi;
    for (int i = 2; i * i <= x; ++i)
        if (x % i == 0) {
            fac.push_back(i);
            while (x % i == 0) x /= i;
        }
    if (x > 1) fac.push_back(x);
    for (int g = 2; g < mod; ++g) {
        bool ok = true;
        for (int p : fac)
            if (qpow(g, phi / p, mod) == 1) { ok = false; break; }
        if (ok) return g;
    }
    return -1;
}

// NTT 核心
void ntt(vector<int>& a, bool invert) {
    int n = a.size();
    for (int i = 1, j = 0; i < n; ++i) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }
    for (int len = 2; len <= n; len <<= 1) {
        int wlen = qpow(G, (MOD - 1) / len);
        if (invert) wlen = qpow(wlen, MOD - 2);
        for (int i = 0; i < n; i += len) {
            ll w = 1;
            for (int j = 0; j < len/2; ++j) {
                int u = a[i+j], v = (int)(a[i+j+len/2] * w % MOD);
                a[i+j] = u + v < MOD ? u + v : u + v - MOD;
                a[i+j+len/2] = u - v >= 0 ? u - v : u - v + MOD;
                w = w * wlen % MOD;
            }
        }
    }
    if (invert) {
        int inv_n = qpow(n, MOD - 2);
        for (int& x : a) x = (int)(x * (ll)inv_n % MOD);
    }
}

// 循環卷積 (模 M-1)
void multiply(vector<int>& a, vector<int>& b, vector<int>& res) {
    int sz = 1;
    while (sz < (int)a.size() + (int)b.size()) sz <<= 1;
    vector<int> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    fa.resize(sz); fb.resize(sz);
    ntt(fa, false); ntt(fb, false);
    for (int i = 0; i < sz; ++i) fa[i] = (int)(fa[i] * (ll)fb[i] % MOD);
    ntt(fa, true);
    res.assign(m, 0);
    for (int i = 0; i < sz; ++i)
        res[i % (m-1)] = (res[i % (m-1)] + fa[i]) % MOD;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    cin >> n >> m >> target >> s;
    int g = primitive_root(m);
    
    // 建立 ind: 數值 -> 指數
    ind[0] = -1;
    int cur = 1;
    for (int i = 0; i < m-1; ++i) {
        ind[cur] = i;
        cur = cur * g % m;
    }
    
    // 初始多項式 B: B[ind[x]] = 1 for x in S
    vector<int> B(m-1, 0);
    for (int i = 0; i < s; ++i) {
        int x; cin >> x;
        if (x != 0) B[ind[x]] = 1;
    }
    
    // 快速冪組合
    vector<int> A(m-1, 0); A[0] = 1; // 單位元
    ll exp = n;
    while (exp) {
        if (exp & 1) {
            vector<int> tmp;
            multiply(A, B, tmp);
            A.swap(tmp);
        }
        if (exp > 1) {
            vector<int> tmp;
            multiply(B, B, tmp);
            B.swap(tmp);
        }
        exp >>= 1;
    }
    
    cout << A[ind[target]] << '\n';
    return 0;
}