Fast Walsh–Hadamard Transform

Introduction

The Walsh Transform¹ is a method used in spectral analysis as an alternative to the Discrete Fourier Transform, widely applied in signal processing. FFT uses double precision, but Walsh decomposes the signal under square waves of different oscillation frequencies, so all coefficients are integers with the same absolute magnitude. This eliminates the need for floating-point multiplication, improving computational speed.

Thus, the core idea of FWT and FFT should be the same: both are transformations of arrays. We denote the result of applying the Fast Walsh Transform to array \(A\) as \(FWT[A]\).

The core idea of FWT is:

We need a new sequence \(C\), obtained from sequences \(A\) and \(B\) through some operation rule, i.e., \(C = A \cdot B\);

We first compute \(FWT[A]\) and \(FWT[B]\) in the forward direction, then compute \(FWT[C]=FWT[A] \cdot FWT[B]\) in \(O(n)\) time complexity, where \(\cdot\) denotes element-wise multiplication of sequences;

Then we perform the inverse operation to obtain the original sequence \(C\). The time complexity is \(O(n \log{n})\).

In competitive programming, FWT is used to solve convolution problems where indices undergo bitwise operations.

Formula: \(C_{i} = \sum_{i=j \oplus k}A_{j} B_{k}\)

(where \(\oplus\) is one of the binary bitwise operations)

Below we take \(\cup\) (bitwise OR), \(\cap\) (bitwise AND), and \(\oplus\) (bitwise XOR) as examples.

Operations of FWT

OR Operation

If \(k=i\cup j\), then the positions where \(i\)'s binary bits are \(1\) and where \(j\)'s binary bits are \(1\) must be subsets of the positions where \(k\)'s binary bits are \(1\).

Now to obtain \(FWT[C] = FWT[A] \cdot FWT[B]\), we need to construct the FWT rule.

By definition, we can construct \(FWT[A]_i = A'_i = \sum_{i=i\cup j}A_{j}\), meaning \(j\) satisfies that its binary \(1\) bits are a subset of \(i\)'s.

Then we have:

\[ \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} \]

Next, let's see how to compute \(FWT[A]\).

Obviously, we cannot enumerate all, as the complexity would be \(O(n^2)\). Since we cannot enumerate entirely, we consider divide and conquer.

We split the entire interval in half. After splitting, writing the indices in binary reveals a pattern.

Let \(A_0\) denote the first half of \(A\), and \(A_1\) denote the second half. Then \(A_0\) corresponds to indices where the highest bit is \(0\), and its subsets are exactly its own subsets (since the highest bit is \(0\)). However, \(A_1\) has highest bit \(1\), and its qualifying subsets include not only itself but also subsets with highest bit \(0\), i.e.,

\[ FWT[A] = merge(FWT[A_0], FWT[A_0] + FWT[A_1]) \]

Here merge means concatenating two arrays like strings, and \(+\) is ordinary addition, meaning adding corresponding binary bits.

Thus, we can complete the concatenation in \(O(\log{n})\) time via binary splitting, with each concatenation requiring one operation, meaning we obtain \(FWT[A]\) in \(O(n\log{n})\) time complexity.

Next is inversion. Actually, inversion is simple. Since we know \(A_0\)'s subsets are itself (\(A_0 = FWT[A_0]\)), and \(A_1\)'s subsets are \(FWT[A_0] + FWT[A_1]\), we easily derive the inversion recurrence:

\[ UFWT[A'] = merge(UFWT[A_0'], UFWT[A_1'] - UFWT[A_0']) \]

Below we provide code implementation. Note that forward and inverse transforms can be merged into one function, with \(\text{type}=1\) for forward transform and \(\text{type}=-1\) for inverse transform.

Implementation

void Or(ll *a, ll type) {  // Iterative implementation, smaller constant
  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;
      }
    }
  }
}

AND Operation

Analogous to OR, we can derive similar conclusions:

\[ FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_1]) \]

\[ UFWT[A'] = merge(UFWT[A_0'] - UFWT[A_1'], UFWT[A_1']) \]

Below we provide code implementation. \(\text{type}=1\) for forward transform, \(\text{type}=-1\) for inverse transform.

Implementation

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;
      }
    }
  }
}

XOR Operation

The convolution for XOR is based on the following principle:

If we let \(x\circ y\) denote the parity of the number of \(1\)s in \(x\cap y\), i.e., \(x\circ y=\text{popcnt}(x\cap y)\bmod 2\), then we easily have \((x\circ y)\oplus (x\circ z)=x\circ(y\oplus z)\).

The operation for \(FWT[A]\) is also easy to derive.

Let \(FWT[A]_i=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j\). We verify the correctness of \(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} \]

Now let's see how to compute \(A,B\) quickly, still using divide and conquer:

For the subsequence \(FWT[A_0]\) where \(i\)'s current bit is \(0\), when performing \(\circ\) operation, we find that computing with \(0\) or \(1\) does not change the result (since \(0\cap 0=0,0\cap1=0\)), so \(\sum_{i\circ j=1}A_j=0\) in \(FWT[A]=\sum_{i\circ j=0}A_j-\sum_{i\circ j=1}A_j\).

For the subsequence \(A_1\) where \(i\)'s current bit is \(1\), when performing \(\circ\) operation, computing with \(0\) yields \(0\), and with \(1\) yields \(1\) (since \(1\cap 0=0,1\cap1=1\)).

In summary:

\[ FWT[A]=merge((FWT[A_0]+FWT[A_1])-0, FWT[A_0]-FWT[A_1]) \]

That is:

\[ FWT[A] = merge(FWT[A_0] + FWT[A_1], FWT[A_0] - FWT[A_1]) \]

Inverse transform is easily obtained:

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

Code is provided. \(\text{type}=1\) for forward transform, \(\text{type}=\frac{1}{2}\) for inverse transform.

Implementation

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;
      }
    }
  }
}

XNOR Operation

Analogous to XOR, we give the formulas:

\(FWT[A]_{i} = \sum_{C_1}A_{j} - \sum_{C_2}A_{j}\) (where \(C_1\) means \(\text{popcnt}(x\cup y)\bmod 2\) is \(0\), \(C_2\) means \(\text{popcnt}(x\cup y)\bmod 2\) is \(1\))

\[ FWT[A] = merge(FWT[A_1] - FWT[A_0], FWT[A_1] + FWT[A_0]) \]

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

Another Perspective on FWT

Let \(c(i,j)\) be the contribution coefficient of \(A_j\) to \(FWT[A]_i\). We can redescribe the FWT transformation process:

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

Since:

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

We can easily prove: \(c(i,j)c(i,k)=c(i,j\odot k)\). Here \(\odot\) is any bitwise operation.

Additionally, the \(c\) function has an important property: it can be processed bitwise.

For example, during transformation:

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

This is inefficient. We split it:

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

Considering the difference between the first and second parts for \(i,j\), we find only the highest bit differs.

So we remove the highest bit of \(i,j\) to get \(i',j'\), and denote the highest bit of \(i\) as \(i_0\). Then:

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

If \(i_0=0\), then:

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

If \(i_0=1\), then:

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

That is, we only need:

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

These four numbers to complete the transformation. We call this matrix the bit matrix.

If we want the inverse transform, we need the inverse of the above bit matrix.

If the inverse matrix is \(c^{-1}\), we can recover the original numbers via a similar operation:

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

The inverse matrix may not exist; for example, if there is a row or column of \(0\)s, the matrix is not invertible. We must be careful when constructing.

Bitwise OR

We can construct:

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

This satisfies \(c(i,j)c(i,k)=c(i,j\cup k)\). We find this matches exactly what we derived earlier: \(FWT[A]=\text{merge}(FWT[A_0], FWT[A_0]+FWT[A_1])\)! Similarly, the following matrix also satisfies \(c(i,j)c(i,k)=c(i,j\cup k)\), but we generally use the one above:

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

Although the following matrix also satisfies \(c(i,j)c(i,k)=c(i,j\cup k)\), it has a row of \(0\)s, so it is not invertible and thus invalid:

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

For inverse transformation, we need to invert the matrix. Taking the topmost matrix as an example, we get:

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

Then, following the forward transform method, substitute the inverse transform matrix.

Bitwise AND

We can construct:

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

This satisfies \(c(i,j)c(i,k)=c(i,j\cap k)\).

Inverse matrix:

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

Bitwise XOR

We can construct:

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

This satisfies \(c(i,j)c(i,k)=c(i,j\oplus k)\).

Inverse matrix:

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

FWT is a Linear Transformation

FWT is a linear transformation. That is, it satisfies:

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

and:

\[ FWT[c\cdot A]=c\cdot FWT[A] \]

K-dimensional FWT

Actually, bitwise operations essentially operate on an \(n\)-dimensional \(\{0,1\}\) vector. OR is taking \(\max\) in each dimension. AND is taking \(\min\) in each dimension. XOR is adding corresponding dimensions modulo \(2\).

Bitwise operations have a characteristic: each bit of the vector is independent.

We extend \(\{0,1\}\) to \([0,K)\cap \mathbf{Z}\), i.e., to base \(K\), and see what we get?

max Operation

We extend \(\cup\) to base \(K\), defining \(i\cup j\) as taking \(\max\) bitwise, with:

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

If \(j=k\), then the above becomes:

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

That is, the \(1\)s in each row must appear before the \(0\)s; if they appear after, it's invalid. Playing around yields a valid construction:

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

Inverting gives:

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

min Operation

We extend \(\cap\) to base \(K\), defining \(i\cap j\) as taking \(\min\) bitwise, with:

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

If \(j=k\), then the above becomes:

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

That is, the \(1\)s in each row must appear after the \(0\)s; if they appear before, it's invalid. Playing around yields a valid construction:

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

Inverting gives:

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

The first two are less commonly used; the following is more common:

Non-carry Addition

We extend \(\oplus\) to base \(K\), defining \(i\oplus j\) as adding bitwise then \(\bmod K\), with:

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

We construct \(c(i,j)=\omega_{K}^j\), which satisfies the requirement:

\[ \omega_{K}^j\omega_{K}^k=\omega_{K}^{j\oplus k} \]

But if every row is the same, the matrix is not invertible. So we can construct \(c(i,j)=\omega_{K}^{(i-1)j}\).

We have the following matrix:

\[ \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} \]

This is the Vandermonde matrix. Its inverse is:

\[ \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} \]

If the modulus given in the problem has primitive roots of unity, we can implement this simply.

But primitive roots of unity may not exist modulo the given modulus, so we consider field extension: artificially define an \(x\) satisfying \(x^K=1\), then directly substitute \(x\) into computations. Each number becomes a polynomial of degree \(k-1\) in \(x\). We only need to compute modulo \(x^K-1\). Then the matrix can be expressed as:

\[ \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} \]

But this may introduce zero divisors, meaning a number may have multiple representations, and we cannot determine a number's true value.

We consider not modulo \(x^K-1\), but modulo the cyclotomic polynomial \(\Phi_{K}(x)\), which satisfies that \(x\) has order \(k\) and is irreducible over \(\mathbb{Q}\). So we define the above computation to be modulo \(\Phi_{K}(x)\).

Another issue is that modulo \(\Phi_{K}(x)\) has large constants (since \(\Phi\) itself is a polynomial). But because \(\Phi_{K}(x)\mid x^k-1\), we can compute modulo \(x^k -1\) and finally reduce modulo \(\Phi_{K}(x)\).

Example Problems

「CF 1103E」Radix sum

Given a sequence \(a_1,a_2,...,a_n\) of length \(n\), for each \(p \in [0,n-1]\), find the number of integer sequences \(i_1,i_2,...,i_n\) satisfying the following conditions, modulo \(2^{58}\):

• \(\forall j \in [1,n] , i_j \in [1,n]\);
• \(\sum\limits_{j=1}^n a_{i_j} = p\), where addition is defined as decimal non-carry addition.

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

Solution

We can think of dp: design state \(f_{i,s}\) meaning considering the \(i\)-th number, current addition state is \(s\). Since FWT transformation is linear, we can first transform to FWT point-value representation, then raise to the \(n\)-th power, and finally transform back.

The above is straightforward, but the problem gives modulus \(2^{58}\). We find no primitive roots of unity, so consider field extension.

Here the cyclotomic polynomial is \(\Phi_{10}(x)=x^4-x^3+x^2-x+1\).

However, during UFWT, we need to divide by base \(10\), but we find \(10\) has no inverse modulo \(2^{58}\). Actually, we find \(5\) does have an inverse modulo \(2^{58}\): \(57646075230342349\). We just need to divide by \(2\) additionally. Let the answer after dividing by \(5\) be \(x\), the true answer be \(y\), i.e., \(2^5y\equiv x\pmod{2^{64}}\). Clearly, we have \(y\equiv \frac{x}{2^5}\pmod{2^{64-5}}\), i.e., \(y\equiv \frac{x}{2^5}\pmod{2^{59}}\), so we can directly divide the final answer by \(2^5\). Although the problem setter didn't know why modulo \(2^{58}\) was chosen, we can take modulo again.

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

Given an ellipse \(E\), where all integer point coordinates are within \([1,4 \cdot 10^6]\). Find \(\sum_{(x,y) \in E} (x \oplus y)^{33}x^{-2}y^{-1} \mod 10^9+7\).

Solution

This is a less straightforward problem. The problem setter provides detailed English solution; see this link.

References

• Peach Sauce's Algorithm Notes
• ZnPdCo's Blog

---

1. Wikipedia ↩

---

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用