经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置:技术经验 » 大数据/云/AI » 人工智能基础 » 查看文章
构造照亮世界——快速沃尔什变换 (FWT)
来源:cnblogs  作者:ZnPdCo  时间:2024/5/8 9:11:47  对本文有异议

博客园

我的博客

快速沃尔什变换解决的卷积问题

快速沃尔什变换(FWT)是解决这样一类卷积问题:

\[c_i=\sum_{i=j\odot k}a_jb_k \]

其中,\(\odot\) 是位运算的一种。举个例子,给定数列 \(a,b\),求:

\[c_i=\sum_{j\oplus k=i} a_jb_k \]

FWT 的思想

看到 FWT 的名字,我们可以联想到之前学过的 FFT(很可惜,我没有写过 FFT 的笔记,所以没有链接),先看看 FFT 的原理:

  1. \(a,b\) 变换为 \(A,B\)\(O(n\log n)\)
  2. 通过 \(C_i=A_iB_i\) 计算,\(O(n)\);
  3. \(C\) 变换回 \(c\)\(O(n\log n)\)

综上,时间复杂度是 \(O(n\log n)\) 的。

在 FFT 中,我们构造了 \(A,B\)\(a,b\) 的点值表示法,这么做满足 \(C_i=A_iB_i\) 且容易变换。

其实 FWT 的思想也是一样的,主要也是需要构造 \(A,B\),使得其满足 \(C_i=A_iB_i\) 且可以快速变换。下面我们举 \(\cup\)(按位或)、\(\cap\)(按位与)和 \(\oplus\)(按位异或)为例。

因为数列长度是 \(2\) 的幂会更好处理,所以下文认为数列长度为 \(2^n\)

按位或

\[c_i=\sum_{j\cup k=i} a_jb_k \]

我们可以构造 \(A_i=\sum_{i\cup j=i} a_i\)。看看为什么需要这么构造。

首先,它满足 \(C_i=A_iB_i\)

\[\begin{align} A_iB_i&=(\sum_{i\cup j=i} a_j)(\sum_{i\cup k=i} b_k) \&=\sum_{i\cup j=i}\sum_{i\cup k=i}a_jb_k \&=\sum_{i\cup j=i}\sum_{i\cup k=i}a_jb_k \&=\sum_{i\cup(j\cup k)=i}a_jb_k \&= C_i \end{align} \]

其次,它可以快速变换。举顺变换的例子。类比 FFT 的步骤,我们采用分治的方法来处理它。假设目前考虑到第 \(i\) 位,其中 \(A_0\)\(A_1\)\(i-1\) 位分治的结果:

\[A=\text{merge}(A_0, A_0+A_1) \]

其中,\(A_0\) 是数列 \(A\) 的左半部分,\(A_1\)\(A\) 的右半部分。\(\text{merge}\) 函数就是把两个数列像拼接字符串一样拼接起来。\(+\) 则是将两个数列对应相加。

这么做为什么是正确的呢?容易发现,\(A_0\) 恰好是当前处理到的二进制位为 \(0\) 的子数列,\(A_1\) 则是当前处理到的二进制位为 \(1\) 的子数列。若当前位为 \(0\),则只能取二进制位为 \(0\) 的子数列 \(A_0\) 才能使得 \(i\cup j=i\)。而若当前位为 \(1\),则两种序列都能取。


考虑逆变换,则是将加上的 \(A_0\) 减回去:

\[a=\text{merge}(a_0, a_1-a_0) \]

下面我们给出代码实现。容易发现顺变换和逆变换可以合并为一个函数,顺变换时 \(\text{type}=1\),逆变换时 \(\text{type}=-1\)

  1. void Or(ll *a, ll type) { // 迭代实现,常数更小
  2. for(ll x = 2; x <= n; x <<= 1) {
  3. ll k = x >> 1;
  4. for(ll i = 0; i < n; i += x) {
  5. for(ll j = 0; j < k; j ++) {
  6. (a[i + j + k] += a[i + j] * type) %= P;
  7. }
  8. }
  9. }
  10. }

按位与

\[c_i=\sum_{j\cap k=i} a_jb_k \]

同理构造 \(A_i=\sum_{i\cap j=i} a_i\)\(C_i=A_iB_i\) 的正确性不证了。

容易发现,\(A_0\) 恰好是当前处理到的二进制位为 \(0\) 的子数列,\(A_1\) 则是当前处理到的二进制位为 \(1\) 的子数列。若当前位为 \(1\),则只能取二进制位为 \(1\) 的子数列 \(A_0\) 才能使得 \(i\cap j=i\)。而若当前位为 \(0\),则两种序列都能取。

\[A=\text{merge}(A_0+A_1, A_1) \]

\[a=\text{merge}(a_0 - a_1, a_1) \]


下面我们给出代码实现。顺变换时 \(\text{type}=1\),逆变换时 \(\text{type}=-1\)

  1. void And(ll *a, ll type) {
  2. for(ll x = 2; x <= n; x <<= 1) {
  3. ll k = x >> 1;
  4. for(ll i = 0; i < n; i += x) {
  5. for(ll j = 0; j < k; j ++) {
  6. (a[i + j] += a[i + j + k] * type) %= P;
  7. }
  8. }
  9. }
  10. }

按位异或

发现异或有点难搞,但这怎么会难倒沃尔什大佬呢?我们引入一个新的运算符 \(\circ\)。定义 \(x\circ y=\text{popcnt}(x\cap y)\bmod 2\),其中 \(\text{popcnt}\) 表示二进制下 \(1\) 的个数,并重申一下 \(\cap\) 表示按位与。

不用慌,我们也不需要你真正实现一个 \(\text{popcnt}\),它仅仅只是作为一个理解的辅助罢了。

我们发现它满足 \((x\circ y)\oplus (x\circ z)=x\circ(y\oplus z)\)。(重申一下 \(\oplus\) 表示按位异或)

感性证明:发现这个新的运算符 \(\circ\) 其实就是 \(x\)\(y\) 相同位数的奇偶性。若 \((x\circ y)\oplus (x\circ z)=0\),则 \(x\)\(y\)\(x\)\(z\) 相同位数个数奇偶性相同,所以 \(y\oplus z\)\(x\) 相同位数个数奇偶性也是相同的 ;若 \((x\circ y)\oplus (x\circ z)=1\),则 \(x\)\(y\)\(x\)\(z\) 相同位数个数奇偶性不同,所以 \(y\oplus z\)\(x\) 相同位数个数奇偶性也是不同的。

\(A_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j\)。我们来证一下 \(C_i=A_iB_i\) 的正确性:

\[\begin{align} A_iB_i&=(\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j)(\sum_{i\circ k=0}b_k-\sum_{i\circ k=1}b_k) \&=(\sum_{i\circ j=0}a_j\sum_{i\circ k=0}b_k+\sum_{i\circ j=1}a_j\sum_{i\circ k=1}b_k)-(\sum_{i\circ j=0}a_j\sum_{i\circ k=1}b_k+\sum_{i\circ j=1}a_j\sum_{i\circ k=0}b_k) \&=\sum_{(j\oplus k)\circ i=0}a_jb_k-\sum_{(j\oplus k)\circ i=1}a_jb_k \&=C_i \end{align} \]

来看看怎么快速计算 \(A,B\) 的值,依旧是分治:

对于 \(i\) 在当前位为 \(0\) 的子数列 \(A_0\),进行 \(\circ\) 运算时发现它和 \(0\) 计算或和 \(1\) 计算结果都不会变(因为 \(0\cap 0=0,0\cap1=0\)),所以 \(A_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j\) 中的 \(\sum_{i\circ j=1}a_j=0\)

对于 \(i\) 在当前位为 \(1\) 的子数列 \(A_1\),进行 \(\circ\) 运算时发现它和 \(0\) 计算结果是 \(0\),和 \(1\) 计算结果是 \(1\)(因为 \(1\cap 0=0,1\cap1=1\))。

综上,有:

\[A=\text{merge}((A_0+A_1)-0, A_0-A_1) \]

也就是:

\[A=\text{merge}(A_0+A_1, A_0-A_1) \]

逆变换易得:

\[a=\text{merge}(\frac{a_0+a_1}{2}, \frac{a_0-a_1}{2}) \]

给出代码,顺变换时 \(\text{type}=1\),逆变换时 \(\text{type}=\frac{1}{2}\)

  1. void Xor(ll *a, ll type) {
  2. for(ll x = 2; x <= n; x <<= 1) {
  3. ll k = x >> 1;
  4. for(ll i = 0; i < n; i += x) {
  5. for(ll j = 0; j < k; j ++) {
  6. (a[i + j] += a[i + j + k]) %= P;
  7. (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
  8. (a[i + j] *= type) %= P;
  9. (a[i + j + k] *= type) %= P;
  10. }
  11. }
  12. }
  13. }

现在大家能去切前两道模板例题,并挑战一下后面的几道题目了。

从另一个角度看待 FWT

我们设 \(c(i,j)\)\(a_j\)\(A_i\) 的贡献系数。我们可以重新描述 FWT 变换的过程:

\[A_i = \sum_{j=0}^{n-1} c(i,j) a_j \]

因为有:

\[A_iB_i=C_i \]

所以我们可以通过简单的证明得到:\(c(i,j)c(i,k)=c(i,j\odot k)\)。其中 \(\odot\) 是任意一种位运算。

同时,\(c\) 函数还有一个重要的性质,它可以按位处理。

举个例子,我们变换的时候:

\[A_i = \sum_{j=0}^{n-1} c(i,j) a_j \]

这么做是比较劣的,我们将其拆分:

\[A_i = \sum_{j=0}^{(n-1)/2} c(i,j) a_j+\sum_{j=(n-1)/2+1}^{n-1} c(i,j) a_j \]

考虑前面的式子和后面的式子 \(i,j\) 的区别,发现只有最高位不同。

所以我们将 \(i,j\) 去除最高位的值为 \(i',j'\),并记 \(i_0\)\(i\) 的最高位。有:

\[A_i = c(i_0,0)\sum_{j=0}^{(n-1)/2} c(i',j') a_j+c(i_0,1)\sum_{j=(n-1)/2+1}^{n-1} c(i',j') a_j \]

如果 \(i_0=0\),则有:

\[A_i = c(0,0)\sum_{j=0}^{(n-1)/2} c(i',j') a_j+c(0,1)\sum_{j=(n-1)/2+1}^{n-1} c(i',j') a_j \]

\(i_0=1\) 则有:

\[A_i = c(1,0)\sum_{j=0}^{(n-1)/2} c(i',j') a_j+c(1,1)\sum_{j=(n-1)/2+1}^{n-1} c(i',j') a_j \]

也就是说,我们只需要:

\[\begin{bmatrix} c(0,0) & c(0,1) \c(1,0) & c(1,1) \end{bmatrix} \]

四个数就可以完成变换了。我们称这个矩阵为位矩阵。


如果我们要进行逆变换,则需要上面的位矩阵的逆矩阵。

若逆矩阵为 \(c^{-1}\),可以通过类似操作得到原数:

\[a_i = \sum_{j=0}^n c^{-1}(i,j) A_j \]

逆矩阵不一定存在,比如如果有一排 \(0\) 或者一列 \(0\) 那么这个矩阵就没有逆,我们在构造时需要格外小心。

按位或

我们可以构造:

\[\begin{bmatrix} 1 & 0 \1 & 1 \end{bmatrix} \]

这样满足 \(c(i,j)c(i,k)=c(i,j\cup k)\)。我们发现,这和我们前面推出的 \(A=\text{merge}(A_0, A_0+A_1)\) 一模一样!同理,下面也是一个满足这个条件的矩阵,但我们一般使用上面这个:

\[\begin{bmatrix} 1 & 1 \1 & 0 \end{bmatrix} \]

虽然下面这个矩阵也满足 \(c(i,j)c(i,k)=c(i,j\cup k)\),但这个矩阵存在一排 \(0\),不存在逆,所以不合法:

\[\begin{bmatrix} 0 & 0 \1 & 1 \end{bmatrix} \]

如果我们要进行逆变换,则需要对矩阵求逆,以最上面这个矩阵为例,得:

\[\begin{bmatrix} 1 & 0 \-1 & 1 \end{bmatrix} \]

然后按照顺变换的方法,把逆变换矩阵代入即可。

按位与

我们可以构造:

\[\begin{bmatrix} 1 & 1 \0 & 1 \end{bmatrix} \]

这样满足 \(c(i,j)c(i,k)=c(i,j\cap k)\)

逆矩阵:

\[\begin{bmatrix} 1 & -1 \0 & 1 \end{bmatrix} \]

按位异或

我们可以构造:

\[\begin{bmatrix} 1 & 1 \1 & -1 \end{bmatrix} \]

这样满足 \(c(i,j)c(i,k)=c(i,j\oplus k)\)

逆矩阵:

\[\begin{bmatrix} 0.5 & 0.5 \0.5 & -0.5 \end{bmatrix} \]

FWT 的性质

FWT 是线性变换。

\(FWT(X)\)\(X\) 的 FWT 变换,则有:

\[FWT(A+B)=FWT(A)+FWT(B) \]

以及:

\[FWT(cA)=cFWT(A) \]

这样就可以实现快速卷积,参考第四道例题。

K 维 FWT

max 运算

我们重新看看我们的 \(\cup\) 运算,发现他实际上就是二进制下的取 \(\max\)。我们将其拓展到 \(K\) 进制,有:

\[c(i,j)c(i,k)=c(i,\max(j,k)) \]

\(j=k\),那么上式又是:

\[c(i,j)c(i,j)=c(i,j) \]

也就是说,每一行的 \(1\) 必定只能在 \(0\) 的前面,如果在后面则不合法了。手玩一下可以发现一组合法构造:

\[\begin{bmatrix} 1 & 0 & 0 & 0 \1 & 1 & 0 & 0 \1 & 1 & 1 & 0 \1 & 1 & 1 & 1 \end{bmatrix} \]

求逆可得:

\[\begin{bmatrix} 1 & 0 & 0 & 0 \-1 & 1 & 0 & 0 \0 & -1 & 1 & 0 \0 & 0 & -1 & 1 \end{bmatrix} \]

min 运算

我们重新看看我们的 \(\cap\) 运算,发现他实际上就是二进制下的取 \(\min\)。我们将其拓展到 \(K\) 进制,有:

\[c(i,j)c(i,k)=c(i,\min(j,k)) \]

\(j=k\),那么上式又是:

\[c(i,j)c(i,j)=c(i,j) \]

也就是说,每一行的 \(1\) 必定只能在 \(0\) 的后面,如果在后面则不合法了。手玩一下可以发现一组合法构造:

\[\begin{bmatrix} 1 & 1 & 1 & 1 \0 & 1 & 1 & 1 \0 & 0 & 1 & 1 \0 & 0 & 0 & 1 \end{bmatrix} \]

求逆可得:

\[\begin{bmatrix} 1 & -1 & 0 & 0 \0 & 1 & -1 & 0 \0 & 0 & 1 & -1 \0 & 0 & 0 & 1 \end{bmatrix} \]

前两者用得比较少,用得比较多的是:

不进位加法

我们重新看看我们的 \(\oplus\) 运算,发现他实际上就是二进制下的不进位加法。我们将其拓展到 \(K\) 进制,有:

\[c(i,j)c(i,k)=c(i,(j+k)\bmod K) \]

我们构造 \(c(i,j)=\omega_{K}^j\),就可以满足要求了:

\[\omega_{K}^j\omega_{k}^k=\omega_{K}^{(j+k)\bmod K} \]

但是每一行都一样矩阵也没有逆,所以我们可以构造 \(c(i,j)=\omega_{K}^{(i-1)j}\) 即可。

有下面这个矩阵:

\[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \1 & \omega_{K}^1 & \omega_{K}^2 & \cdots & \omega_{K}^{k-1} \1 & \omega_{K}^2 & \omega_{K}^4 & \cdots & \omega_{K}^{2(k-1)} \1 & \omega_{K}^3 & \omega_{K}^6 & \cdots & \omega_{K}^{3(k-1)} \\vdots & \vdots & \vdots & \ddots & \vdots \1 & \omega_{K}^{k-1} & \omega_{K}^{2(k-1)} & \cdots & \omega_{K}^{k(k-1)} \end{bmatrix} \]

这不就是我们熟悉的范德蒙德矩阵吗?

现在我们也知道矩阵的逆了:

\[\frac{1}{K}\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \1 & \omega_{K}^{-1} & \omega_{K}^{-2} & \cdots & \omega_{K}^{-(k-1)} \1 & \omega_{K}^{-2} & \omega_{K}^{-4} & \cdots & \omega_{K}^{-2(k-1)} \1 & \omega_{K}^{-3} & \omega_{K}^{-6} & \cdots & \omega_{K}^{-3(k-1)} \\vdots & \vdots & \vdots & \ddots & \vdots \1 & \omega_{K}^{-(k-1)} & \omega_{K}^{-2(k-1)} & \cdots & \omega_{K}^{-k(k-1)} \end{bmatrix} \]

如果我们题目给出的模数是存在单位根的,我们就可以简单实现,可以参考第六道例题。


但是单位根在模意义下可能不存在,所以我们考虑扩域,就是人为地定义一个 \(x\),满足 \(x^K=1\),然后直接把 \(x\) 代入计算,这样每个数都是一个关于 \(x\)\(k-1\) 次多项式。我们只需要在 \(\pmod {x^K-1}\) 下计算即可。那么矩阵可以这么表示:

\[\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \1 & x^1 & x^2 & \cdots & x^{k-1} \1 & x^2 & x^4 & \cdots & x^{2(k-1)} \1 & x^3 & x^6 & \cdots & x^{3(k-1)} \\vdots & \vdots & \vdots & \ddots & \vdots \1 & x^{k-1} & x^{2(k-1)} & \cdots & x^{k(k-1)} \end{bmatrix} \]

但是这么做可能会存在零因子,也就是一个数有多种表示方法,我们无法确定一个数的真实值。

我们考虑不 \(\pmod {x^K-1}\) 了,我们 \(\bmod\) 分圆多项式 \(\Phi_{K}(x)\),他满足 \(x\) 的阶为 \(k\),且在 \(Q\) 上不可约。所以我们定义上面的计算是在 \(\pmod {\Phi_{K}(x)}\) 下进行即可。

另一方面,如何求分圆多项式,这一点可以在因式分解这道题的题解区里了解。这里给出分圆多项式的表:

还有一个问题是,\(\bmod \Phi_{K}(x)\) 常数大(因为 \(\Phi\) 本身就是一个多项式)。但是因为 \(\Phi_{K}(x)\mid x^k-1\),我们只需要在计算时 \(\bmod x^k -1\),最后再 \(\bmod \Phi_{K}(x)\) 即可。

具体实现参考第七道例题。

例题

「洛谷 P4717」 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)

\(\cup\)\(\cap\)\(\oplus\) 的三种卷积。

\(n\le17\)

这题也就是模板题了,下文直接给出代码:

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll long long
  4. #define P 998244353
  5. const ll N = 1 << 18;
  6. ll n;
  7. ll A[N], B[N];
  8. ll a[N], b[N];
  9. void init() {
  10. for(ll i = 0; i < n; i ++) a[i] = A[i], b[i] = B[i];
  11. }
  12. void Or(ll *a, ll type) {
  13. for(ll x = 2; x <= n; x <<= 1) {
  14. ll k = x >> 1;
  15. for(ll i = 0; i < n; i += x) {
  16. for(ll j = 0; j < k; j ++) {
  17. (a[i + j + k] += a[i + j] * type) %= P;
  18. }
  19. }
  20. }
  21. }
  22. void And(ll *a, ll type) {
  23. for(ll x = 2; x <= n; x <<= 1) {
  24. ll k = x >> 1;
  25. for(ll i = 0; i < n; i += x) {
  26. for(ll j = 0; j < k; j ++) {
  27. (a[i + j] += a[i + j + k] * type) %= P;
  28. }
  29. }
  30. }
  31. }
  32. void Xor(ll *a, ll type) {
  33. for(ll x = 2; x <= n; x <<= 1) {
  34. ll k = x >> 1;
  35. for(ll i = 0; i < n; i += x) {
  36. for(ll j = 0; j < k; j ++) {
  37. (a[i + j] += a[i + j + k]) %= P;
  38. (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
  39. (a[i + j] *= type) %= P;
  40. (a[i + j + k] *= type) %= P;
  41. }
  42. }
  43. }
  44. }
  45. void calc() {
  46. for(ll i = 0; i < n; i ++) (a[i] *= b[i]) %= P;
  47. }
  48. void print() {
  49. for(ll i = 0; i < n; i ++) printf("%lld ", (a[i] % P + P) % P);
  50. printf("\n");
  51. }
  52. int main() {
  53. scanf("%lld", &n);
  54. n = 1 << n;
  55. for(ll i = 0; i < n; i ++) scanf("%lld", &A[i]);
  56. for(ll i = 0; i < n; i ++) scanf("%lld", &B[i]);
  57. init(); Or(a, 1); Or(b, 1); calc(); Or(a, P - 1); print();
  58. init(); And(a, 1); And(b, 1); calc(); And(a, P - 1); print();
  59. init(); Xor(a, 1); Xor(b, 1); calc(); Xor(a, 499122177); print();
  60. }

「洛谷 P6097」 【模板】子集卷积

求:

\[c_k=\sum_{\substack{{i \cap j=0}\\{i\cup j=k}}} a_i b_j \]

\(n\le20\)

首先,下半部分是我们喜闻乐见的 FWT 常见形式,而上半部分我们可以看成是 \(i\)\(j\) 不交。有:

\[i\cup j=0\Rightarrow \text{popcnt}(i)+\text{popcnt}(j)=\text{popcnt}(i\cup j) \]

所以我们可以构造:

\[A_{i,j}=\sum_{\substack{{i\cup k=i}\\{\text{popcnt}(j)=k}}} a_i \]

可以枚举 \(\text{popcnt}\) 的值,分开考虑。

那么求 \(C\) 的时候有 \(C_{i,j}=\sum_{j=0}^n A_{i,k}B_{i,j-k}\)

然后就可以做了。

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define popcnt(x) __builtin_popcountll(x)
  4. #define ll long long
  5. const ll M = 20, N = 1 << M, P = 1e9 + 9;
  6. ll n, m;
  7. ll a[M + 1][N], b[M + 1][N], c[M + 1][N];
  8. void Or(ll *a, ll type) {
  9. for(ll x = 2; x <= n; x <<= 1) {
  10. ll k = x >> 1;
  11. for(ll i = 0; i < n; i += x) {
  12. for(ll j = 0; j < k; j ++) {
  13. (a[i + j + k] += a[i + j] * type) %= P;
  14. }
  15. }
  16. }
  17. }
  18. int main() {
  19. scanf("%lld" ,&m);
  20. n = 1 << m;
  21. for(ll i = 0; i < n; i ++) {
  22. scanf("%lld", &a[popcnt(i)][i]);
  23. }
  24. for(ll i = 0; i < n; i ++) {
  25. scanf("%lld", &b[popcnt(i)][i]);
  26. }
  27. for(ll i = 0; i <= m; i ++) {
  28. Or(a[i], 1);
  29. Or(b[i], 1);
  30. }
  31. for(ll i = 0; i <= m; i ++) {
  32. for(ll j = 0; j <= i; j ++) {
  33. for(ll k = 0; k < n; k ++) {
  34. (c[i][k] += a[j][k] * b[i - j][k]) %= P;
  35. }
  36. }
  37. }
  38. for(ll i = 0; i <= m; i ++) {
  39. Or(c[i], -1);
  40. }
  41. for(ll i = 0; i < n; i ++) {
  42. printf("%lld ", (c[popcnt(i)][i] % P + P) % P);
  43. }
  44. }

「牛客 881D」Parity of Tuples

给定 \(n\times m\) 的矩阵 \(a\),定义 \(\text{cnt}(x)\) 为矩阵中有多少行对于 \(x\) 是合法的,合法的定义为这一行中每一个数 \(a_{i,j}\cap x\) 的二进制值中都有奇数个 \(1\)

你需要求出对于所有的 \(x\)\(\text{cnt}\) 的取值。

\(n\le10^5,m\le10,x\le2^{20}\)

再次重申,\(\cap\) 是按位与的意思。

首先我们用数学公式定义一下 \(\text{cnt}\)(因为公式复杂,所以加了 \(\tt large\)):

\[\large \text{cnt}(x)=\frac{1}{2^m}\sum_{i=1}^n\prod_{j=1}^m (1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}) \]

说明一下正确性。如果 \(\text{popcnt}(a_{i,j}\cap x)\) 是奇数的话,那么 \((-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 的结果就是 \(-1\)。最后 \(1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 就是 \(2\),最后会被 \(\frac{1}{2^m}\) 除去;如果 \(\text{popcnt}(a_{i,j}\cap x)\) 是偶数的话,那么 \((-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 的结果就是 \(1\)。最后 \(1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}\) 就是 \(0\),那么整行的结果都是 \(0\)

然后我们发现它是可以展开的:

\[\large \begin{align} \prod_{j=1}^m (1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}) &= (1-(-1)^{\text{popcnt}(a_{i,1}\cap x)})(1-(-1)^{\text{popcnt}(a_{i,2}\cap x)})\cdots(1-(-1)^{\text{popcnt}(a_{i,m}\cap x)}) \&= 1 - \sum_{a=1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)} + \sum_{a=1}^m\sum_{b=a+1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)+\text{popcnt}(a_{i,b}\cap x)} - \& \sum_{a=1}^m\sum_{b=a+1}^m\sum_{c=b+1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)+\text{popcnt}(a_{i,b}\cap x)+\text{popcnt}(a_{i,c}\cap x)} + \cdots \end{align} \]

然后我们有一个性质:

\[\large (-1)^{\sum_{i=1}^n\text{popcnt}(a_i\cap x)}=(-1)^{\text{popcnt}((\oplus_{i=1}^na_i)\cap x)} \]

也就是 \(\sum_{i=1}^n\text{popcnt}(a_i\cap x)\) 的奇偶性与 \(\text{popcnt}((\oplus_{i=1}^na_i)\cap x)\) 的相同。这点在上面的新的运算符 \(\circ\) 的性质中有类似的体现。

容易得到:

\[\large\begin{align} &=1 - \sum_{a=1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)} + \sum_{a=1}^m\sum_{b=a+1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)+\text{popcnt}(a_{i,b}\cap x)} - \& \sum_{a=1}^m\sum_{b=a+1}^m\sum_{c=b+1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)+\text{popcnt}(a_{i,b}\cap x)+\text{popcnt}(a_{i,c}\cap x)} + \cdots \ &=1 - \sum_{a=1}^m (-1)^{\text{popcnt}(a_{i,a}\cap x)} + \sum_{a=1}^m\sum_{b=a+1}^m (-1)^{\text{popcnt}((a_{i,a}\oplus a_{i,b})\cap x)} - \& \sum_{a=1}^m\sum_{b=a+1}^m\sum_{c=b+1}^m (-1)^{\text{popcnt}((a_{i,a}\oplus a_{i,b}\oplus a_{i,c})\cap x)} + \cdots \\end{align} \]

我们发现一加一减的可以容斥,我们容斥计算 \(f_i\) 表示 \(n\) 行的所有式子中 \((-1)^i\) 前面的系数和。

  1. // num 处理到第几列
  2. // x 当前的指数
  3. // mu 当前的系数(+1 or -1)
  4. void dfs(ll *a, ll num, ll x, ll mu) {
  5. if(num > m) {
  6. f[x] += mu;
  7. return;
  8. }
  9. dfs(a, num + 1, x, mu); // 不加入第 num 列,系数不变
  10. dfs(a, num + 1, x ^ a[num], -mu);
  11. }

这样我们就可以进一步化简:

\[\begin{align} &= \sum_{i=0}^{2^k-1} f_{x\cap i}(-1)^{x\cap i} \end{align} \]

我们突然发现后面这个 \((-1)^i\) 取值只有两种,当 \(x\cap i\) 是奇数时取值为 \(-1\),否则为 \(1\)

好了,现在我们的问题转换为了求出:

\[\sum_{\text{popcnt}(x\cap i)\bmod2=0} f_i-\sum_{\text{popcnt}(x\cap i)\bmod2=1} f_i \]

这不就是 FWT 中的异或变换吗:

\[A_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j \]

综上,我们发现这题就是推式子容斥之后得到 FWT 的形式。


原题需要将输出加密:

\[\bigoplus\limits_{x = 0}^{2^k - 1} (\text{cnt}(x) \times 3^x \bmod (10^9 + 7)) \]

  1. #include <bits/stdc++.h>
  2. #define ll long long
  3. using namespace std;
  4. const ll P = 1e9 + 7;
  5. #define N 100010
  6. #define M 20
  7. #define K 21
  8. ll n, m, k;
  9. ll f[1 << K];
  10. ll a[N][M];
  11. // num 处理到第几列
  12. // x 当前的指数
  13. // mu 当前的系数(+1 or -1)
  14. void dfs(ll *a, ll num, ll x, ll mu) {
  15. if(num > m) {
  16. f[x] += mu;
  17. return;
  18. }
  19. dfs(a, num + 1, x, mu); // 不加入第 num 列,系数不变
  20. dfs(a, num + 1, x ^ a[num], -mu);
  21. }
  22. void Xor(ll *a, ll type) {
  23. for(ll x = 2; x <= (1 << k); x <<= 1) {
  24. ll z = x >> 1;
  25. for(ll i = 0; i < (1 << k); i += x) {
  26. for(ll j = 0; j < z; j ++) {
  27. (a[i + j] += a[i + j + z]) %= P;
  28. (a[i + j + z] = a[i + j] - 2 * a[i + j + z]) %= P;
  29. (a[i + j] *= type) %= P;
  30. (a[i + j + z] *= type) %= P;
  31. }
  32. }
  33. }
  34. }
  35. ll qpow(ll x, ll y) {
  36. if(y == 0) return 1;
  37. if(y % 2 == 1) return x * qpow(x, y - 1) % P;
  38. ll tmp = qpow(x, y / 2);
  39. return tmp * tmp % P;
  40. }
  41. int main() {
  42. while(scanf("%lld %lld %lld", &n, &m, &k) != EOF) {
  43. for(ll i = 0; i < (1 << k); i ++) f[i] = 0;
  44. for(ll i = 1; i <= n; i ++) {
  45. for(ll j = 1; j <= m; j ++) {
  46. scanf("%lld", &a[i][j]);
  47. }
  48. dfs(a[i], 1, 0, 1);
  49. }
  50. Xor(f, 1);
  51. ll pw = 1, inv = qpow(1 << m, P - 2), ans = 0;
  52. for(ll i = 0; i < (1 << k); i ++) {
  53. ans ^= f[i] * pw % P * inv % P;
  54. (pw *= 3) %= P;
  55. }
  56. printf("%lld\n", ans);
  57. }
  58. }

「AT ABC212H」 Nim Counting

给定两个数 \(N,K\),以及一个长度为 \(K\) 的整数数组 \((A_1,A_2,\cdots, A_K)\)

两个人玩 Nim 游戏。

现在通过以下方式生成一个游戏:

任意选择一个 \(1\le M\le N\)\(M\) 表示石子堆数。

对于每一堆,其石子数是 \(A\) 中任意一个数。

对于 \(\sum_{i=1}^N K^i\) 种游戏,求先手获胜的游戏数,答案对 \(998244353\) 取模。

\(n\le2\times10^5,K\le2^{16},a_i\le2^{16}\)

根据玩 Nim 游戏的经验,可以发现先手获胜当且仅当 \(\bigoplus_{i=0}^n A_i\neq 0\)

所以我们定义 dp 式子 \(f_{i,j}\) 表示有 \(i\) 个石堆,且石堆异或和为 \(j\) 的获胜方案数,有:

\[f_{i-1,j}\to \sum_{k=1}^Kf_{i,j\oplus a_k} \]

答案就是 \(\sum_{i=1}^n\sum_{j\neq0} f_{i,j}\)

直接转移是朴素的,发现上面的式子刚好是 FWT 异或操作,也就是:

\[f_{i,j}=\sum_{k\oplus x=j} f_{i-1,k}a_x \]

我们定义 \(a\) 是一个全是 \(1\) 的数组即可。

同时,我们发现其实不需要真的进行 \(n\) 次卷积,其实只需要将 FWT 变换过之后的结果 \(A\),求出 \(A+A^2+A^3+\cdots+A^n\) 即可。

上面的可以通过等比数列求和公式计算。

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll long long
  4. #define P 998244353
  5. const ll K = 1 << 20;
  6. ll n, k, ans;
  7. ll f[K];
  8. void FWT(ll *a, ll type) {
  9. for(ll x = 2; x <= K; x <<= 1) {
  10. ll k = x >> 1;
  11. for(ll i = 0; i < K; i += x) {
  12. for(ll j = 0; j < k; j ++) {
  13. (a[i + j] += a[i + j + k]) %= P;
  14. (a[i + j + k] = a[i + j] - 2 * a[i + j + k]) %= P;
  15. (a[i + j] *= type) %= P;
  16. (a[i + j + k] *= type) %= P;
  17. }
  18. }
  19. }
  20. }
  21. ll qpow(ll x, ll y) {
  22. if(y == 0) return 1;
  23. if(y % 2 == 1) return x * qpow(x, y - 1) % P;
  24. ll tmp = qpow(x, y / 2);
  25. return tmp * tmp % P;
  26. }
  27. int main() {
  28. scanf("%lld %lld", &n, &k);
  29. for(ll i = 1; i <= k; i ++) {
  30. ll x;
  31. scanf("%lld", &x);
  32. f[x] ++;
  33. }
  34. FWT(f, 1);
  35. for(ll i = 0; i < K; i ++) {
  36. if(f[i] == 1) f[i] = n;
  37. else {
  38. f[i] = f[i] * (qpow(f[i], n) - 1) % P * qpow(f[i] - 1, P - 2) % P;
  39. }
  40. }
  41. FWT(f, 499122177);
  42. for(ll i = 1; i < K; i ++) {
  43. (ans += f[i]) %= P;
  44. }
  45. printf("%lld", (ans % P + P) % P);
  46. }

「AT ARC100E」 Or Plus Max

给你一个长度为 \(2^n\) 的序列 \(a\),每个\(1\le K\le 2^n-1\),找出最大的 \(a_i+a_j\)\(i \cup j \le K\)\(0 \le i < j < 2^n\))并输出。

\(n\le18\)

就是要求 \(\max_{i\cup j=k} a_i+a_j\)

我们维护 \(f_{i}\) 表示 \(\max_{i\cup j=i} a_i\)\(g_i=\text{max2}_{i\cup j=i} a_i\)\(\text{max2}\) 表示次大值。

然后就像 FWT 的或变换一样了。

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll long long
  4. const ll N = 1 << 18;
  5. ll n;
  6. struct node {
  7. ll mx1, mx2;
  8. node(ll a = 0, ll b = 0):mx1(a), mx2(b) {}
  9. friend node operator +(const node &x, const node &y) {
  10. if(x.mx1 > y.mx1) {
  11. return node(x.mx1, max(x.mx2, y.mx1));
  12. }
  13. return node(y.mx1, max(y.mx2, x.mx1));
  14. }
  15. } a[N];
  16. void FWT(node *a) {
  17. for(ll x = 2; x <= n; x <<= 1) {
  18. ll k = x >> 1;
  19. for(ll i = 0; i < n; i += x) {
  20. for(ll j = 0; j < k; j ++) {
  21. a[i + j + k] = a[i + j] + a[i + j + k];
  22. }
  23. }
  24. }
  25. }
  26. int main() {
  27. scanf("%lld", &n);
  28. n = 1 << n;
  29. for(ll i = 0; i < n; i ++) {
  30. scanf("%lld", &a[i].mx1);
  31. }
  32. FWT(a);
  33. for(ll i = 0; i < n; i ++) {
  34. a[i].mx1 = a[i].mx1 + a[i].mx2;
  35. }
  36. ll ans = 0;
  37. for(ll i = 1; i < n; i ++) {
  38. ans = max(ans, a[i].mx1);
  39. printf("%lld\n", ans);
  40. }
  41. }

「HDU 6618」 Good Numbers

定义一个正整数 \(n\) 是好数当且仅当 \(n\) 在 8 进制表示下所有的数码出现的次数为 3 的倍数(出现 0 次亦可)。

有多少个 \(k\) 位的 8 进制数(不含前导 0),满足这个数是好的,且是 \(p\) 的倍数。对 \(10^9+9\) 取模。

例如:当 \(k=3,p=2\) 时,好数有 \(222,444,666\) 三个。

\(k\le10^{18},p<8\)

考虑状压 dp,设 \(f_{i,s,j}\) 表示第 \(i\) 位,\(8\) 种数出现次数对 \(3\) 取模的状压情况,以及数对 \(p\) 取模的结果为 \(j\)

答案就是 \(f_{k,0,0}\)

直接暴力枚举位数转移是朴素的,瓶颈在于 \(k\),考虑优化掉 \(k\)

发现我们可以使用像快速幂一样的方法,也就是倍增 dp。

转移公式就是:

\[f_{2i,s_1\oplus s_2,j_1+t\times j_2}\gets f_{i,s_1,j_1}f_{i,s_2,j_2} \]

其中 \(t\) 是转移的位数,而 \(\oplus\) 在这里是不进位三进制加法。

发现这样多了瓶颈——我们需要枚举 \(s_1\)\(s_2\)

但是我们发现这不就是 FWT 中异或的形式吗:\(c_{i\oplus j}\gets a_ib_j\)。考虑三进制 FWT 加速。下面给出 FWT 的代码,w1 是原根的一次方,w2 是原根的二次方:

  1. void FWT(ll *a, ll type) {
  2. for (ll x = 3; x <= N; x *= 3) {
  3. ll k = x / 3;
  4. for (ll i = 0; i < N; i += x) {
  5. for (ll j = 0; j < k; j ++) {
  6. for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k];
  7. if (type == 1) {
  8. tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
  9. tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
  10. tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
  11. } else {
  12. tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
  13. tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
  14. tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
  15. for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P;
  16. }
  17. for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l];
  18. }
  19. }
  20. }
  21. }

因为 \(1e9+9\) 存在原根 \(2\),然后就朴素实现了,注意位矩阵:

\[\begin{bmatrix} 1 & 1 & 1 \1 & \omega_{3}^1 & \omega_{3}^2 \1 & \omega_{3}^2 & \omega_{3}^4 \\end{bmatrix} \]

代码:

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll long long
  4. const ll P = 1e9 + 9;
  5. ll qpow(ll x, ll y) {
  6. if(y == 0) return 1;
  7. if(y % 2 == 1) return x * qpow(x, y - 1) % P;
  8. ll tmp = qpow(x, y / 2);
  9. return tmp * tmp % P;
  10. }
  11. const ll G = 2;
  12. const ll w1 = qpow(G, (P - 1) / 3);
  13. const ll w2 = qpow(G, (P - 1) / 3 * 2);
  14. const ll inv3 = qpow(3, P - 2);
  15. const ll N = 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3;
  16. ll n, p;
  17. ll tmp[8][N], res[8][N], one[8][N];
  18. ll a[8][N], b[8][N];
  19. ll pw3[8];
  20. ll tmp1[3], tmp2[3];
  21. void FWT(ll *a, ll type) {
  22. for (ll x = 3; x <= N; x *= 3) {
  23. ll k = x / 3;
  24. for (ll i = 0; i < N; i += x) {
  25. for (ll j = 0; j < k; j ++) {
  26. for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k];
  27. if (type == 1) {
  28. tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
  29. tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
  30. tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
  31. } else {
  32. tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P;
  33. tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P;
  34. tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P;
  35. for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P;
  36. }
  37. for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l];
  38. }
  39. }
  40. }
  41. }
  42. ll base = 1;
  43. void fun(ll x) {
  44. if(x == 1) {
  45. memset(res, 0, sizeof res);
  46. memset(tmp, 0, sizeof tmp);
  47. memset(one, 0, sizeof one);
  48. for(ll i = 1; i < 8; i ++) res[i % p][pw3[i]] = 1;
  49. for(ll i = 0; i < 8; i ++) tmp[i % p][pw3[i]] = 1;
  50. for(ll i = 0; i < 8; i ++) one[i % p][pw3[i]] = 1;
  51. for (int i = 0; i < p; i ++) {
  52. FWT(tmp[i], 1);
  53. FWT(res[i], 1);
  54. FWT(one[i], 1);
  55. }
  56. base = 8 % p;
  57. return;
  58. }
  59. if(x % 2 == 1) {
  60. fun(x - 1);
  61. memset(a, 0, sizeof a);
  62. memset(b, 0, sizeof b);
  63. for (ll i = 0; i < p; i ++) {
  64. for (ll j = 0; j < p; j ++) {
  65. ll k = (i * 8 + j) % p;
  66. for (ll x = 0; x < N; x ++)
  67. (a[k][x] += tmp[i][x] * one[j][x]) %= P,
  68. (b[k][x] += res[i][x] * one[j][x]) %= P;
  69. }
  70. }
  71. memcpy(tmp, a, sizeof a);
  72. memcpy(res, b, sizeof b);
  73. (base *= 8) %= P;
  74. return;
  75. }
  76. fun(x / 2);
  77. memset(a, 0, sizeof a);
  78. memset(b, 0, sizeof b);
  79. for (ll i = 0; i < p; i ++) {
  80. for (ll j = 0; j < p; j ++) {
  81. ll k = (i * base + j) % p;
  82. for (ll x = 0; x < N; x ++)
  83. (a[k][x] += tmp[i][x] * tmp[j][x]) %= P,
  84. (b[k][x] += res[i][x] * tmp[j][x]) %= P;
  85. }
  86. }
  87. memcpy(tmp, a, sizeof a);
  88. memcpy(res, b, sizeof b);
  89. (base *= base) %= p;
  90. }
  91. int main() {
  92. pw3[0] = 1;
  93. for(ll i = 1; i <= 8; i ++) {
  94. pw3[i] = pw3[i - 1] * 3;
  95. }
  96. while(scanf("%lld %lld", &n, &p) != EOF) {
  97. fun(n);
  98. FWT(res[0], -1);
  99. printf("%lld\n", res[0][0]);
  100. }
  101. }

「CF 1103E」Radix sum

给定一个长度为 \(n\) 的序列 \(a_1,a_2,...,a_n\),对于每一个 \(p \in [0,n-1]\),求满足下列条件的整数序列 \(i_1,i_2,...,i_n\) 的方案数,对 \(2^{58}\) 取模:

  • \(\forall j \in [1,n] , i_j \in [1,n]\)
  • \(\sum\limits_{j=1}^n a_{i_j} = p\),这里的加法定义为十进制不进位加法。

\(n\le10^5,a_i\le10^5\)

我们可以想到 dp:设计状态 \(f_{i,s}\) 表示考虑到第 \(i\) 个数,当前加法状态为 \(s\)。因为 FWT 变换时线性的,可以先变换为 FWT 点值表示法,然后变成自己的 \(n\) 次幂,最后再变换回来。

上面是平凡的,但是题目给出了模数 \(2^{58}\)。发现没有单位根,所以考虑扩域。

这里的分圆多项式 \(\Phi_{10}(x)=x^4-x^3+x^2-x+1\)

然而我们发现 IFWT 时,需要除去进制 \(10\),然而我们发现 \(10\)\(2^{58}\) 下没有逆元。实际上我们发现 \(5\)\(2^{58}\) 下是有逆元的:\(57646075230342349\),我们只需要再除去一个 \(2\) 就可以了。设已经除以了 \(5\) 的答案为 \(x\),真正的答案为 \(y\),也就是 \(2^5y\equiv x\pmod{2^{64}}\),显然,我们有 \(y\equiv \frac{x}{2^5}\pmod{2^{64-5}}\),也就是 \(y\equiv \frac{x}{2^5}\pmod{2^{59}}\),所以直接将最后的答案除以 \(2^5\) 即可。虽然出题人不知道为什么要模 \(2^{58}\),但再取下模即可。

然后就是平凡实现了:

  1. #include <bits/stdc++.h>
  2. using namespace std;
  3. #define ll unsigned long long
  4. const ll P = 1ull << 58, N = 1e5 + 10;
  5. const ll m = 5, K = 10;
  6. ll inv5;
  7. ll n;
  8. ll pw[m + 1];
  9. ll qpow(ll x, ll y) {
  10. if(y == 0) return 1;
  11. if(y % 2 == 1) return x * qpow(x, y - 1);
  12. ll tmp = qpow(x, y / 2);
  13. return tmp * tmp;
  14. }
  15. struct poly {
  16. ll a[30];
  17. poly() {memset(a, 0, sizeof a);}
  18. ll operator [](ll x) const {return a[x];}
  19. ll& operator [](ll x) {return a[x];}
  20. friend poly operator *(const poly &x, const poly &y) {
  21. poly z;
  22. for(ll i = 0; i < K; i ++) {
  23. for(ll j = 0; j < K; j ++) {
  24. z[(i + j) % K] += x[i] * y[j];
  25. }
  26. }
  27. return z;
  28. }
  29. friend poly operator *(const poly &x, const ll &y) {
  30. poly z;
  31. for(ll i = 0; i < K; i ++) {
  32. z[i] += x[i] * y;
  33. }
  34. return z;
  35. }
  36. friend poly operator +(const poly &x, const poly &y) {
  37. poly z;
  38. for(ll i = 0; i < K; i ++) {
  39. z[i] += x[i] + y[i];
  40. }
  41. return z;
  42. }
  43. poly w(ll x) {
  44. poly res;
  45. for(ll i = 0; i < K; i ++) {
  46. res[(i + x) % K] += a[i];
  47. }
  48. return res;
  49. }
  50. } T, f[N], one;
  51. poly qpow(poly x, ll y) {
  52. if(y == 0) return one;
  53. if(y % 2 == 1) return x * qpow(x, y - 1);
  54. poly tmp = qpow(x, y / 2);
  55. return tmp * tmp;
  56. }
  57. poly tmp1[30], tmp2[30];
  58. void FWT(poly *a, ll type) {
  59. for(ll x = K; x <= pw[m]; x *= K) {
  60. ll k = x / K;
  61. for(ll i = 0; i < pw[m]; i += x) {
  62. for(ll j = 0; j < k; j ++) {
  63. for(ll l = 0; l < K; l ++) tmp1[l] = a[i + j + l * k], tmp2[l] = poly();
  64. if(type == 1) {
  65. for(ll l = 0; l < K; l ++) {
  66. for(ll v = 0; v < K; v ++) {
  67. tmp2[l] = tmp2[l] + tmp1[v].w(l * v % K);
  68. }
  69. }
  70. for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l];
  71. } else {
  72. for(ll l = 0; l < K; l ++) {
  73. for(ll v = 0; v < K; v ++) {
  74. tmp2[l] = tmp2[l] + tmp1[v].w((K - (l * v % K)) % K);
  75. }
  76. }
  77. for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l] * inv5;
  78. }
  79. }
  80. }
  81. }
  82. }
  83. ll mod(poly x){
  84. ll n = 4;
  85. for(ll i = K - 1; i >= n; i --){
  86. ll u = x[i];
  87. for(ll j = 1; j <= n; j ++) x[i - j] -= u * T[n - j];
  88. }
  89. ll u = x[0];
  90. u >>= m;
  91. return u % P;
  92. }
  93. int main() {
  94. pw[0] = 1;
  95. for(ll i = 1; i <= m; i ++) pw[i] = pw[i - 1] * K;
  96. T[0] = 1, T[1] = -1, T[2] = 1, T[3] = -1, T[4] = 1; // 分圆多项式phi10
  97. one[0] = 1;
  98. inv5 = 57646075230342349ull;
  99. scanf("%llu", &n);
  100. for(ll i = 1; i <= n; i ++) {
  101. ll x;
  102. scanf("%llu", &x);
  103. f[x][0] ++;
  104. }
  105. FWT(f, 1);
  106. for(ll i = 0; i < pw[m]; i ++) f[i] = qpow(f[i], n);
  107. FWT(f, -1);
  108. for(ll i = 0; i < n; i ++) cout<<mod(f[i])<<'\n';
  109. }

原文链接:https://www.cnblogs.com/znpdco/p/18172429

 友情链接:直通硅谷  点职佳  北美留学生论坛

本站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号