# Fast Walsh-Hadamard Transform in Competitive Programming

This should be the very first English post I write on my blog and I expect there to be some minor errors. This is a popular technique in the Chinese competitive programming community but there doesn't seem to be a lot of documentation about its application in the English CP community. The posts I found on Codeforces doesn't seem to be very clear to me...

# Prerequisites

- A decent proficiency in competitive programming.
- A basic understanding of the Cooley-Tukey FFT and its application in competitive programming.
- A decent understanding of bitwise operations.

# Why do we need FWHT?

Recall what we would do if we are to quickly calculate the following convolution of two sequences and , each of length :

We use FFT, which applies the following transformation to the input sequence: Since the calculation of this transformation (and its inverse) can be done in a divide-and-conquer manner in and the element wise product of the transformation is equivalent to the convolution on the original series, we are able to calculate the convolution in .Now we try to generalize our findings to a more general case:

where is some binary operation. The convolution we see at the beginning is a special case where .FWHT is an algorithm that borrows similar notions from FFT and is able to compute the convolution in time for (bitwise OR, bitwise AND, and bitwise XOR). Why do the convolutions of these bitwise operations matter? Observe that binary representation is a way of encoding sets and these three operations correspond to set union, set intersection and set symmetric difference respectively, therefore, FWHT can be used to accelerate set-based DPs.

# Bitwise OR convolution

Let's start with the convolution with respect to bitwise OR:

We start by exploiting an interesting property of bitwise OR: or its clearer equivalent in set-based notations:**Claim:**The following transformation can turn OR convolutions to element-wise multiplications:

**Proof:**Then how are we able to compute quickly? A trivial implementation still takes time.

Recall what we did in FFT: we divide
into two subsequences based on parity of indices, a.k.a, the last bit of
indices. We did this because the root of unity has such amazing property
as .
We could do that here as well, but a limitation of dividing based on the
last bit is that the order of elements changes in the process, so an
efficient in-place implementation has to do a pre-shuffle to cope with
that. Since OR is a bitwise operation, which bit based on which we
divide doesn't really matter much. Why not simply **divide based
on the first, or the most significant bit**, such that the order
of elements is preserved in the process? Dividing based on the highest
bit of indices, simply put, is to split
into the first half,
,
and the second half,
,
in their natural order.

Here I introduce a notation, or . In the context where the length of the sequence is (and is a power of ), where , and is just . In other words, has as the highest bit and has as the highest bit.

(Note using this notation, and )

To make our writing clearer, denote

We want to express each element of as a combination of some element in and .We first look at the first half of . Using the notation I defined above, these elements can be expressed as .

We know that the highest bit of should always be , so the condition in the second summation is never satisfied, and we can simply throw the second term away. And simplifies to , so we get, by definition of : What about the second half, ? Together we get: with the trivial recursion boundary when .(here I use the tuple notation to denote concatenation, and to denote element-wise addition).

This is something we can write an in-place implementation for with ease:

```
void fwht_or(int n, int *a, int *A) {
copy(a, a + n, A);
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1)
for (int l = 0; l < n; l += s)
for (int i = 0; i < h; i++)
A[l + h + i] += A[l + i]
}
```

Its time complexity is obviously with a really small constant factor.

Its reverse transform turns out to be simple as well, suppose we know and let

(Assuming is a power of and and each have length )Then we can recover and :

Implementation:```
void ifwht_or(int n, int *a, int *A) {
std::copy(A, A + n, a);
// If n is guaranteed to be a power of 2 then we don't need n_ and the min(...) in the inner loop.
int n_ = 1; while (n_ < n) n_ <<= 1;
// n_ = 1 << (32 - __builtin_clz(n - 1));
for (int s = n_, h = n_ / 2; h; s >>= 1, h >>= 1)
for (int l = 0; l < n; l += s)
for (int i = 0; i < std::min(i, n - l - h); i++)
a[l + h + i] -= a[l + i]
}
```

And an amazing thing about this, which I haven't quite figured out why, is that the order of the outermost loop above can be reversed and both functions can be merged into one:

```
void fwht_or(int n, int *a, int *A, int dir = 1) {
std::copy(a, a + n, A);
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1)
for (int l = 0; l < n; l += s)
for (int i = 0; i < h; i++)
A[l + h + i] += dir * A[l + i]
}
```

(Fast bitwise OR / set union convolution is sometimes aliased "Fast Mobius Transform" in Chinese CP community. Both are essentially the same.)

# Bitwise AND convolution

The bitwise AND convolution

can be accelerated in a similar way.(Actually, by de Morgan's Law we can always reduce an AND convolution to an OR convolution)

Note that AND also has this interesting property:

or in set notations: Thus, we can prove in a way similar to what we did in OR convolution that the transform can turn convolutions to element-wise multiplications.We still adopt the same divide-and-conquer approach and continue to use the notations .

Consider the first half of , which can be expressed as :

And by the properties of AND, both and simplify to . So by definition we get Then consider the other half of : Together we have: with the trivial recursion boundary when .This gives an efficient implementation very similar to
`fwht_or`

above:

```
void fwht_and(int n, int *a, int *A) {
std::copy(a, a + n, A);
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1)
for (int l = 0; l < n; l += s)
for (int i = 0; i < h; i++)
A[l + i] += A[l + h + i]
}
```

The inverse transform is simple as well. Let , then

The code:```
void ifwht_and(int n, int *a, int *A) {
std::copy(A, A + n, a);
// If n is guaranteed to be a power of 2 then we don't need n_ and the min(...) in the inner loop.
int n_ = 1; while (n_ < n) n_ <<= 1;
for (int s = n_, h = n_ / 2; h; s >>= 1, h >>= 1)
for (int l = 0; l < n; l += s)
for (int i = 0; i < std::min(i, n - l - h); i++)
a[l + i] -= a[l + h + i]
}
```

The order of the outermost loop can be reversed and we can also merge the functions above together.

# Bitwise XOR convolution

The XOR operation does **not** have such nice property
as

We first introduce an auxiliary operation, define , where denotes the number of s in the binary representation of .

**Claim:** The transformation below turns convolutions
to element-wise multiplications:

**Proof:**How to simplify those terms?

Observe that by the definition of XOR we have

So if we apply modulo on both sides, Plug in and we get We are almost there. Apply the identity below, whose proof I simply omit here, (This is something good about bitwise operations: if you cannot prove an identity the smart way you can always fall back on the dumb method -- making a truth table)We finally get

(We are actually quite familiar with this if we remove the circles outside s and s)With this conclusion we can simplify the four terms above:

which completes the proof.We then explore how to compute efficiently. Divide and conquer is still our friend, and dividing based on the highest bit works here so we continue to use those notations.

Consider the first half of ...

and the other half: So together we get and the inverse transform The code for both transforms are a bit longer than those for OR and AND, but not by too much:```
void fwht_xor(int n, int *a, int *A) {
std::copy(a, a + n, A);
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1) {
for (int l = 0; l < n; l += s) {
for (int i = 0; i < h; i++) {
int t = A[l + h + i];
A[l + h + i] = A[l + i] - t;
A[l + i] += t;
}
}
}
}
void ifwht_xor(int n, int *a, int *A) {
std::copy(A, A + n, a);
// If n is guaranteed to be a power of 2 then we don't need n_ and the min(...) in the inner loop.
int n_ = 1; while (n_ < n) n_ <<= 1;
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1) {
for (int l = 0; l < n; l += s) {
for (int i = 0; i < std::min(i, n - l - h); i++) {
int t = a[l + h + i];
a[l + h + i] = (a[l + i] - t) / 2;
a[l + i] = (a[l + i] + t) / 2;
}
}
}
}
```

They can be merged as well:

```
void fwht_xor(int n, int *a, int *A, bool inv = false) {
std::copy(a, a + n, A);
for (int s = 2, h = 1; s <= n; s <<= 1, h <<= 1) {
for (int l = 0; l < n; l += s) {
for (int i = 0; i < h; i++) {
int t = A[l + h + i];
A[l + h + i] = A[l + i] - t;
A[l + i] += t;
if (inv) A[l + h + i] /= 2, A[l + i] /= 2;
}
}
}
}
```

**This code above is what Wikipedia refers to as the authentic
Fast Walsh-Hadamard Transform**.

# Some sidenotes

Note that though FFT and FWHT shares the same idea of divide and conquer, FWHT does not require to be a power of whereas FFT does. (Well actually neither of them "require" to be a power of , but to apply FFT when is not a power of you either need to pad with s or you have to make your implementation really complicated).

Also, I just came to know that if we express WHT in the language of matrices and vectors, the matrix is called a Hadamard Matrix.

Another fact that I didn't quite understand is why the order of the inverse FWHT can be reversed.

For instance, when , after fully dividing the sequence into individual elements, we first merge , then we merge and finally . Naturally when we do the inverse transform we have to start with , recover , then recover and then recover the individual elements. But the popular implementation seems to suggest that the inverse transformation algorithm works in another direction as well. I am now puzzled why this is true and currently I'm just taking this for granted. Perhaps I derived the inversion in a different way than others did? If you have a simple explanation please leave a comment :)