快速沃尔什变换

简介

沃尔什转换（Walsh Transform）是在频谱分析上作为离散傅立叶变换的替代方案的一种方法。——维基百科

其实这个变换在信号处理中应用很广泛，FFT 是 double 类型的，但是 Walsh 把信号在不同震荡频率方波下拆解，因此所有的系数都是绝对值大小相同的整数，这使得不需要作浮点数的乘法运算，提高了运算速度。

所以，FWT 和 FFT 的核心思想应该是相同的，都是对数组的变换。我们记对数组进行快速沃尔什变换后得到的结果为。

那么 FWT 核心思想就是：

我们需要一个新序列，由序列和序列经过某运算规则得到，即 $C = A \cdot B$ ；

我们先正向得到序列，再根据 $FWT[C]=FWT[A] \cdot FWT[B]$ 在的时间复杂度内求出，其中 $\cdot$ 是序列对应位置相乘；

然后逆向运算得到原序列。时间复杂度为 $O(n \log{n})$ 。

在算法竞赛中，FWT 是用于解决对下标进行位运算卷积问题的方法。

公式： $C_{i} = \sum_{i=j \oplus k}A_{j} B_{k}$

（其中 $\oplus$ 是二元位运算中的某一种）

下面我们举 $\cup$ （按位或）、 $\cap$ （按位与）和 $\oplus$ （按位异或）为例。

FWT 的运算

或运算

如果有 $k=i\cup j$ ，那么的二进制位为的位置和的二进制位为的位置肯定是的二进制位为的位置的子集。

现在要得到 $FWT[C] = FWT[A] \cdot FWT[B]$ ，我们就要构造这个 FWT 的规则。

我们按照定义，显然可以构造 $FWT[A]_i = A'_i = \sum_{i=i\cup j}A_{j}$ ，来表示满足二进制中为的子集。

那么有：

\begin{aligned} FWT[A]_i\cdot FWT[B]_i&=\left(\sum_{i\cup j=i} A_j\right)\left(\sum_{i\cup k=i} B_k\right) \\ &=\sum_{i\cup j=i}\sum_{i\cup k=i}A_jB_k \\ &=\sum_{i\cup(j\cup k)=i}A_jB_k \\ &= FWT[C]_i \end{aligned}

那么我们接下来看怎么求。

首先肯定不能枚举了，复杂度为。既然不能整体枚举，我们就考虑分治。

我们把整个区间二分，其实二分区间之后，下标写成二进制形式是有规律可循的。

我们令表示的前一半，表示区间的后一半，那么就是 A 下标最大值的最高位为，他的子集就是他本身的子集（因为最高位为了），但是的最高位是，他满足条件的子集不仅仅是他本身，还包最高位为的子集，即

其中 merge 表示像字符串拼接一样把两个数组拼起来，就是普通加法，表示对应二进制位相加。

这样我们就通过二分能在 $O(\log{n})$ 的时间复杂度内完成拼接，每次拼接的时候要完成一次运算，也就是说在 $O(n\log{n})$ 的时间复杂度得到了。

接下来就是反演了，其实反演是很简单的，既然知道了的本身的子集是他自己（），的子集是，那就很简单的得出反演的递推式了：

下面我们给出代码实现。容易发现顺变换和逆变换可以合并为一个函数，顺变换时 $\text{type}=1$ ，逆变换时 $\text{type}=-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;
      }
    }
  }
}

与运算

与运算类比或运算可以得到类似结论

下面我们给出代码实现。顺变换时 $\text{type}=1$ ，逆变换时 $\text{type}=-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;
      }
    }
  }
}

异或运算

异或的卷积是基于如下原理：

若我们令 $x\circ y$ 表示 $x\cap y$ 中数量的奇偶性，即 $x\circ y=\text{popcnt}(x\cap y)\bmod 2$ ，那么容易有 $(x\circ y)\oplus (x\circ z)=x\circ(y\oplus z)$ 。

对于的运算其实也很好得到。

设 $FWT[A]_i=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j$ 。我们来证一下 $FWT[C] = FWT[A] \cdot FWT[B]$ 的正确性：

\begin{aligned} FWT[A]_iFWT[B]_i&=\left(\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j\right)\left(\sum_{i\circ k=0}B_k-\sum_{i\circ k=1}B_k\right) \\ &=\left(\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\right)-\left(\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\right) \\ &=\sum_{(j\oplus k)\circ i=0}A_jB_k-\sum_{(j\oplus k)\circ i=1}A_jB_k \\ &=FWT[C]_i \end{aligned}

来看看怎么快速计算的值，依旧是分治：

对于在当前位为的子数列，进行 $\circ$ 运算时发现它和计算或和计算结果都不会变（因为 $0\cap 0=0,0\cap1=0$ ），所以 $FWT[A]=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j$ 中的 $\sum_{i\circ j=1}A_j=0$ 。

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

综上，有：

也就是：

逆变换易得：

UFWT[A'] = merge(\frac{UFWT[A_0'] + UFWT[A_1']}{2}, \frac{UFWT[A_0'] - UFWT[A_1']}{2})

给出代码，顺变换时 $\text{type}=1$ ，逆变换时 $\text{type}=\frac{1}{2}$ 。

实现

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[A]_{i} = \sum_{C_1}A_{j} - \sum_{C_2}A_{j}$ （表示 $\text{popcnt}(x\cup y)\bmod 2$ 为，表示 $\text{popcnt}(x\cup y)\bmod 2$ 为）

UFWT[A'] = merge(\frac{UFWT[A_1'] - UFWT[A_0']}{2}, \frac{UFWT[A_1'] + UFWT[A_0']}{2})

另一个角度的 FWT

我们设是对的贡献系数。我们可以重新描述 FWT 变换的过程：

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

因为有：

FWT[A]_i\cdot FWT[B]_i=FWT[C]_i

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

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

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

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

这么做是比较劣的，我们将其拆分：

FWT[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

考虑前面的式子和后面的式子的区别，发现只有最高位不同。

所以我们将去除最高位的值为，并记为的最高位。有：

FWT[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

如果，则有：

FWT[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

则有：

FWT[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) FWT[A]_j

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

按位或

我们可以构造：

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

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

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

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

\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[c\cdot A]=c\cdot FWT[A]

K 维 FWT

其实位运算的本质是对一个维 $\{0,1\}$ 向量的运算。或运算就是每一维取 $\max$ 。且运算就是每一维取 $\min$ 。异或运算则是每一维对应相加再 $\bmod 2$ 。

位运算有个特点：向量的每一位都是独立的。

我们把 $\{0,1\}$ 扩展到 $[0,K)\cap \mathbf{Z}$ 也就是扩展到进制，看看会得到什么？

max 运算

我们将 $\cup$ 运算拓展到进制，定义 $i\cup j$ 表示按位取 $\max$ ，有：

c(i,j)c(i,k)=c(i,j\cup k)

若，那么上式又是：

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

\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$ 运算拓展到进制，定义 $i\cap j$ 表示按位取 $\min$ ，有：

c(i,j)c(i,k)=c(i,j\cap k)

若，那么上式又是：

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

\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$ 运算拓展到进制，定义 $i\oplus j$ 表示按位相加再 $\bmod K$ ，有：

c(i,j)c(i,k)=c(i,j\oplus k)

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

\omega_{K}^j\omega_{K}^k=\omega_{K}^{j\oplus 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-1)(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-1)(k-1)} \end{bmatrix}

如果我们题目给出的模数是存在单位根的，我们就可以简单实现。

但是 单位根在模意义下可能不存在，所以我们考虑扩域，就是人为地定义一个，满足，然后直接把代入计算，这样每个数都是一个关于的次多项式。我们只需要在 $\bmod {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-1)(k-1)} \end{bmatrix}

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

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

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

例题

「CF 1103E」Radix sum

给定一个长度为的序列，对于每一个 $p \in [0,n-1]$ ，求满足下列条件的整数序列的方案数，对 $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}$ 表示考虑到第个数，当前加法状态为。因为 FWT 变换是线性的，可以先变换为 FWT 点值表示法，然后变成自己的次幂，最后再变换回来。

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

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

然而我们发现 UFWT 时，需要除去进制，然而我们发现在 $2^{58}$ 下没有逆元。实际上我们发现在 $2^{58}$ 下是有逆元的：，我们只需要再除去一个就可以了。设已经除以了的答案为，真正的答案为，也就是 $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^{58}$ ，但再取下模即可。

【CF103329F】【XXII Opencup, Grand Prix of XiAn】The Struggle

给出一个椭圆，其中所有整点的坐标均在 $[1,4 \cdot 10^6]$ 之间。求 $\sum_{(x,y) \in E} (x \oplus y)^{33}x^{-2}y^{-1} \mod 10^9+7$ 的值。

题解

这是一道比较不裸的题，出题人提供了详细的英文题解，具体请见此链接。

参考资料

本页面最近更新：2024/6/22 15:45:13，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：Xeonacid, ZnPdCo, BinDir0, Enter-tainer, Ir1d, lxlonlyn, nocriz, nocrizwang, Tiphereth-A, TrisolarisHD, xyf007
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用