题目地址：洛谷：【P3321】[SDOI2015]序列统计 – 洛谷、BZOJ：Problem 3992. — [SDOI2015]序列统计

题目描述

小C有一个集合S，里面的元素都是小于M的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为N的数列，数列中的每个数都属于集合S。小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数x，求所有可以生成出的，且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个。小C认为，两个数列{Ai}和{Bi}不同，当且仅当至少存在一个整数i，满足Ai≠Bi。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案mod 1004535809的值就可以了。

输入输出格式

输入格式：
一行，四个整数，N、M、x、|S|，其中|S|为集合S中元素个数。第二行，|S|个整数，表示集合S中的所有元素。

输出格式：
一行，一个整数，表示你求出的种类数mod 1004535809的值。

输入输出样例

输入样例#1：

4 3 1 2
1 2

输出样例#1：

8

题解

参考资料：bzoj 3992: [SDOI2015]序列统计 （NTT+快速幂+DP） – CSDN博客
本题需要的数学姿势：数学笔记：数论（欧拉函数、阶、原根） | KSkun’s Blog
首先，计数自然是一个动态规划，我们考虑设计状态dp[i][j]为选了i个数且结果mod M的值是j，转移如下
dp[i][j] = \sum_{k \in S} dp[i-1][j'] \ (j'k \bmod M = j)
直观上感觉直接做是一个O(n|S|)的，好像只能得10分。
接下来我们要考虑一个问题，模数很特殊，并不是常见的1000000007，而是1004535809，它的原根为3，是常用的NTT模数，难道出题人在指引我们使用NTT？
想用NTT做，首先得搞出卷积，卷积下标是加减的事情没办法像上面那样乘吧。那么考虑利用原根把乘法转换为加法，dp[i][j]表示选到第i个数结果mod M的值为g^j，改变以后的转移方程就成这样了
dp[i][j] = \sum_{g^k \bmod M \in S} dp[i-1][j'] \ ((j' + k) \bmod (M - 1) = j)
或者这样？用f[i]表示g^i是否在集合中。
dp[i][j] = \sum_{k=1}^{M-1} dp[i-1][j-k] \cdot f[k]
咦，后面那部分现在可以卷积了，而且你会发现每一层的转移中f数组都是不变的，也就是说这是个线性变换，参照矩阵快速幂的操作，我们可以做多项式的快速幂，即dp[n] = dp[0] \times f^n。答案便是 dp[n][i] \ (g^i=x) 。
于是我们的思路就很明确了：求原根、算出f数组、多项式快速幂。
没完呢！这里有一个坑点，集合内的元素可以有0，0显然对答案没贡献，所以你要跳过它处理。

代码

// Code by KSkun, 2018/4
#include <cstdio>
#include <cmath>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

inline LL fpow(LL n, LL k, LL p) {
    n %= p;
    LL res = 1;
    while(k) {
        if(k & 1) res = res * n % p;
        n = n * n % p;
        k >>= 1;
    }
    return res;
}


const int MAXN = 1 << 15, G = 3, MO = 1004535809;

int n, m, x, s, t;
LL f[MAXN], g[MAXN];
bool vis[MAXN];

LL rev[MAXN], dft1[MAXN], dft2[MAXN], dftr[MAXN];

inline void calrev(int n, int len) {
    for(int i = 0; i < n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));
    }
}

inline void dft(int n, LL *arr, int f) {
    for(int i = 0; i < n; i++) {
        if(i < rev[i]) std::swap(arr[i], arr[rev[i]]);
    }
    for(int i = 1; i < n; i <<= 1) {
        LL gn = fpow(G, (MO - 1) / (i << 1), MO);
        if(f == -1) gn = fpow(gn, MO - 2, MO);
        for(int j = 0; j < n; j += i << 1) {
            LL g = 1;
            for(int k = 0; k < i; k++) {
                LL x = arr[j + k], y = g * arr[j + k + i] % MO;
                arr[j + k] = (x + y) % MO;
                arr[j + k + i] = ((x - y) % MO + MO) % MO;
                g = g * gn % MO;
            }
        }
    }
    if(f == -1) {
        LL invn = fpow(n, MO - 2, MO);
        for(int i = 0; i < n; i++) {
            arr[i] = arr[i] * invn % MO;
        }
    }
}

inline void ntt(LL *a, LL *b) {
    memset(dft1, 0, sizeof(dft1));
    memset(dft2, 0, sizeof(dft2));
    memset(dftr, 0, sizeof(dftr));
    int n, len = 0;
    for(n = 1; n <= (m - 1) << 1; n <<= 1) len++;
    for(int i = 0; i < n; i++) {
        dft1[i] = a[i];
        dft2[i] = b[i];
    }
    calrev(n, len);
    dft(n, dft1, 1); dft(n, dft2, 1);
    for(int i = 0; i < n; i++) {
        dftr[i] = dft1[i] * dft2[i] % MO;
    }
    dft(n, dftr, -1);
    for(int i = 0; i < m - 1; i++) {
        a[i] = (dftr[i] + dftr[i + m - 1]) % MO;
    }
}

inline LL calg(LL n) {
    int fact[8005], tot = 0, x = n - 1;
    for(int i = 2; i * i <= x; i++) {
        if(x % i == 0) {
            fact[tot++] = (n - 1) / i;
            while(x % i == 0) x /= i;
        }
    }
    if(x > 1) fact[tot++] = (n - 1) / x;
    bool success;
    for(int i = 2;; i++) {
        success = true;
        for(int j = 0; j < tot; j++) {
            if(fpow(i, fact[j], n) == 1) {
                success = false;
                break;
            }
        }
        if(success) return i;
    }
}

int main() {
    n = readint(); m = readint(); x = readint(); s = readint();
    for(int i = 1; i <= s; i++) {
        t = readint();
        vis[t] = 1;
    }
    int xpow, G = calg(m), gn = 1;
    for(int i = 0; i < m - 1; i++) {
        if(vis[gn]) f[i] = 1;
        if(gn == x) xpow = i;
        gn = gn * G % m;
    }
    g[0] = 1;
    for(int i = n; i; i >>= 1) {
        if(i & 1) ntt(g, f);
        ntt(f, f);
    }
    printf("%lld", g[xpow]);
    return 0;
}