经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置:技术经验 » 程序设计 » C++ » 查看文章
BZOJ3992: [SDOI2015]序列统计(NTT 原根 生成函数)
来源:cnblogs  作者:自为风月马前卒  时间:2018/11/28 9:54:26  对本文有异议

题意

题目链接

给出大小为\(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\)项的系数加起来

至于为啥只统计一项。

  1. #include<bits/stdc++.h>
  2. using namespace std;
  3. const int mod = 1004535809, G = 3, Gi = 334845270, MAXN = 1e5 + 10;
  4. inline int read() {
  5. char c = getchar(); int x = 0, f = 1;
  6. while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
  7. while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
  8. return x * f;
  9. }
  10. int N, M, X, S;
  11. int r[MAXN], lim, L, ind[MAXN], s[MAXN], f[MAXN], a[MAXN], b[MAXN];
  12. int mul(int a, int b) {
  13. return 1ll * a * b % mod;
  14. }
  15. int add(int x, int y) {
  16. if(x + y < 0) return x + y + mod;
  17. return x + y >= mod ? x + y - mod : x + y;
  18. }
  19. int dec(int x, int y) {
  20. return x - y < 0 ? x - y + mod : x - y;
  21. }
  22. int fp(int a, int p, int mod) {
  23. int base = 1;
  24. while(p) {
  25. if(p & 1) base = 1ll * base * a % mod;
  26. a = 1ll * a * a % mod; p >>= 1;
  27. }
  28. return base;
  29. }
  30. int GetG(int x) {
  31. static int q[MAXN]; int tot = 0, tp = x - 1;
  32. for(int i = 2; i * i <= tp; i++) {
  33. if(!(tp % i)) {
  34. q[++tot] = i;
  35. while(!(tp % i)) tp /= i;
  36. }
  37. }
  38. if(tp > 1) q[++tot] = tp;
  39. for(int i = 2, j; i <= x - 1; i++) {
  40. for(j = 1; j <= tot; j++) if(fp(i, (x - 1) / q[j], x) == 1) break;
  41. if(j == tot + 1) return i;
  42. }
  43. }
  44. void NTT(int *a, int N, int type) {
  45. for(int i = 1; i < N; i++) if(i < r[i]) swap(a[i], a[r[i]]);
  46. for(int mid = 1; mid < N; mid <<= 1) {
  47. int R = mid << 1, Wn = fp(type == 1 ? G : Gi, (mod - 1) / R, mod);
  48. for(int j = 0; j < lim; j += R) {
  49. for(int w = 1, k = 0; k < mid; k++, w = mul(w, Wn)) {
  50. int x = a[j + k], y = mul(w, a[j + k + mid]);
  51. a[j + k] = add(x, y);
  52. a[j + k + mid] = dec(x, y);
  53. }
  54. }
  55. }
  56. if(type == -1) {
  57. for(int i = 0, inv = fp(lim, mod - 2, mod); i < N; i++) a[i] = mul(a[i], inv);
  58. }
  59. }
  60. void mul(int *a1, int *b1, int *c) {
  61. memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b));//tag
  62. for(int i = 0; i < M - 1; i++) a[i] = a1[i], b[i] = b1[i];
  63. NTT(a, lim, 1); NTT(b, lim, 1);
  64. for(int i = 0; i < lim; i++) a[i] = mul(a[i], b[i]);
  65. NTT(a, lim, -1);
  66. for(int i = 0; i < M - 1; i++) c[i] = add(a[i], a[i + M - 1]);
  67. }
  68. void Pre() {
  69. lim = 1;
  70. while(lim <= 2 * (M - 2)) lim <<= 1, L++;
  71. for(int i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | (i & 1) << (L - 1);
  72. int d = GetG(M);
  73. for(int i = 0; i < M - 1; i++) ind[fp(d, i, M)] = i;
  74. }
  75. int main() {
  76. N = read(); M = read(); X = read(); S = read();
  77. Pre();
  78. for(int i = 1; i <= S; i++) {
  79. int x = read();
  80. if(x) f[ind[x]]++;
  81. }
  82. s[ind[1]] = 1;
  83. while(N) {
  84. if(N & 1) mul(s, f, s);
  85. mul(f, f, f); N >>= 1;
  86. }
  87. printf("%d", s[ind[X]]);
  88. return 0;
  89. }
  90. /*
  91. 40000000 3 1 2
  92. 1 2
  93. 4 3 1 2
  94. 1 2
  95. */
 友情链接:直通硅谷  点职佳  北美留学生论坛

本站QQ群:前端 618073944 | Java 606181507 | Python 626812652 | C/C++ 612253063 | 微信 634508462 | 苹果 692586424 | C#/.net 182808419 | PHP 305140648 | 运维 608723728

W3xue 的所有内容仅供测试,对任何法律问题及风险不承担任何责任。通过使用本站内容随之而来的风险与本站无关。
关于我们  |  意见建议  |  捐助我们  |  报错有奖  |  广告合作、友情链接(目前9元/月)请联系QQ:27243702 沸活量
皖ICP备17017327号-2 皖公网安备34020702000426号