1. 程式人生 > >BZOJ3992: [SDOI2015]序列統計(NTT 原根 生成函數)

BZOJ3992: [SDOI2015]序列統計(NTT 原根 生成函數)

static sina har sts 原根 bre printf 個數 bzoj3

題意

題目鏈接

給出大小為\(S\)的集合,從中選出\(N\)個數,滿足他們的乘積\(\% M = X\)的方案數

Sol

神仙題Orz

首先不難列出最裸的dp方程,設\(f[i][j]\)表示選了\(i\)個數,他們的乘積為\(j\)的方案數

\(g[k] = [\exists a_i = k]\)

轉移的時候

\[f[i + 1][(j * k) \% M] += f[i][j] * g[k]\]

不難發現每次的轉移都是相同的,因此可以直接矩陣快速冪,時間復雜度變為\(logN M^2\)

觀察上面的式子,如果我們能把\((j * k) \% M\),變成\((j + k) \% M\)

的話,就是一個循環卷積的形式了

這裏可以用原根來實現,設\(g\)表示\(M\)的原根,\(mp[i] = j\)表示\(g^i = j\)

直接對每個物品構造生成函數,利用mp轉移即可

因為轉移是個循環卷積,所以統計答案的時候應該把第\(i\)項和第\(i+m-1\)項的系數加起來

至於為啥只統計一項。

技術分享圖片

#include<bits/stdc++.h>
using namespace std;
const int mod = 1004535809, G = 3, Gi = 334845270, MAXN = 1e5 + 10; 
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, X, S;
int r[MAXN], lim, L, ind[MAXN], s[MAXN], f[MAXN], a[MAXN], b[MAXN];
int mul(int a, int b) {
    return 1ll * a * b % mod;
}
int add(int x, int y) {
    if(x + y < 0) return x + y + mod;
    return x + y >= mod ? x + y - mod : x + y;
}
int dec(int x, int y) {
    return x - y < 0 ? x - y + mod : x - y;
}
int fp(int a, int p, int mod) {
    int base = 1;
    while(p) {
        if(p & 1) base = 1ll * base * a % mod; 
        a = 1ll * a * a % mod; p >>= 1;
    }
    return base;
}
int GetG(int x) {
    static int q[MAXN]; int tot = 0, tp = x - 1;
    for(int i = 2; i * i <= tp; i++) {
        if(!(tp % i)) {
            q[++tot] = i;
            while(!(tp % i)) tp /= i;
        }
    }
    if(tp > 1) q[++tot] = tp;
    for(int i = 2, j; i <= x - 1; i++) {
        for(j = 1; j <= tot; j++) if(fp(i, (x - 1) / q[j], x) == 1) break;
        if(j == tot + 1) return i;
    }
}
void NTT(int *a, int N, int type) {
    for(int i = 1; i < N; i++) if(i < r[i]) swap(a[i], a[r[i]]);
    for(int mid = 1; mid < N; mid <<= 1) {
        int R = mid << 1, Wn = fp(type == 1 ? G : Gi, (mod - 1) / R, mod);
        for(int j = 0; j < lim; j += R) {
            for(int w = 1, k = 0; k < mid; k++, w = mul(w, Wn)) {
                int x = a[j + k], y = mul(w, a[j + k + mid]);
                a[j + k] = add(x, y);
                a[j + k + mid] = dec(x, y);
            }
        }
    }
    if(type == -1) {
        for(int i = 0, inv = fp(lim, mod - 2, mod); i < N; i++) a[i] = mul(a[i], inv);
    }
}
void mul(int *a1, int *b1, int *c) {
    memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b));//tag
    for(int i = 0; i < M - 1; i++) a[i] = a1[i], b[i] = b1[i];
    NTT(a, lim, 1); NTT(b, lim, 1);
    for(int i = 0; i < lim; i++) a[i] = mul(a[i], b[i]);
    NTT(a, lim, -1);
    for(int i = 0; i < M - 1; i++) c[i] = add(a[i], a[i + M - 1]);
}
void Pre() {
    lim = 1;
    while(lim <= 2 * (M - 2)) lim <<= 1, L++;
    for(int i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | (i & 1) << (L - 1);
    int d = GetG(M);
    for(int i = 0; i < M - 1; i++) ind[fp(d, i, M)] = i;
}
int main() {
    N = read(); M = read(); X = read(); S = read();
    Pre();
    for(int i = 1; i <= S; i++) {
        int x = read();
        if(x) f[ind[x]]++;
    }
    s[ind[1]] = 1;
    while(N) {
        if(N & 1) mul(s, f, s);
        mul(f, f, f); N >>= 1;
    }
    printf("%d", s[ind[X]]);
    return 0;
}
/*
40000000 3 1 2
1 2

4 3 1 2
1 2
*/

BZOJ3992: [SDOI2015]序列統計(NTT 原根 生成函數)