[BZOJ3992][SDOI2015]序列统计（DP+原根+NTT）

考虑一个DP方案：
$f[i][j]$表示到序列的第$i$个数，乘积模$M$等于$j$的序列个数。
然而，$O(M^2N)$的朴素DP和$O(M^3\log N)$的矩阵快速幂都过不去。
考虑到如果不是乘积模$M$等于$j$，而是和模$M$等于$j$，
那么可以利用快速幂+卷积实现$O(M\log M\log N)$的复杂度。
也就是说，如果把一个状态看作一个多项式，即
$f[i][1]+f[i][2]x+f[i][3]x^2+...+f[i][M-1]x^{M-2}$
把转移也看作一个多项式，即如果给定的集合内包含$x$，则$in[x]=1$，否则等于$0$，则转移对应的多项式为：
$in[1]+in[2]x+in[3]x^2+...+in[M-1]x^{M-2}$
则此时可以得出：$f[i+1]=f[i]*in$。（每一次多项式乘法之后，所有次数为$k$（$k>M-2$）的项都要加到次数为$k-(M-1)$的项上面去，并把次数为$k$的项清零）
最后结果就是$in^N[x]$。可以用FFT或NTT求出。
回到问题，可以发现乘积在一定条件下可以转换成和：由于此题中的$M$是质数，所以只需要求出$M$的原根$g$，并用$g^i$形式的式子代替问题中所有$[1,M)$的数，这样就能把乘积转换成和了。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int N = 6e5 + 5, ZZQ = 1004535809;
int n, PYZ, X, m, cnt[N], G, tot, cyx[N], maxn, orz, rev[N], ans[N],
to[N], res[N], A[N], B[N];
inline int qpow(int a, int b, const int &MX) {
    int res = 1;
    while (b) {
        if (b & 1) res = 1ll * res * a % MX;
        a = 1ll * a * a % MX;
        b >>= 1;
    }
    return res;
}
inline int orzPyzDalao() {
    int i, j, S = sqrt(PYZ - 1), tmp = PYZ - 1;
    for (i = 2; i <= S; i++) if (!(tmp % i)) {
        cyx[++tot] = i;
        while (!(tmp % i)) tmp /= i;
    }
    if (tmp > 1) cyx[++tot] = tmp;
    for (i = 2; i < PYZ; i++) {
        bool flag = 1;
        for (j = 1; j <= tot; j++)
            if (qpow(i, (PYZ - 1) / cyx[j], PYZ) == 1)
                {flag = 0; break;}
        if (flag) return i;
    }
    return -19370707;
}
inline void NTT(const int &n, int *a, const int &op) {
    int i, j, k, r = op == 1 ? 3 : 334845270;
    for (i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (k = 1; k < n; k <<= 1) {
        int x = qpow(r, (ZZQ - 1) / (k << 1), ZZQ);
        for (i = 0; i < n; i += (k << 1)) {
            int w = 1;
            for (j = 0; j < k; j++) {
                int u = a[i + j], v = 1ll * w * a[i + j + k] % ZZQ;
                a[i + j] = (u + v) % ZZQ; a[i + j + k] = (u - v + ZZQ) % ZZQ;
                w = 1ll * w * x % ZZQ;
            }
        }
    }
}
inline void prod(const int &n, int *a, int *b) {
    int i; NTT(n, a, 1); NTT(n, b, 1);
    for (i = 0; i <= n; i++) ans[i] = 1ll * a[i] * b[i] % ZZQ, a[i] = b[i] = 0;
    NTT(n, ans, -1); for (i = 0; i <= n; i++)
        ans[i] = 1ll * ans[i] * qpow(n, ZZQ - 2, ZZQ) % ZZQ;
    for (i = PYZ - 1; i <= n; i++)
        (ans[i % (PYZ - 1)] += ans[i]) %= ZZQ, ans[i] = 0;
}
int main() {
    int i; n = read(); PYZ = read(); X = read(); m = read();
    G = orzPyzDalao(); orz = 14; maxn = 1 << 14;
    for (i = 0; i <= maxn; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << orz - 1);
    int x = 1; for (i = 0; i < PYZ - 1; i++)
        to[x] = i, x = 1ll * x * G % PYZ;
    while (m--) {
        x = read(); if (x) cnt[to[x % PYZ]] = 1;
    }
    res[0] = 1; while (n) {
        if (n & 1) {
            for (i = 0; i <= maxn; i++) A[i] = res[i], B[i] = cnt[i];
            prod(maxn, A, B);
            for (i = 0; i <= maxn; i++) res[i] = ans[i];
        }
        for (i = 0; i <= maxn; i++) A[i] = B[i] = cnt[i];
        prod(maxn, A, B);
        for (i = 0; i <= maxn; i++) cnt[i] = ans[i];
        n >>= 1;
    }
    cout << res[to[X]] << endl;
    return 0;
}