构造照亮世界——快速沃尔什变换 (FWT)

博客园

我的博客

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

快速沃尔什变换（FWT）是解决这样一类卷积问题：

ci=∑i=j⊙kajbkc_i=\sum_{i=j\odot k}a_jb_k
ci=i=j⊙k∑ajbk其中，⊙\odot⊙ 是位运算的一种。举个例子，给定数列 a,ba,ba,b，求：

ci=∑j⊕k=iajbkc_i=\sum_{j\oplus k=i} a_jb_k
ci=j⊕k=i∑ajbk## FWT 的思想

看到 FWT 的名字，我们可以联想到之前学过的 FFT（很可惜，我没有写过 FFT 的笔记，所以没有链接），先看看 FFT 的原理：

把 a,ba,ba,b 变换为 A,BA,BA,B，O(nlog⁡n)O(n\log n)O(nlogn)；
通过 Ci=AiBiC_i=A_iB_iCi=AiBi 计算，O(n)O(n)O(n);
把 CCC 变换回 ccc，O(nlog⁡n)O(n\log n)O(nlogn)。

综上，时间复杂度是 O(nlog⁡n)O(n\log n)O(nlogn) 的。

在 FFT 中，我们构造了 A,BA,BA,B 为 a,ba,ba,b 的点值表示法，这么做满足 Ci=AiBiC_i=A_iB_iCi=AiBi 且容易变换。

其实 FWT 的思想也是一样的，主要也是需要构造 A,BA,BA,B，使得其满足 Ci=AiBiC_i=A_iB_iCi=AiBi 且可以快速变换。下面我们举 ∪\cup∪（按位或）、∩\cap∩（按位与）和 ⊕\oplus⊕（按位异或）为例。

因为数列长度是 222 的幂会更好处理，所以下文认为数列长度为 2n2^n2n。

按位或

ci=∑j∪k=iajbkc_i=\sum_{j\cup k=i} a_jb_k
ci=j∪k=i∑ajbk我们可以构造 Ai=∑i∪j=iaiA_i=\sum_{i\cup j=i} a_iAi=∑i∪j=iai。看看为什么需要这么构造。

首先，它满足 Ci=AiBiC_i=A_iB_iCi=AiBi：

AiBi=(∑i∪j=iaj)(∑i∪k=ibk)=∑i∪j=i∑i∪k=iajbk=∑i∪j=i∑i∪k=iajbk=∑i∪(j∪k)=iajbk=Ci\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}
AiBi=(i∪j=i∑aj)(i∪k=i∑bk)=i∪j=i∑i∪k=i∑ajbk=i∪j=i∑i∪k=i∑ajbk=i∪(j∪k)=i∑ajbk=Ci其次，它可以快速变换。举顺变换的例子。类比 FFT 的步骤，我们采用分治的方法来处理它。假设目前考虑到第 iii 位，其中 A0A_0A0 和 A1A_1A1 是 i−1i-1i−1 位分治的结果：

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

这么做为什么是正确的呢？容易发现，A0A_0A0 恰好是当前处理到的二进制位为 000 的子数列，A1A_1A1 则是当前处理到的二进制位为 111 的子数列。若当前位为 000，则只能取二进制位为 000 的子数列 A0A_0A0 才能使得 i∪j=ii\cup j=ii∪j=i。而若当前位为 111，则两种序列都能取。

考虑逆变换，则是将加上的 A0A_0A0 减回去：

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

|  | void Or(ll *a, ll type) { // 迭代实现，常数更小 |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j + k] += a[i + j] * type) %= P;  |
|  |  } |
|  |  } |
|  |  } |
|  | } |

按位与

ci=∑j∩k=iajbkc_i=\sum_{j\cap k=i} a_jb_k
ci=j∩k=i∑ajbk同理构造 Ai=∑i∩j=iaiA_i=\sum_{i\cap j=i} a_iAi=∑i∩j=iai。Ci=AiBiC_i=A_iB_iCi=AiBi 的正确性不证了。

容易发现，A0A_0A0 恰好是当前处理到的二进制位为 000 的子数列，A1A_1A1 则是当前处理到的二进制位为 111 的子数列。若当前位为 111，则只能取二进制位为 111 的子数列 A0A_0A0 才能使得 i∩j=ii\cap j=ii∩j=i。而若当前位为 000，则两种序列都能取。

A=merge(A0+A1,A1)A=\text{merge}(A_0+A_1, A_1)
A=merge(A0+A1,A1)a=merge(a0−a1,a1)a=\text{merge}(a_0 - a_1, a_1)
a=merge(a0−a1,a1)

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

|  | void And(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j] += a[i + j + k] * type) %= P;  |
|  |  } |
|  |  } |
|  |  } |
|  | } |

按位异或

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

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

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

感性证明：发现这个新的运算符 ∘\circ∘ 其实就是 xxx 与 yyy 相同位数的奇偶性。若 (x∘y)⊕(x∘z)=0(x\circ y)\oplus (x\circ z)=0(x∘y)⊕(x∘z)=0，则 xxx 与 yyy、xxx 与 zzz 相同位数个数奇偶性相同，所以 y⊕zy\oplus zy⊕z 和 xxx 相同位数个数奇偶性也是相同的；若 (x∘y)⊕(x∘z)=1(x\circ y)\oplus (x\circ z)=1(x∘y)⊕(x∘z)=1，则 xxx 与 yyy、xxx 与 zzz 相同位数个数奇偶性不同，所以 y⊕zy\oplus zy⊕z 和 xxx 相同位数个数奇偶性也是不同的。

设 Ai=∑i∘j=0aj−∑i∘j=1ajA_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_jAi=∑i∘j=0aj−∑i∘j=1aj。我们来证一下 Ci=AiBiC_i=A_iB_iCi=AiBi 的正确性：

AiBi=(∑i∘j=0aj−∑i∘j=1aj)(∑i∘k=0bk−∑i∘k=1bk)=(∑i∘j=0aj∑i∘k=0bk+∑i∘j=1aj∑i∘k=1bk)−(∑i∘j=0aj∑i∘k=1bk+∑i∘j=1aj∑i∘k=0bk)=∑(j⊕k)∘i=0ajbk−∑(j⊕k)∘i=1ajbk=Ci\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}
AiBi=(i∘j=0∑aj−i∘j=1∑aj)(i∘k=0∑bk−i∘k=1∑bk)=(i∘j=0∑aji∘k=0∑bk+i∘j=1∑aji∘k=1∑bk)−(i∘j=0∑aji∘k=1∑bk+i∘j=1∑aji∘k=0∑bk)=(j⊕k)∘i=0∑ajbk−(j⊕k)∘i=1∑ajbk=Ci来看看怎么快速计算 A,BA,BA,B 的值，依旧是分治：

对于 iii 在当前位为 000 的子数列 A0A_0A0，进行 ∘\circ∘ 运算时发现它和 000 计算或和 111 计算结果都不会变（因为 0∩0=0,0∩1=00\cap 0=0,0\cap1=00∩0=0,0∩1=0），所以 Ai=∑i∘j=0aj−∑i∘j=1ajA_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_jAi=∑i∘j=0aj−∑i∘j=1aj 中的 ∑i∘j=1aj=0\sum_{i\circ j=1}a_j=0∑i∘j=1aj=0。

对于 iii 在当前位为 111 的子数列 A1A_1A1，进行 ∘\circ∘ 运算时发现它和 000 计算结果是 000，和 111 计算结果是 111（因为 1∩0=0,1∩1=11\cap 0=0,1\cap1=11∩0=0,1∩1=1）。

综上，有：

A=merge((A0+A1)−0,A0−A1)A=\text{merge}((A_0+A_1)-0, A_0-A_1)
A=merge((A0+A1)−0,A0−A1)也就是：

A=merge(A0+A1,A0−A1)A=\text{merge}(A_0+A_1, A_0-A_1)
A=merge(A0+A1,A0−A1)逆变换易得：

a=merge(a0+a12,a0−a12)a=\text{merge}(\frac{a_0+a_1}{2}, \frac{a_0-a_1}{2})
a=merge(2a0+a1,2a0−a1)给出代码，顺变换时 type=1\text{type}=1type=1，逆变换时 type=12\text{type}=\frac{1}{2}type=21。

|  | void Xor(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j] += a[i + j + k]) %= P;  |
|  |  (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;  |
|  |  (a[i + j] *= type) %= P; |
|  |  (a[i + j + k] *= type) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | } |

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

从另一个角度看待 FWT

我们设 c(i,j)c(i,j)c(i,j) 是 aja_jaj 对 AiA_iAi 的贡献系数。我们可以重新描述 FWT 变换的过程：

Ai=∑j=0n−1c(i,j)ajA_i = \sum_{j=0}^{n-1} c(i,j) a_j
Ai=j=0∑n−1c(i,j)aj因为有：

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

同时，ccc 函数还有一个重要的性质，它可以按位处理。

举个例子，我们变换的时候：

Ai=∑j=0n−1c(i,j)ajA_i = \sum_{j=0}^{n-1} c(i,j) a_j
Ai=j=0∑n−1c(i,j)aj这么做是比较劣的，我们将其拆分：

Ai=∑j=0(n−1)/2c(i,j)aj+∑j=(n−1)/2+1n−1c(i,j)ajA_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
Ai=j=0∑(n−1)/2c(i,j)aj+j=(n−1)/2+1∑n−1c(i,j)aj考虑前面的式子和后面的式子 i,ji,ji,j 的区别，发现只有最高位不同。

所以我们将 i,ji,ji,j 去除最高位的值为 i′,j′i’,j’i′,j′，并记 i0i_0i0 为 iii 的最高位。有：

Ai=c(i0,0)∑j=0(n−1)/2c(i′,j′)aj+c(i0,1)∑j=(n−1)/2+1n−1c(i′,j′)ajA_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
Ai=c(i0,0)j=0∑(n−1)/2c(i′,j′)aj+c(i0,1)j=(n−1)/2+1∑n−1c(i′,j′)aj如果 i0=0i_0=0i0=0，则有：

Ai=c(0,0)∑j=0(n−1)/2c(i′,j′)aj+c(0,1)∑j=(n−1)/2+1n−1c(i′,j′)ajA_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
Ai=c(0,0)j=0∑(n−1)/2c(i′,j′)aj+c(0,1)j=(n−1)/2+1∑n−1c(i′,j′)aji0=1i_0=1i0=1 则有：

Ai=c(1,0)∑j=0(n−1)/2c(i′,j′)aj+c(1,1)∑j=(n−1)/2+1n−1c(i′,j′)ajA_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
Ai=c(1,0)j=0∑(n−1)/2c(i′,j′)aj+c(1,1)j=(n−1)/2+1∑n−1c(i′,j′)aj也就是说，我们只需要：

[c(0,0)c(0,1)c(1,0)c(1,1)]\begin{bmatrix}
c(0,0) & c(0,1) \
c(1,0) & c(1,1)
\end{bmatrix}
[c(0,0)c(1,0)c(0,1)c(1,1)]四个数就可以完成变换了。我们称这个矩阵为位矩阵。

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

若逆矩阵为 c−1c^{-1}c−1，可以通过类似操作得到原数：

ai=∑j=0nc−1(i,j)Aja_i = \sum_{j=0}^n c^{-1}(i,j) A_j
ai=j=0∑nc−1(i,j)Aj逆矩阵不一定存在，比如如果有一排 000 或者一列 000 那么这个矩阵就没有逆，我们在构造时需要格外小心。

按位或

我们可以构造：

[1011]\begin{bmatrix}
1 & 0 \
1 & 1
\end{bmatrix}
[1101]这样满足 c(i,j)c(i,k)=c(i,j∪k)c(i,j)c(i,k)=c(i,j\cup k)c(i,j)c(i,k)=c(i,j∪k)。我们发现，这和我们前面推出的 A=merge(A0,A0+A1)A=\text{merge}(A_0, A_0+A_1)A=merge(A0,A0+A1) 一模一样！同理，下面也是一个满足这个条件的矩阵，但我们一般使用上面这个：

[1110]\begin{bmatrix}
1 & 1 \
1 & 0
\end{bmatrix}
[1110]虽然下面这个矩阵也满足 c(i,j)c(i,k)=c(i,j∪k)c(i,j)c(i,k)=c(i,j\cup k)c(i,j)c(i,k)=c(i,j∪k)，但这个矩阵存在一排 000，不存在逆，所以不合法：

[0011]\begin{bmatrix}
0 & 0 \
1 & 1
\end{bmatrix}
[0101]如果我们要进行逆变换，则需要对矩阵求逆，以最上面这个矩阵为例，得：

[10−11]\begin{bmatrix}
1 & 0 \
-1 & 1
\end{bmatrix}
[1−101]然后按照顺变换的方法，把逆变换矩阵代入即可。

按位与

我们可以构造：

[1101]\begin{bmatrix}
1 & 1 \
0 & 1
\end{bmatrix}
[1011]这样满足 c(i,j)c(i,k)=c(i,j∩k)c(i,j)c(i,k)=c(i,j\cap k)c(i,j)c(i,k)=c(i,j∩k)。

逆矩阵：

[1−101]\begin{bmatrix}
1 & -1 \
0 & 1
\end{bmatrix}
[10−11]### 按位异或

我们可以构造：

[111−1]\begin{bmatrix}
1 & 1 \
1 & -1
\end{bmatrix}
[111−1]这样满足 c(i,j)c(i,k)=c(i,j⊕k)c(i,j)c(i,k)=c(i,j\oplus k)c(i,j)c(i,k)=c(i,j⊕k)。

逆矩阵：

[0.50.50.5−0.5]\begin{bmatrix}
0.5 & 0.5 \
0.5 & -0.5
\end{bmatrix}
[0.50.50.5−0.5]## FWT 的性质

FWT 是线性变换。

若 FWT(X)FWT(X)FWT(X) 是 XXX 的 FWT 变换，则有：

FWT(A+B)=FWT(A)+FWT(B)FWT(A+B)=FWT(A)+FWT(B)
FWT(A+B)=FWT(A)+FWT(B)以及：

FWT(cA)=cFWT(A)FWT(cA)=cFWT(A)
FWT(cA)=cFWT(A)这样就可以实现快速卷积，参考第四道例题。

K 维 FWT

max 运算

我们重新看看我们的 ∪\cup∪ 运算，发现他实际上就是二进制下的取 max⁡\maxmax。我们将其拓展到 KKK 进制，有：

c(i,j)c(i,k)=c(i,max⁡(j,k))c(i,j)c(i,k)=c(i,\max(j,k))
c(i,j)c(i,k)=c(i,max(j,k))若 j=kj=kj=k，那么上式又是：

c(i,j)c(i,j)=c(i,j)c(i,j)c(i,j)=c(i,j)
c(i,j)c(i,j)=c(i,j)也就是说，每一行的 111 必定只能在 000 的前面，如果在后面则不合法了。手玩一下可以发现一组合法构造：

[1000110011101111]\begin{bmatrix}
1 & 0 & 0 & 0 \
1 & 1 & 0 & 0 \
1 & 1 & 1 & 0 \
1 & 1 & 1 & 1
\end{bmatrix}
1111011100110001求逆可得：

[1000−11000−11000−11]\begin{bmatrix}
1 & 0 & 0 & 0 \
-1 & 1 & 0 & 0 \
0 & -1 & 1 & 0 \
0 & 0 & -1 & 1
\end{bmatrix}
1−10001−10001−10001### min 运算

我们重新看看我们的 ∩\cap∩ 运算，发现他实际上就是二进制下的取 min⁡\minmin。我们将其拓展到 KKK 进制，有：

c(i,j)c(i,k)=c(i,min⁡(j,k))c(i,j)c(i,k)=c(i,\min(j,k))
c(i,j)c(i,k)=c(i,min(j,k))若 j=kj=kj=k，那么上式又是：

c(i,j)c(i,j)=c(i,j)c(i,j)c(i,j)=c(i,j)
c(i,j)c(i,j)=c(i,j)也就是说，每一行的 111 必定只能在 000 的后面，如果在后面则不合法了。手玩一下可以发现一组合法构造：

[1111011100110001]\begin{bmatrix}
1 & 1 & 1 & 1 \
0 & 1 & 1 & 1 \
0 & 0 & 1 & 1 \
0 & 0 & 0 & 1
\end{bmatrix}
1000110011101111求逆可得：

[1−10001−10001−10001]\begin{bmatrix}
1 & -1 & 0 & 0 \
0 & 1 & -1 & 0 \
0 & 0 & 1 & -1 \
0 & 0 & 0 & 1
\end{bmatrix}
1000−11000−11000−11前两者用得比较少，用得比较多的是：

不进位加法

我们重新看看我们的 ⊕\oplus⊕ 运算，发现他实际上就是二进制下的不进位加法。我们将其拓展到 KKK 进制，有：

c(i,j)c(i,k)=c(i,(j+k) mod K)c(i,j)c(i,k)=c(i,(j+k)\bmod K)
c(i,j)c(i,k)=c(i,(j+k)modK)我们构造 c(i,j)=ωKjc(i,j)=\omega_{K}^jc(i,j)=ωKj，就可以满足要求了：

ωKjωkk=ωK(j+k) mod K\omega_{K}^j\omega_{k}k=\omega_{K}^{(j+k)\bmod K}
ωKjωkk=ωK(j+k)modK但是每一行都一样矩阵也没有逆，所以我们可以构造 c(i,j)=ωK(i−1)jc(i,j)=\omega_{K}^{(i-1)j}c(i,j)=ωK(i−1)j 即可。

有下面这个矩阵：

[111⋯11ωK1ωK2⋯ωKk−11ωK2ωK4⋯ωK2(k−1)1ωK3ωK6⋯ωK3(k−1)⋮⋮⋮⋱⋮1ωKk−1ωK2(k−1)⋯ωKk(k−1)]\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}
1111⋮11ωK1ωK2ωK3⋮ωKk−11ωK2ωK4ωK6⋮ωK2(k−1)⋯⋯⋯⋯⋱⋯1ωKk−1ωK2(k−1)ωK3(k−1)⋮ωKk(k−1)这不就是我们熟悉的范德蒙德矩阵吗？

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

1K[111⋯11ωK−1ωK−2⋯ωK−(k−1)1ωK−2ωK−4⋯ωK−2(k−1)1ωK−3ωK−6⋯ωK−3(k−1)⋮⋮⋮⋱⋮1ωK−(k−1)ωK−2(k−1)⋯ωK−k(k−1)]\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}
K11111⋮11ωK−1ωK−2ωK−3⋮ωK−(k−1)1ωK−2ωK−4ωK−6⋮ωK−2(k−1)⋯⋯⋯⋯⋱⋯1ωK−(k−1)ωK−2(k−1)ωK−3(k−1)⋮ωK−k(k−1)如果我们题目给出的模数是存在单位根的，我们就可以简单实现，可以参考第六道例题。

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

[111⋯11x1x2⋯xk−11x2x4⋯x2(k−1)1x3x6⋯x3(k−1)⋮⋮⋮⋱⋮1xk−1x2(k−1)⋯xk(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}
1111⋮11x1x2x3⋮xk−11x2x4x6⋮x2(k−1)⋯⋯⋯⋯⋱⋯1xk−1x2(k−1)x3(k−1)⋮xk(k−1)但是这么做可能会存在零因子，也就是一个数有多种表示方法，我们无法确定一个数的真实值。

我们考虑不 (modxK−1)\pmod {x^K-1}(modxK−1) 了，我们 mod \bmodmod 分圆多项式 ΦK(x)\Phi_{K}(x)ΦK(x)，他满足 xxx 的阶为 kkk，且在 QQQ 上不可约。所以我们定义上面的计算是在 (modΦK(x))\pmod {\Phi_{K}(x)}(modΦK(x)) 下进行即可。

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

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

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

例题

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

求 ∪\cup∪、∩\cap∩、⊕\oplus⊕ 的三种卷积。

n≤17n\le17n≤17

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

|  | #include  |
|  | using namespace std; |
|  | #define ll long long |
|  | #define P 998244353 |
|  | const ll N = 1 << 18; |
|  | ll n; |
|  | ll A[N], B[N]; |
|  | ll a[N], b[N]; |
|  | void init() { |
|  | for(ll i = 0; i < n; i ++) a[i] = A[i], b[i] = B[i]; |
|  | } |
|  | void Or(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j + k] += a[i + j] * type) %= P;  |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | void And(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j] += a[i + j + k] * type) %= P;  |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | void Xor(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j] += a[i + j + k]) %= P;  |
|  |  (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;  |
|  |  (a[i + j] *= type) %= P; |
|  |  (a[i + j + k] *= type) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | void calc() { |
|  | for(ll i = 0; i < n; i ++) (a[i] *= b[i]) %= P; |
|  | } |
|  | void print() { |
|  | for(ll i = 0; i < n; i ++) printf("%lld ", (a[i] % P + P) % P); |
|  | printf("\n"); |
|  | } |
|  | int main() { |
|  | scanf("%lld", &n); |
|  |  n = 1 << n; |
|  | for(ll i = 0; i < n; i ++) scanf("%lld", &A[i]); |
|  | for(ll i = 0; i < n; i ++) scanf("%lld", &B[i]); |
|  |  |
|  | init(); Or(a, 1); Or(b, 1); calc(); Or(a, P - 1); print(); |
|  | init(); And(a, 1); And(b, 1); calc(); And(a, P - 1); print(); |
|  | init(); Xor(a, 1); Xor(b, 1); calc(); Xor(a, 499122177); print(); |
|  | } |

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

求：

ck=∑i∩j=0i∪j=kaibjc_k=\sum_{\substack{{i \cap j=0}\{i\cup j=k}}} a_i b_j
ck=i∩j=0i∪j=k∑aibjn≤20n\le20n≤20

首先，下半部分是我们喜闻乐见的 FWT 常见形式，而上半部分我们可以看成是 iii 与 jjj 不交。有：

i∪j=0⇒popcnt(i)+popcnt(j)=popcnt(i∪j)i\cup j=0\Rightarrow \text{popcnt}(i)+\text{popcnt}(j)=\text{popcnt}(i\cup j)
i∪j=0⇒popcnt(i)+popcnt(j)=popcnt(i∪j)所以我们可以构造：

Ai,j=∑i∪k=ipopcnt(j)=kaiA_{i,j}=\sum_{\substack{{i\cup k=i}\{\text{popcnt}(j)=k}}} a_i
Ai,j=i∪k=ipopcnt(j)=k∑ai可以枚举 popcnt\text{popcnt}popcnt 的值，分开考虑。

那么求 CCC 的时候有 Ci,j=∑j=0nAi,kBi,j−kC_{i,j}=\sum_{j=0}^n A_{i,k}B_{i,j-k}Ci,j=∑j=0nAi,kBi,j−k。

然后就可以做了。

|  | #include  |
|  | using namespace std; |
|  | #define popcnt(x) \_\_builtin\_popcountll(x) |
|  | #define ll long long |
|  | const ll M = 20, N = 1 << M, P = 1e9 + 9; |
|  | ll n, m; |
|  | ll a[M + 1][N], b[M + 1][N], c[M + 1][N]; |
|  | void Or(ll *a, ll type) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j + k] += a[i + j] * type) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | int main() { |
|  | scanf("%lld" ,&m); |
|  |  n = 1 << m; |
|  | for(ll i = 0; i < n; i ++) { |
|  | scanf("%lld", &a[popcnt(i)][i]); |
|  |  } |
|  | for(ll i = 0; i < n; i ++) { |
|  | scanf("%lld", &b[popcnt(i)][i]); |
|  |  } |
|  | for(ll i = 0; i <= m; i ++) { |
|  | Or(a[i], 1); |
|  | Or(b[i], 1); |
|  |  } |
|  | for(ll i = 0; i <= m; i ++) { |
|  | for(ll j = 0; j <= i; j ++) { |
|  | for(ll k = 0; k < n; k ++) { |
|  |  (c[i][k] += a[j][k] * b[i - j][k]) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | for(ll i = 0; i <= m; i ++) { |
|  | Or(c[i], -1); |
|  |  } |
|  | for(ll i = 0; i < n; i ++) { |
|  | printf("%lld ", (c[popcnt(i)][i] % P + P) % P); |
|  |  } |
|  | } |

「牛客 881D」Parity of Tuples:西部官网

给定 n×mn\times mn×m 的矩阵 aaa，定义 cnt(x)\text{cnt}(x)cnt(x) 为矩阵中有多少行对于 xxx 是合法的，合法的定义为这一行中每一个数 ai,j∩xa_{i,j}\cap xai,j∩x 的二进制值中都有奇数个 111。

你需要求出对于所有的 xxx，cnt\text{cnt}cnt 的取值。

n≤105,m≤10,x≤220n\le10^{5,m\le10,x\le2}{20}n≤105,m≤10,x≤220

再次重申，∩\cap∩ 是按位与的意思。

首先我们用数学公式定义一下 cnt\text{cnt}cnt（因为公式复杂，所以加了 large\tt largelarge）：

cnt(x)=12m∑i=1n∏j=1m(1−(−1)popcnt(ai,j∩x))\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)})
cnt(x)=2m1i=1∑nj=1∏m(1−(−1)popcnt(ai,j∩x))说明一下正确性。如果 popcnt(ai,j∩x)\text{popcnt}(a_{i,j}\cap x)popcnt(ai,j∩x) 是奇数的话，那么 (−1)popcnt(ai,j∩x)(-1)^{\text{popcnt}(a_{i,j}\cap x)}(−1)popcnt(ai,j∩x) 的结果就是 −1-1−1。最后 1−(−1)popcnt(ai,j∩x)1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}1−(−1)popcnt(ai,j∩x) 就是 222，最后会被 12m\frac{1}{2^m}2m1 除去；如果 popcnt(ai,j∩x)\text{popcnt}(a_{i,j}\cap x)popcnt(ai,j∩x) 是偶数的话，那么 (−1)popcnt(ai,j∩x)(-1)^{\text{popcnt}(a_{i,j}\cap x)}(−1)popcnt(ai,j∩x) 的结果就是 111。最后 1−(−1)popcnt(ai,j∩x)1-(-1)^{\text{popcnt}(a_{i,j}\cap x)}1−(−1)popcnt(ai,j∩x) 就是 000，那么整行的结果都是 000。

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

∏j=1m(1−(−1)popcnt(ai,j∩x))=(1−(−1)popcnt(ai,1∩x))(1−(−1)popcnt(ai,2∩x))⋯(1−(−1)popcnt(ai,m∩x))=1−∑a=1m(−1)popcnt(ai,a∩x)+∑a=1m∑b=a+1m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)−∑a=1m∑b=a+1m∑c=b+1m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)+popcnt(ai,c∩x)+⋯\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}
j=1∏m(1−(−1)popcnt(ai,j∩x))=(1−(−1)popcnt(ai,1∩x))(1−(−1)popcnt(ai,2∩x))⋯(1−(−1)popcnt(ai,m∩x))=1−a=1∑m(−1)popcnt(ai,a∩x)+a=1∑mb=a+1∑m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)−a=1∑mb=a+1∑mc=b+1∑m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)+popcnt(ai,c∩x)+⋯然后我们有一个性质：

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

容易得到：

=1−∑a=1m(−1)popcnt(ai,a∩x)+∑a=1m∑b=a+1m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)−∑a=1m∑b=a+1m∑c=b+1m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)+popcnt(ai,c∩x)+⋯=1−∑a=1m(−1)popcnt(ai,a∩x)+∑a=1m∑b=a+1m(−1)popcnt((ai,a⊕ai,b)∩x)−∑a=1m∑b=a+1m∑c=b+1m(−1)popcnt((ai,a⊕ai,b⊕ai,c)∩x)+⋯\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}
=1−a=1∑m(−1)popcnt(ai,a∩x)+a=1∑mb=a+1∑m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)−a=1∑mb=a+1∑mc=b+1∑m(−1)popcnt(ai,a∩x)+popcnt(ai,b∩x)+popcnt(ai,c∩x)+⋯=1−a=1∑m(−1)popcnt(ai,a∩x)+a=1∑mb=a+1∑m(−1)popcnt((ai,a⊕ai,b)∩x)−a=1∑mb=a+1∑mc=b+1∑m(−1)popcnt((ai,a⊕ai,b⊕ai,c)∩x)+⋯我们发现一加一减的可以容斥，我们容斥计算 fif_ifi 表示 nnn 行的所有式子中 (−1)i(-1)^i(−1)i 前面的系数和。

|  | // num 处理到第几列 |
|  | // x 当前的指数 |
|  | // mu 当前的系数（+1 or -1） |
|  | void dfs(ll *a, ll num, ll x, ll mu) { |
|  | if(num > m) { |
|  |  f[x] += mu; |
|  | return; |
|  |  } |
|  | dfs(a, num + 1, x, mu); // 不加入第 num 列，系数不变 |
|  | dfs(a, num + 1, x ^ a[num], -mu); |
|  | } |

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

=∑i=02k−1fx∩i(−1)x∩i\begin{align}
&= \sum_{i=0}^{2k-1} f_{x\cap i}(-1)^{x\cap i}
\end{align}
=i=0∑2k−1fx∩i(−1)x∩i我们突然发现后面这个 (−1)i(-1)^i(−1)i 取值只有两种，当 x∩ix\cap ix∩i 是奇数时取值为 −1-1−1，否则为 111。

好了，现在我们的问题转换为了求出：

∑popcnt(x∩i) mod 2=0fi−∑popcnt(x∩i) mod 2=1fi\sum_{\text{popcnt}(x\cap i)\bmod2=0} f_i-\sum_{\text{popcnt}(x\cap i)\bmod2=1} f_i
popcnt(x∩i)mod2=0∑fi−popcnt(x∩i)mod2=1∑fi这不就是 FWT 中的异或变换吗：

Ai=∑i∘j=0aj−∑i∘j=1ajA_i=\sum_{i\circ j=0}a_j-\sum_{i\circ j=1}a_j
Ai=i∘j=0∑aj−i∘j=1∑aj综上，我们发现这题就是推式子容斥之后得到 FWT 的形式。

原题需要将输出加密：

⨁x=02k−1(cnt(x)×3x mod (109+7))\bigoplus\limits_{x = 0}^{2k - 1} (\text{cnt}(x) \times 3^x \bmod (10^9 + 7))
x=0⨁2k−1(cnt(x)×3xmod(109+7))

|  | #include  |
|  | #define ll long long |
|  | using namespace std; |
|  | const ll P = 1e9 + 7; |
|  | #define N 100010 |
|  | #define M 20 |
|  | #define K 21 |
|  | ll n, m, k; |
|  | ll f[1 << K]; |
|  | ll a[N][M]; |
|  | // num 处理到第几列 |
|  | // x 当前的指数 |
|  | // mu 当前的系数（+1 or -1） |
|  | void dfs(ll *a, ll num, ll x, ll mu) { |
|  | if(num > m) { |
|  |  f[x] += mu; |
|  | return; |
|  |  } |
|  | dfs(a, num + 1, x, mu); // 不加入第 num 列，系数不变 |
|  | dfs(a, num + 1, x ^ a[num], -mu); |
|  | } |
|  | void Xor(ll *a, ll type) { |
|  | for(ll x = 2; x <= (1 << k); x <<= 1) { |
|  |  ll z = x >> 1; |
|  | for(ll i = 0; i < (1 << k); i += x) { |
|  | for(ll j = 0; j < z; j ++) { |
|  |  (a[i + j] += a[i + j + z]) %= P; |
|  |  (a[i + j + z] = a[i + j] - 2 * a[i + j + z]) %= P; |
|  |  (a[i + j] *= type) %= P; |
|  |  (a[i + j + z] *= type) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | ll qpow(ll x, ll y) { |
|  | if(y == 0) return 1; |
|  | if(y % 2 == 1) return x * qpow(x, y - 1) % P; |
|  |  ll tmp = qpow(x, y / 2); |
|  | return tmp * tmp % P; |
|  | } |
|  | int main() { |
|  | while(scanf("%lld %lld %lld", &n, &m, &k) != EOF) { |
|  | for(ll i = 0; i < (1 << k); i ++) f[i] = 0; |
|  | for(ll i = 1; i <= n; i ++) { |
|  | for(ll j = 1; j <= m; j ++) { |
|  | scanf("%lld", &a[i][j]); |
|  |  } |
|  | dfs(a[i], 1, 0, 1); |
|  |  } |
|  | Xor(f, 1); |
|  |  ll pw = 1, inv = qpow(1 << m, P - 2), ans = 0; |
|  | for(ll i = 0; i < (1 << k); i ++) { |
|  |  ans ^= f[i] * pw % P * inv % P; |
|  |  (pw *= 3) %= P; |
|  |  } |
|  | printf("%lld\n", ans); |
|  |  } |
|  | } |

「AT ABC212H」 Nim Counting

给定两个数 N,KN,KN,K，以及一个长度为 KKK 的整数数组 (A1,A2,⋯ ,AK)(A_1,A_2,\cdots, A_K)(A1,A2,⋯,AK)。

两个人玩 Nim 游戏。

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

任意选择一个 1≤M≤N1\le M\le N1≤M≤N，MMM 表示石子堆数。

对于每一堆，其石子数是 AAA 中任意一个数。

对于 ∑i=1NKi\sum_{i=1}^N K^i∑i=1NKi 种游戏，求先手获胜的游戏数，答案对 998244353998244353998244353 取模。

n≤2×105,K≤216,ai≤216n\le2\times10^5,K\le2{16},a_i\le2^{16}n≤2×105,K≤216,ai≤216

根据玩 Nim 游戏的经验，可以发现先手获胜当且仅当 ⨁i=0nAi≠0\bigoplus_{i=0}^n A_i\neq 0⨁i=0nAi=0。

所以我们定义 dp 式子 fi,jf_{i,j}fi,j 表示有 iii 个石堆，且石堆异或和为 jjj 的获胜方案数，有：

fi−1,j→∑k=1Kfi,j⊕akf_{i-1,j}\to \sum_{k=1}^Kf_{i,j\oplus a_k}
fi−1,j→k=1∑Kfi,j⊕ak答案就是 ∑i=1n∑j≠0fi,j\sum_{i=1}^n\sum_{j\neq0} f_{i,j}∑i=1n∑j=0fi,j。

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

fi,j=∑k⊕x=jfi−1,kaxf_{i,j}=\sum_{k\oplus x=j} f_{i-1,k}a_x
fi,j=k⊕x=j∑fi−1,kax我们定义 aaa 是一个全是 111 的数组即可。

同时，我们发现其实不需要真的进行 nnn 次卷积，其实只需要将 FWT 变换过之后的结果 AAA，求出 A+A2+A3+⋯+AnA+A^2+A3+\cdots+A^nA+A2+A3+⋯+An 即可。

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

|  | #include  |
|  | using namespace std; |
|  | #define ll long long |
|  | #define P 998244353 |
|  | const ll K = 1 << 20;  |
|  | ll n, k, ans; |
|  | ll f[K]; |
|  | void FWT(ll *a, ll type) { |
|  | for(ll x = 2; x <= K; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < K; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  (a[i + j] += a[i + j + k]) %= P; |
|  |  (a[i + j + k] = a[i + j] - 2 * a[i + j + k]) %= P; |
|  |  (a[i + j] *= type) %= P; |
|  |  (a[i + j + k] *= type) %= P; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | ll qpow(ll x, ll y) { |
|  | if(y == 0) return 1; |
|  | if(y % 2 == 1) return x * qpow(x, y - 1) % P; |
|  |  ll tmp = qpow(x, y / 2); |
|  | return tmp * tmp % P; |
|  | } |
|  | int main() { |
|  | scanf("%lld %lld", &n, &k); |
|  | for(ll i = 1; i <= k; i ++) { |
|  |  ll x; |
|  | scanf("%lld", &x); |
|  |  f[x] ++; |
|  |  } |
|  | FWT(f, 1); |
|  | for(ll i = 0; i < K; i ++) { |
|  | if(f[i] == 1) f[i] = n; |
|  | else { |
|  |  f[i] = f[i] * (qpow(f[i], n) - 1) % P * qpow(f[i] - 1, P - 2) % P; |
|  |  } |
|  |  } |
|  | FWT(f, 499122177); |
|  | for(ll i = 1; i < K; i ++) { |
|  |  (ans += f[i]) %= P; |
|  |  } |
|  | printf("%lld", (ans % P + P) % P); |
|  | } |

「AT ARC100E」 Or Plus Max

给你一个长度为 2n2^n2n 的序列 aaa，每个1≤K≤2n−11\le K\le 2^n-11≤K≤2n−1，找出最大的 ai+aja_i+a_jai+aj（i∪j≤Ki \cup j \le Ki∪j≤K，0≤i<j<2n0 \le i < j < 2^n0≤i<j<2n）并输出。

n≤18n\le18n≤18

就是要求 max⁡i∪j=kai+aj\max_{i\cup j=k} a_i+a_jmaxi∪j=kai+aj。

我们维护 fif_{i}fi 表示 max⁡i∪j=iai\max_{i\cup j=i} a_imaxi∪j=iai，gi=max2i∪j=iaig_i=\text{max2}_{i\cup j=i} a_igi=max2i∪j=iai，max2\text{max2}max2 表示次大值。

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

|  | #include  |
|  | using namespace std; |
|  | #define ll long long |
|  | const ll N = 1 << 18; |
|  | ll n; |
|  | struct node { |
|  |  ll mx1, mx2; |
|  | node(ll a = 0, ll b = 0):mx1(a), mx2(b) {} |
|  | friend node operator +(const node &x, const node &y) { |
|  | if(x.mx1 > y.mx1) { |
|  | return node(x.mx1, max(x.mx2, y.mx1)); |
|  |  } |
|  | return node(y.mx1, max(y.mx2, x.mx1)); |
|  |  } |
|  | } a[N]; |
|  | void FWT(node *a) { |
|  | for(ll x = 2; x <= n; x <<= 1) { |
|  |  ll k = x >> 1; |
|  | for(ll i = 0; i < n; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  |  a[i + j + k] = a[i + j] + a[i + j + k]; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | int main() { |
|  | scanf("%lld", &n); |
|  |  n = 1 << n; |
|  | for(ll i = 0; i < n; i ++) { |
|  | scanf("%lld", &a[i].mx1); |
|  |  } |
|  | FWT(a); |
|  | for(ll i = 0; i < n; i ++) { |
|  |  a[i].mx1 = a[i].mx1 + a[i].mx2; |
|  |  } |
|  |  ll ans = 0; |
|  | for(ll i = 1; i < n; i ++) { |
|  |  ans = max(ans, a[i].mx1); |
|  | printf("%lld\n", ans); |
|  |  } |
|  | } |

「HDU 6618」 Good Numbers

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

有多少个 kkk 位的 8 进制数（不含前导 0），满足这个数是好的，且是 ppp 的倍数。对 109+910^9+9109+9 取模。

例如：当 k=3,p=2k=3,p=2k=3,p=2 时，好数有 222,444,666222,444,666222,444,666 三个。

k≤1018,p<8k\le10^{18},p<8k≤1018,p<8

考虑状压 dp，设 fi,s,jf_{i,s,j}fi,s,j 表示第 iii 位，888 种数出现次数对 333 取模的状压情况，以及数对 ppp 取模的结果为 jjj。

答案就是 fk,0,0f_{k,0,0}fk,0,0。

直接暴力枚举位数转移是朴素的，瓶颈在于 kkk，考虑优化掉 kkk。

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

转移公式就是：

f2i,s1⊕s2,j1+t×j2←fi,s1,j1fi,s2,j2f_{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}
f2i,s1⊕s2,j1+t×j2←fi,s1,j1fi,s2,j2其中 ttt 是转移的位数，而 ⊕\oplus⊕ 在这里是不进位三进制加法。

发现这样多了瓶颈——我们需要枚举 s1s_1s1 和 s2s_2s2。

但是我们发现这不就是 FWT 中异或的形式吗：ci⊕j←aibjc_{i\oplus j}\gets a_ib_jci⊕j←aibj。考虑三进制 FWT 加速。下面给出 FWT 的代码，w1 是原根的一次方，w2 是原根的二次方：

|  | void FWT(ll *a, ll type) { |
|  | for (ll x = 3; x <= N; x *= 3) { |
|  |  ll k = x / 3; |
|  | for (ll i = 0; i < N; i += x) { |
|  | for (ll j = 0; j < k; j ++) { |
|  | for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k]; |
|  | if (type == 1) { |
|  |  tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; |
|  |  tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; |
|  |  tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; |
|  |  } else { |
|  |  tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; |
|  |  tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; |
|  |  tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; |
|  | for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P; |
|  |  } |
|  | for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l]; |
|  |  } |
|  |  } |
|  |  } |
|  | } |

因为 1e9+91e9+91e9+9 存在原根 222，然后就朴素实现了，注意位矩阵：

[1111ω31ω321ω32ω34]\begin{bmatrix}
1 & 1 & 1 \
1 & \omega_{3}^1 & \omega_{3}^2 \
1 & \omega_{3}^2 & \omega_{3}^4 \
\end{bmatrix}
1111ω31ω321ω32ω34代码：

|  | #include  |
|  | using namespace std; |
|  | #define ll long long |
|  | const ll P = 1e9 + 9; |
|  | ll qpow(ll x, ll y) { |
|  | if(y == 0) return 1; |
|  | if(y % 2 == 1) return x * qpow(x, y - 1) % P; |
|  |  ll tmp = qpow(x, y / 2); |
|  | return tmp * tmp % P;  |
|  | } |
|  | const ll G = 2; |
|  | const ll w1 = qpow(G, (P - 1) / 3); |
|  | const ll w2 = qpow(G, (P - 1) / 3 * 2); |
|  | const ll inv3 = qpow(3, P - 2); |
|  | const ll N = 3 * 3 * 3 * 3 * 3 * 3 * 3 * 3; |
|  | ll n, p; |
|  | ll tmp[8][N], res[8][N], one[8][N]; |
|  | ll a[8][N], b[8][N]; |
|  | ll pw3[8]; |
|  | ll tmp1[3], tmp2[3]; |
|  | void FWT(ll *a, ll type) { |
|  | for (ll x = 3; x <= N; x *= 3) { |
|  |  ll k = x / 3; |
|  | for (ll i = 0; i < N; i += x) { |
|  | for (ll j = 0; j < k; j ++) { |
|  | for (ll l = 0; l < 3; l++) tmp1[l] = a[i + j + l * k]; |
|  | if (type == 1) { |
|  |  tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; |
|  |  tmp2[1] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; |
|  |  tmp2[2] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; |
|  |  } else { |
|  |  tmp2[0] = (tmp1[0] + tmp1[1] + tmp1[2]) % P; |
|  |  tmp2[1] = (tmp1[0] + tmp1[1] * w2 + tmp1[2] * w1) % P; |
|  |  tmp2[2] = (tmp1[0] + tmp1[1] * w1 + tmp1[2] * w2) % P; |
|  | for (ll l = 0; l < 3; l++) (tmp2[l] *= inv3) %= P; |
|  |  } |
|  | for (ll l = 0; l < 3; l++) a[i + j + l * k] = tmp2[l]; |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | ll base = 1; |
|  | void fun(ll x) { |
|  | if(x == 1) { |
|  | memset(res, 0, sizeof res); |
|  | memset(tmp, 0, sizeof tmp); |
|  | memset(one, 0, sizeof one); |
|  | for(ll i = 1; i < 8; i ++) res[i % p][pw3[i]] = 1; |
|  | for(ll i = 0; i < 8; i ++) tmp[i % p][pw3[i]] = 1; |
|  | for(ll i = 0; i < 8; i ++) one[i % p][pw3[i]] = 1; |
|  | for (int i = 0; i < p; i ++) { |
|  | FWT(tmp[i], 1); |
|  | FWT(res[i], 1); |
|  | FWT(one[i], 1); |
|  |  } |
|  |  base = 8 % p; |
|  | return; |
|  |  } |
|  | if(x % 2 == 1) { |
|  | fun(x - 1); |
|  | memset(a, 0, sizeof a); |
|  | memset(b, 0, sizeof b); |
|  | for (ll i = 0; i < p; i ++) { |
|  | for (ll j = 0; j < p; j ++) { |
|  |  ll k = (i * 8 + j) % p; |
|  | for (ll x = 0; x < N; x ++) |
|  |  (a[k][x] += tmp[i][x] * one[j][x]) %= P, |
|  |  (b[k][x] += res[i][x] * one[j][x]) %= P; |
|  |  } |
|  |  } |
|  | memcpy(tmp, a, sizeof a); |
|  | memcpy(res, b, sizeof b); |
|  |  (base *= 8) %= P; |
|  | return; |
|  |  } |
|  | fun(x / 2); |
|  | memset(a, 0, sizeof a); |
|  | memset(b, 0, sizeof b); |
|  | for (ll i = 0; i < p; i ++) { |
|  | for (ll j = 0; j < p; j ++) { |
|  |  ll k = (i * base + j) % p; |
|  | for (ll x = 0; x < N; x ++) |
|  |  (a[k][x] += tmp[i][x] * tmp[j][x]) %= P, |
|  |  (b[k][x] += res[i][x] * tmp[j][x]) %= P; |
|  |  } |
|  |  } |
|  | memcpy(tmp, a, sizeof a); |
|  | memcpy(res, b, sizeof b); |
|  |  (base *= base) %= p; |
|  | } |
|  | int main() { |
|  |  pw3[0] = 1; |
|  | for(ll i = 1; i <= 8; i ++) { |
|  |  pw3[i] = pw3[i - 1] * 3; |
|  |  } |
|  | while(scanf("%lld %lld", &n, &p) != EOF) { |
|  | fun(n); |
|  | FWT(res[0], -1); |
|  | printf("%lld\n", res[0][0]); |
|  |  } |
|  | } |

「CF 1103E」Radix sum

给定一个长度为 nnn 的序列 a1,a2,…,ana_1,a_2,…,a_na1,a2,…,an，对于每一个 p∈[0,n−1]p \in [0,n-1]p∈[0,n−1]，求满足下列条件的整数序列 i1,i2,…,ini_1,i_2,…,i_ni1,i2,…,in 的方案数，对 2582^{58}258 取模：

∀j∈[1,n],ij∈[1,n]\forall j \in [1,n] , i_j \in [1,n]∀j∈[1,n],ij∈[1,n]；
∑j=1naij=p\sum\limits_{j=1}^n a_{i_j} = pj=1∑naij=p，这里的加法定义为十进制不进位加法。

n≤105,ai≤105n\le10^5,a_i\le105n≤105,ai≤105

我们可以想到 dp：设计状态 fi,sf_{i,s}fi,s 表示考虑到第 iii 个数，当前加法状态为 sss。因为 FWT 变换时线性的，可以先变换为 FWT 点值表示法，然后变成自己的 nnn 次幂，最后再变换回来。

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

这里的分圆多项式 Φ10(x)=x4−x3+x2−x+1\Phi_{10}(x)=x^4-x3+x^2-x+1Φ10(x)=x4−x3+x2−x+1。

然而我们发现 IFWT 时，需要除去进制 101010，然而我们发现 101010 在 2582^{58}258 下没有逆元。实际上我们发现 555 在 2582^{58}258 下是有逆元的：576460752303423495764607523034234957646075230342349，我们只需要再除去一个 222 就可以了。设已经除以了 555 的答案为 xxx，真正的答案为 yyy，也就是 25y≡x(mod264)2^5y\equiv x\pmod{2^{64}}25y≡x(mod264)，显然，我们有 y≡x25(mod264−5)y\equiv \frac{x}{2^5}\pmod{2{64-5}}y≡25x(mod264−5)，也就是 y≡x25(mod259)y\equiv \frac{x}{2^5}\pmod{2{59}}y≡25x(mod259)，所以直接将最后的答案除以 252^525 即可。虽然出题人不知道为什么要模 2582^{58}258，但再取下模即可。

然后就是平凡实现了：

|  | #include  |
|  | using namespace std; |
|  | #define ll unsigned long long |
|  | const ll P = 1ull << 58, N = 1e5 + 10; |
|  | const ll m = 5, K = 10; |
|  | ll inv5; |
|  | ll n; |
|  | ll pw[m + 1]; |
|  | ll qpow(ll x, ll y) { |
|  | if(y == 0) return 1; |
|  | if(y % 2 == 1) return x * qpow(x, y - 1); |
|  |  ll tmp = qpow(x, y / 2); |
|  | return tmp * tmp; |
|  | } |
|  | struct poly { |
|  |  ll a[30]; |
|  | poly() {memset(a, 0, sizeof a);} |
|  |  ll operator [](ll x) const {return a[x];} |
|  |  ll& operator [](ll x) {return a[x];} |
|  | friend poly operator *(const poly &x, const poly &y) { |
|  |  poly z; |
|  | for(ll i = 0; i < K; i ++) { |
|  | for(ll j = 0; j < K; j ++) { |
|  |  z[(i + j) % K] += x[i] * y[j]; |
|  |  } |
|  |  } |
|  | return z; |
|  |  } |
|  | friend poly operator *(const poly &x, const ll &y) { |
|  |  poly z; |
|  | for(ll i = 0; i < K; i ++) { |
|  |  z[i] += x[i] * y; |
|  |  } |
|  | return z; |
|  |  } |
|  | friend poly operator +(const poly &x, const poly &y) { |
|  |  poly z; |
|  | for(ll i = 0; i < K; i ++) { |
|  |  z[i] += x[i] + y[i]; |
|  |  } |
|  | return z; |
|  |  } |
|  | poly w(ll x) { |
|  |  poly res; |
|  | for(ll i = 0; i < K; i ++) { |
|  |  res[(i + x) % K] += a[i]; |
|  |  } |
|  | return res; |
|  |  } |
|  | } T, f[N], one; |
|  | poly qpow(poly x, ll y) { |
|  | if(y == 0) return one; |
|  | if(y % 2 == 1) return x * qpow(x, y - 1); |
|  |  poly tmp = qpow(x, y / 2); |
|  | return tmp * tmp; |
|  | } |
|  | poly tmp1[30], tmp2[30]; |
|  | void FWT(poly *a, ll type) { |
|  | for(ll x = K; x <= pw[m]; x *= K) { |
|  |  ll k = x / K; |
|  | for(ll i = 0; i < pw[m]; i += x) { |
|  | for(ll j = 0; j < k; j ++) { |
|  | for(ll l = 0; l < K; l ++) tmp1[l] = a[i + j + l * k], tmp2[l] = poly(); |
|  | if(type == 1) { |
|  | for(ll l = 0; l < K; l ++) { |
|  | for(ll v = 0; v < K; v ++) { |
|  |  tmp2[l] = tmp2[l] + tmp1[v].w(l * v % K); |
|  |  } |
|  |  } |
|  | for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l]; |
|  |  } else { |
|  | for(ll l = 0; l < K; l ++) { |
|  | for(ll v = 0; v < K; v ++) { |
|  |  tmp2[l] = tmp2[l] + tmp1[v].w((K - (l * v % K)) % K); |
|  |  } |
|  |  } |
|  | for(ll l = 0; l < K; l ++) a[i + j + l * k] = tmp2[l] * inv5; |
|  |  } |
|  |  } |
|  |  } |
|  |  } |
|  | } |
|  | ll mod(poly x){ |
|  |  ll n = 4; |
|  | for(ll i = K - 1; i >= n; i --){ |
|  |  ll u = x[i]; |
|  | for(ll j = 1; j <= n; j ++) x[i - j] -= u * T[n - j]; |
|  |  } |
|  |  ll u = x[0]; |
|  |  u >>= m; |
|  | return u % P; |
|  | } |
|  | int main() { |
|  |  pw[0] = 1; |
|  | for(ll i = 1; i <= m; i ++) pw[i] = pw[i - 1] * K; |
|  |  T[0] = 1, T[1] = -1, T[2] = 1, T[3] = -1, T[4] = 1; // 分圆多项式phi10 |
|  |  one[0] = 1; |
|  |  inv5 = 57646075230342349ull; |
|  | scanf("%llu", &n); |
|  | for(ll i = 1; i <= n; i ++) { |
|  |  ll x; |
|  | scanf("%llu", &x); |
|  |  f[x][0] ++; |
|  |  } |
|  | FWT(f, 1); |
|  | for(ll i = 0; i < pw[m]; i ++) f[i] = qpow(f[i], n); |
|  | FWT(f, -1); |
|  | for(ll i = 0; i < n; i ++) cout<<mod(f[i])<<'\n'; |
|  | } |