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