学习笔记：快速沃尔什变换

参考文献

快速变换之Fast Walsh Hadamard Transform

前置知识

1. FFT（快速傅里叶变换）
2. 卷积定理：函数卷积的傅里叶变换是函数傅里叶变换的乘积
3. 正交：是线性代数的概念，是垂直这一直观概念的推广。 作为一个形容词，只有在一个确定的内积空间中才有意义。 若内积空间中两向量的内积为0，则称它们是正交的。 如果能够定义向量间的夹角，则正交可以直观的理解为垂直。
4. 二次剩余：在数论中，特别在同余理论里，一个整数 $X$ 对另一个整数 $P$ 的二次剩余（英语：Quadratic residue）指 $X$ 的平方 $X^2$ 除以 $P$ 得到的余数。

引入

首先，我们知道多项式乘法可以通过 FFT 降低复杂度。

即通过 FFT，能够在 $O(n\log{n})$ 的时间复杂度求下面的循环卷积：

$$c_k=\sum_{i+j=k}a_ib_j$$

由卷积定理保证了：

$$DFT(c)=DFT(a)·DFT(b)$$

因此如果要求 $a·b$，可以先利用 FFT进行转换，再点积，进行逆转换：

$$c=IDFT(DFT(a)·DFT(b))$$

但是，假如我们需要将卷积转化为与/或/异或卷积，也就是求解：

$$c_k=\sum_{i*j=k}a_ib_j$$

其中 $*$ 是与/或/异或运算之一，是否存在一种线性变换满足如下相同的卷积定理呢？

$$A(c)=A(a)·A(b)$$

并且，正如 FFT ，该变换也要有快速算法。答案是肯定的，即沃尔什变换。

每种运算对应一种线性变换，其中异或运算就对应 WHT，沃尔什-阿达玛变换。

简介

沃尔什-阿达玛变换（Walsh-Hadamard Transform）是一种广义傅里叶变换，是信号处理，集成电路和图像处理中频谱分析的一种变换方法，用于替换离散傅里叶变换。快速沃尔什-阿达玛变换（FWHT，或者通常在国内算法竞赛界被称为 FWT）是一种进行 WHT 的快速算法，和 FFT 类似。

傅里叶变换是浮点类型（还是复数）的，把信号在谐波下进行分解，而 Walsh 把信号在不同震荡频率方波下分解，因此所有的系数都是绝对值大小相同的整数，这使得不需要作浮点数的乘法运算，提高了运算速度。

这里谐波指的就是单位正弦波或者余弦波的倍数频率的波，比如 $\sin(kx)$。而方波就是矩形的波。

如果用数学公式描述，就是使用 Hadamard 矩阵对信号向量进行变换，以下是四阶和八阶的 Hadamard 矩阵：

$$H_4=\begin{bmatrix}1&1&1&1\\1&-1&1&-1\\1&1&-1&-1\\1&-1&-1&1\end{bmatrix}$$ $$H_8=\begin{bmatrix}1&1&1&1&1&1&1&1\\1&-1&1&-1&1&-1&1&-1\\1&1&-1&-1&1&1&-1&-1\\1&-1&-1&1&1&-1&-1&1\\1&1&1&1&-1&-1&-1&-1\\1&-1&1&-1&-1&1&-1&1\\1&1&-1&-1&-1&-1&1&1\\1&-1&-1&1&-1&1&1&-1\end{bmatrix}$$

Hadamard 矩阵

Hadamard 矩阵即一个只包含 $1$ 和 $-1$ 的矩阵，满足每一行/列都相互正交。

$1$，$2$ 以及 $4^k$，$2*4^k$ 阶的 Hadamard 矩阵一定存在，可证明 Hadamard 矩阵的阶数除了 $1$，$2$ 外都是 $4$ 的倍数。

因而有 Hadamard 猜想：$4n$ 阶矩阵一定存在。

一般 $4n$ 阶的 Hadarmard 矩阵存在性问题以及构造，依旧是个未解决难题。关于 Hadarmard 矩阵，可以看张贤达的《矩阵分析》一书以及冯荣权宋春伟的《组合数学》有一些结论和比较旧的相关文献。

西尔维斯特（Sylvester）给出了由一个 Hadarmard 矩阵构造一个新的 Hadarmard 矩阵的方法，若 $H$ 是一个 Hadarmard 矩阵，则容易证明下面矩阵也是 Hadarmard 矩阵

$$\begin{bmatrix}H&H\\H&-H\end{bmatrix}$$

Arasu 等人的 Hadamard and Conference Matrices 中提出了构造 Hadamard 的另一种方法。

Paley 在 1933 年利用二次剩余的理论发现了一种更有效的构造 Hadamard 矩阵。

沃尔什·阿达玛变换（Walsh-Hadamard Transform）

$2^k$ 阶的 Walsh-Hadamard Transform 定义为

$$X=\frac{1}{2^k}H_{2^k}x\,\,\,\text{或者}\,\,\,X=H_{2^k}x$$

显然逆变换就是

$$x=H_{2^k}^{T}X\,\,\,\text{或者}\,\,\,x=\frac{1}{2^k}H_{2^k}^{T}X$$

此处使用后者的定义，$H_{2^k}$ 是通过西尔维斯特构造法得到的 Hadarmard 矩阵。

$$H_{2^k}=\begin{bmatrix}H_{2^{k-1}}&H_{2^{k-1}}\\H_{2^{k-1}}&-H_{2^{k-1}}\end{bmatrix}$$

容易发现 $H_{2^k}$ 是一个对称正交阵，所以 $H_{2^k}^T=H_{2^k}$。

沃尔什变换可以展开成显式的求和表达式。注意到 $H_1=\begin{bmatrix}1\end{bmatrix}$，则 $H_{2^k}$ 的 i 行 j 列为

$$H_{2^k}=(-1)^{[i\geq2^{k-1}\text{and}j\geq2^{k-1}]}H_{2^{k-1}(i\%2^{k-1},j\%2^{k-1})}$$

记 $pc(a)$ 为 $a$ 中二进制 $1$ 的个数，$[expression]$ 当 $expression$ 为真时为 $1$，否则为 $0$。则可以得到

$$H_{2^k}=(-1)^{pc(i\&j)}$$

因此可以将 Walsh 变换写成

$$H_nx=\left(\sum(-1)^{pc(i\&0)}x_i,\sum(-1)^{pc(i\&1)}x_i,\cdots,\sum(-1)^{pc(i\&(n-1))}x_i \right)$$

其中 $n=2^k$ 的形式较为常见，或者简写成

$$WT[x]_i=(H_{2^k}x)_i=\sum_{i\oplus j=0}x_j-\sum_{i\oplus j=1}$$

$WT[x]_i$ 表示 $x$ 经过 Walsh 变换后的第 $i$ 位。

快速沃尔什变换

我们可以将 $x$ 进行拆分，$x=[x_0,x_1]$。因而

$$H_{2^k}x=\begin{bmatrix}H_{2^{k-1}}&H_{2^{k-1}}\\H_{2^{k-1}}&-H_{2^{k-1}}\end{bmatrix}=\begin{bmatrix}H_{2^{k-1}}x_0+H_{2^{k-1}}x_1\\H_{2^{k-1}}x_0-H_{2^{k-1}}x_1\end{bmatrix}$$

故只需要计算 $H_{2^{k-1}}x_0$ 和 $H_{2^{k-1}}x_1$，按照这个过程继续分解就能得到快速沃尔什变换。

逆变换同理：

$$\frac{1}{2^k}H_{2^k}^TX=\frac{1}{2^k}\begin{bmatrix}H_{2^{k-1}}^TX_0+H_{2^{k-1}}^TX_1\\H_{2^{k-1}}^TX_0-H_{2^{k-1}}^TX_1\end{bmatrix}=\frac{1}{2}\begin{bmatrix}\frac{1}{2^{k-1}}H_{2^{k-1}}^TX_0+\frac{1}{2^{k-1}}H_{2^{k-1}}^TX_1\\\frac{1}{2^{k-1}}H_{2^{k-1}}^TX_0-\frac{1}{2^{k-1}}H_{2^{k-1}}^TX_1\end{bmatrix}$$

也可以写成如下形式

$$FWT(X)=\text{merge}(FWT(X_0)+FWT(X_1),FWT(X_0)-FWT(X_1))\\ IFWT(X)=\text{merge}\left(\frac{IFWT(X_0)+IFWT(X_1)}{2},\frac{IFWT(X_0)-IFWT(X_1)}{2}\right)$$

其中 FWT , IFWT 是正逆 Walsh 变换，merge 是向量拼接。这种形式在网上比较常见。

Xor(dyadic) 卷积和卷积定理

对两个长度为 $2^k$ 的向量 $x,y$，定义 Xor/dyadic 卷积如下

$$z(n)=\sum_{d=0}^{2^k-1}x(d)y(n\oplus d)$$

其中 $n<2^k$，或者等价形式

$$z(n)=\sum_{a\oplus b=n}x(a)y(b)$$

其中$\oplus$为异或，$x(a)$ 表示 $x$ 第 $a$ 位。显然这个卷积满足交换律，乘法分配律，结合律。

按照上面定义的 Xor 卷积有如下的卷积定理

$$z=x*y\Rightarrow Z=XY$$

对所有 $2^k$ 阶的 WHT 变换都成立。其中 $X,Y,Z$ 是 $x,y,z$ 的 Walsh-Hadarmard 变换。

证明：采用归纳法，当 $k=0$ 时

$$Z(0)=z(0)=x(0)y(0)=X(0)Y(0)$$

假设当 $n=k$ 时成立，现在令 $n=k+1$，首先我们将 $x,y$ 拆分成两半 $X_i=H_{2^k}x_i,Y_i=H_{2^k}y_i$，有

$$X_i=H_{2^k}x_i=\begin{bmatrix}H_{2^k}&H_{2^k}\\H_{2^k}&-H_{2^k}\end{bmatrix}\begin{bmatrix}x_0\\x_1\end{bmatrix}$$ $$X=\begin{bmatrix}X_0+X_1\\X_0-X_1\end{bmatrix}\text{and}\,Y=\begin{bmatrix}Y_0+Y_1\\Y_0-Y_1\end{bmatrix}$$ $$z=x*y=\frac{1}{2}\begin{bmatrix}x_0*y_0+x_1*y_1\\x_0*y_1+x_1*y_0\end{bmatrix}$$

因此

$$\begin{align} Z &= H_{2^{k+1}}z\\ &= \begin{bmatrix}H_{2^k}&H_{2^k}\\H_{2^k}&-H_{2^k}\end{bmatrix}\begin{bmatrix}x_0*y_0+x_1*y_1\\x_0*y_1+x_1*y_0\end{bmatrix}\\ &= \begin{bmatrix}X_0Y_0+X_1Y_1+X_0Y_1+X_1Y_0\\X_0Y_0+X_1Y_1-X_0Y_1-X_1Y_0\end{bmatrix}\\ &= \begin{bmatrix}(X_0+X_1)(Y_0+Y_1)\\(X_0-X_1)(Y_0-Y_1)\end{bmatrix}\\ &= XY \end{align}$$

证毕。

下面证明 Walsh-Hadarmard 变换是满足这个卷积定理的唯一可逆线性变换（在忽略行顺序的意义下，因为行之间的交换不会影响这个性质）。

$$z(n)=\sum_{a\oplus b=n}x(a)y(b)$$

当这个变换是一阶时，显然，变换为单位变换。故下面设阶大于1。

引理1：

$$A_{:,i\oplus j}=A_{:,i}·A_{:,j}$$

其中 $A_{:,i}$ 表示 $A$ 的第 $i$ 列。

证明：首先，令 $x,y$ 为第 $i,j$ 位为 $1$，其他为 $0$，可以得到

$Az$ 为 $A$ 的 $i\oplus j$列。而 $Ax$ 为 $A$ 的第 $i$ 列，$Ay$ 为 $A$ 的第 $j$ 列。因此

$$A_{:,i\oplus j}=A_{:,i}·A_{:,j}$$

推论2：由引理1，令 $i=0$ 知道第一列全为 $1$，其他列都只有 $-1$ 和 $1$，对一个固定的 $t$，任意 $i$ 和 $i\oplus 2^t$ 两列的按位乘都相等。并且 $A$ 可逆所以所有列都不一样。

引理3（用处不大）： $A$ 的某一行全 $1$，其他行 $-1$ 和 $1$ 各占一半

证明：令 $y$ 为 $1$ 向量。则 $z=x*1=11^T$ 得

$$Az=Ax·A1=1^TxA1$$

因此除了某一行全为 $1$，其他行的和都为 $0$。这是因为，至多有一行为全 $1$，其他为 $0$，否则 $A$ 不可逆。另外，必然有一行为全 $1$，因为 $1$ 向量的补空间为 $2^k-1$ 维，如果没有一行全 $1$ ，则

$$A1=0\Rightarrow1=0$$

定理4：Walsh-Hadarmard 变换是满足这个卷积定理的唯一可逆线性变换（忽略行序的情况）

证明：由推论 1，要确定一行，只需确定其中的 $k$ 个位置的值，这 $k$ 个位置为 $2^0,2^1,\cdots,2^{k-1}$。由于推论 2，每个位置只能取 $\pm1$。因此，可能有 $2^k$ 种情况。

而 $A$ 共有 $2^k$ 行且可逆，因此每一种情况对应 $A$ 中的一行。故 $A$ 必然就是 Hadarmard 矩阵。（一开始想用归纳证明，把自己绕进去好久）

或卷积和与卷积

对于或以及与，沃尔什阿达玛变换就不满足条件了。然而有了上面证明唯一性的经验，我们就可以照猫画虎，来找满足条件的线性变换。

仿照引理 1，推论 2，可以得到，对于或卷积，这个变换 $A$ 必须满足 $$

$$A_{:,i|j}=A_{:,i}·A_{:,j}$$

并且第一列必须为 $1$。并且元素只能为 $0$ 或 $1$。令 $j=n-1=2^k-1$，得

$$A_{:,j}=A_{:,i}·A_{:,j},\forall i$$

所以最后一列，如果某个位置为 $1$，则该行为 $1$。因此，最后一列有且只有一个 $1$。其中 $2^k=n$。同理 $2^0,2^1,\cdots,2^{k-1}$ 这些列决定了整个矩阵。

同样 $A$ 只能有 $0/1$ 。所以所有行包含了 $2^k$ 种情况。我们可以让从下往上让第 $i$ 行这 $k$ 个位对应二进制 $i$。通过分析知道 $A$ 具有如下形式：

$$A=\begin{bmatrix}H&H\\H&0\end{bmatrix}$$

其中 $H$ 必须满足 $2^{k-1}$ 阶的卷积定理。因此满足或卷积定理的线性变可以通过如下方法构造（$H_1=[1]$）

$$H_{2^k}=\begin{bmatrix}H_{2^{k-1}}&H_{2^{k-1}}\\H_{2^{k-1}}&0\end{bmatrix}$$

或者（随便改变行序都行，不过为了变换方便，我们采用这两种方式之一）

$$H_{2^k}=\begin{bmatrix}H_{2^{k-1}}&0\\H_{2^{k-1}}&H_{2^{k-1}}\end{bmatrix}$$

同理可以得到满足与卷积定理的线性变换为可以通过下面方法构造

$$H_{2^k}=\begin{bmatrix}H_{2^{k-1}}&H_{2^{k-1}}\\0&H_{2^{k-1}}\end{bmatrix}$$

可以得到或形式快速变换可以写成

$$FWT_{or}(X)=\text{merge}(FWT_{or}(X_0),FWT_{or}(X_0)+FWT_{or}(X_1))$$ $$IFWT_{or}(X)=\text{merge}(IFWT_{or}(X_0),IFWT_{or}(X_0)-IFWT_{or}(X_1))$$

与形式的变换写成

$$FWT_{and}(X)=\text{merge}(FWT_{and}(X_0)+FWT_{and}(X_1),FWT_{and}(X_1))$$ $$IFWT_{and}(X)=\text{merge}(IFWT_{and}(X_0)-IFWT_{and}(X_1),IFWT_{and}(X_1))$$

C++实现

重新回顾一下三种快速变换

• OR

$$FWT_{or}(X)=\text{merge}(FWT_{or}(X_0),FWT_{or}(X_0)+FWT_{or}(X_1))$$ $$IFWT_{or}(X)=\text{merge}(IFWT_{or}(X_0),IFWT_{or}(X_0)-IFWT_{or}(X_1))$$

• AND

$$FWT_{and}(X)=\text{merge}(FWT_{and}(X_0)+FWT_{and}(X_1),FWT_{and}(X_1))$$ $$IFWT_{and}(X)=\text{merge}(IFWT_{and}(X_0)-IFWT_{and}(X_1),IFWT_{and}(X_1))$$

• XOR

$$FWT(X)=\text{merge}(FWT(X_0)+FWT(X_1),FWT(X_0)-FWT(X_1))\\ IFWT(X)=\text{merge}\left(\frac{IFWT(X_0)+IFWT(X_1)}{2},\frac{IFWT(X_0)-IFWT(X_1)}{2}\right)$$

实现起来比 FFT 简单得多，从底至上循环递推一下即可。可以把正负变换一起实现，因为只有系数不一样。

此处除了 Fwt 类，还实现了两个类，整数取模类ModInt, 多项式类Poly。为了复用性，牺牲了一点性能，一般不至于把时间卡那么死。这个多项式类和ModInt在其地方也都用到。还有个快速幂qpow也写成模板的形式。fwt部分的代码也就30行。

#include <bits/stdc++.h>
using namespace std;
constexpr int P = 998244353;

template<class T>
T qpow(T n, int k) {
    T res = 1;
    for (; k; k >>= 1, n *= n) {
    if (k & 1) {
        res *= n;
    }
    }
    return res;
}

template<int P = P>
struct ModInt {
    int x;
    ModInt(int x_ = 0) : x(mod(x_)) {}
    ModInt(const ModInt& that) : x(that.x) {}
    int val() const {
    return x;
    }

    // assume -P <= x < 2P
    int mod(int x) {
    if (x < 0) {
        x += P;
    }
    if (x >= P) {
        x -= P;
    }
    return x;
    }
    ModInt operator-() const {
    return ModInt(P - x);
    }
    ModInt inv() const {
    // assert(x != 0);
    return qpow(*this, P - 2);
    }
    ModInt& operator*=(const  ModInt& rhs) {
    x = int64_t(x) * rhs.x % P;
    return *this;
    }
    ModInt& operator+=(const  ModInt& rhs) {
    x = mod(x + rhs.x);
    return *this;
    }
    ModInt& operator-=(const  ModInt& rhs) {
    x = mod(x - rhs.x);
    return *this;
    }
    ModInt& operator/=(const  ModInt& rhs) {
    return *this *= rhs.inv();
    }
    friend  ModInt operator*(const  ModInt& lhs, const  ModInt& rhs) {
    ModInt res = lhs;
    res *= rhs;
    return res;
    }
    friend  ModInt operator+(const  ModInt& lhs, const  ModInt& rhs) {
    ModInt res = lhs;
    res += rhs;
    return res;
    }
    friend  ModInt operator-(const  ModInt& lhs, const  ModInt& rhs) {
    ModInt res = lhs;
    res -= rhs;
    return res;
    }
    friend  ModInt operator/(const  ModInt& lhs, const  ModInt& rhs) {
    ModInt res = lhs;
    res /= rhs;
    return res;
    }
};

template<typename T>
struct Poly {
    vector<T> v;
    Poly(int n = 0) {
    v.resize(n);
    }

    Poly(const Poly& p) {
    v.resize(p.size());
    for (int i = 0; i < size(); i++) {
        v[i] = p[i];
    }
    }

    size_t size() const {
    return v.size();
    }

    T& operator[](int i) {
    return v[i];
    }

    T operator[](int i) const {
    return v[i];
    }

    Poly operator-() const {
    Poly p = v;
    for (int i = 0; i < p.size(); i++) {
        p[i] = -p[i];
    }
    return p;
    }

    Poly& operator*=(Poly& rhs) {
    for (int i = 0; i < size(); i++) {
        v[i] *= rhs[i];
    }
    return *this;
    }
    Poly& operator+=(const Poly& rhs) {
    for (int i = 0; i < size(); i++) {
        v[i] += rhs[i];
    }
    return *this;
    }
    Poly& operator-=(const Poly& rhs) {
    for (int i = 0; i < size(); i++) {
        v[i] -= rhs[i];
    }
    return *this;
    }
    Poly& operator/=(const Poly& rhs) {
    for (int i = 0; i < size(); i++) {
        v[i] /= rhs[i];
    }
    return *this;
    }
    friend Poly operator*(Poly lhs, Poly rhs) {
    return lhs *= rhs;
    }
    friend Poly operator+(Poly lhs, Poly rhs) {
    return lhs += rhs;
    }
    friend Poly operator-(Poly lhs, const Poly rhs) {
    return lhs -= rhs;
    }
    friend Poly operator/(Poly lhs, Poly rhs) {
    return lhs /= rhs;
    }
};

struct Fwt {
    template<typename T>
    static Poly<T> OR(Poly<T> f, bool inv = false) {
    T x = inv ? -1 : 1;
    int n = f.size();
    for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for (int i = 0; i < n; i += o)
        for (int j = 0; j < k; j++)
            f[i + j + k] += f[i + j] * x;
    return f;
    }
    template<typename T>
    static Poly<T> AND(Poly<T> f, bool inv = false) {
    T x = inv ? -1 : 1;
    int n = f.size();
    for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for (int i = 0; i < n; i += o)
        for (int j = 0; j < k; j++)
            f[i + j] += f[i + j + k] * x;
    return f;
    }

    template<typename T>
    static Poly<T> XOR(Poly<T> f, bool inv = false) {
    T x = inv ? T(1) / 2 : 1;
    int n = f.size();
    for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for (int i = 0; i < n; i += o)
        for (int j = 0; j < k; j++)
            f[i + j] += f[i + j + k],
            f[i + j + k] = f[i + j] - f[i + j + k] - f[i + j + k],
            f[i + j] *= x, f[i + j + k] *= x;
    return f;
    }

    template<typename T>
    static Poly<T> IOR(Poly<T> f) {
    return Fwt::OR(f, true);
    }
    template<typename T>
    static Poly<T> IAND(Poly<T> f) {
    return Fwt::AND(f, true);
    }
    template<typename T>
    static Poly<T> IXOR(Poly<T> f) {
    return Fwt::XOR(f, true);
    }

};
template<int P>
istream& operator>>(istream& ins,
    ModInt<P>& v) {
    ins >> v.x;
    return ins;
}

template<typename T>
istream& operator>>(istream& ins,
    vector<T>& v) {
    for (int i = 0; i < int(v.size()); i++) ins >> v[i];
    return ins;
}

template <typename T>
istream& operator>>(istream& ins,
    Poly<T>& p) {
    ins >> p.v;
    return ins;
}

template<int P>
ostream& operator<<(ostream& outs,
    const ModInt<P>& v) {
    outs << v.val();
    return outs;
}

template<typename T>
ostream& operator<<(ostream& outs,
    const vector<T>& v) {
    for (int i = 0; i < int(v.size()); i++) outs << v[i] << " ";
    return outs;
}

template <typename T>
ostream& operator<<(ostream& outs,
    Poly<T>& p) {
    outs << p.v;
    return outs;
}

int main() {
    int n, m;
    cin >> m;
    n = 1 << m;

    Poly<ModInt<>> a(n), b(n);
    cin >> a >> b;
    Poly<ModInt<>> c = Fwt::IOR(Fwt::OR(a) * Fwt::OR(b));
    Poly<ModInt<>> d = Fwt::IAND(Fwt::AND(a) * Fwt::AND(b));
    Poly<ModInt<>> e = Fwt::IXOR(Fwt::XOR(a) * Fwt::XOR(b));

    cout << c << "\n" << d << "\n" << e << "\n";

    return 0;
}