FWHT (Fast Walsh-Hadamard Transform)

다음 글을 정리한 것이다. https://alan20210202.github.io/2021/07/25/FWHT/

 

문제

binary operator $\star$과 수열 $a, b$에 대해 연산 $\circledast$를 아래와 같이 정의하자

$$\{a \circledast b \}_i = \sum_{j \star k = i} {a_j b_k} $$

$\star$가 $+$일 때 이 convolution은 FFT를 이용해 해결할 수 있다. 

이제 $\star$가 bitwise or, and, xor인 경우를 해결해보자. 편의상 $a$와 $b$의 길이 $n$은 2의 거듭제곱이라 하자.

 

or convolution은 Fast Mobius Transform이라고도 부르고, xor convolution만을 FWHT라 부르기도 한다.

표기

우선 풀이에 들어가기 앞서, 몇 가지 표기를 정의하자.

$0|x$와 $1|x$는 각각 $x$라는 수의 앞에 bit 0과 1을 추가한 것이다. 즉, $0|x=x$이고 $1|x=x+n/2$이다.

수열 $a$에 대해 $a^0$은 $a$의 앞쪽 절반이고 $a^1$은 $a$의 뒤쪽 절반이다. $a^0_i=a_i=a_{0|i}$이고, $a^1_i=a_{i+n/2}=a_{1|i}$이다.

Bitwise or

아래와 같은 문제를 풀고 싶은 상황이다.

$$\{a \circledast b \}_i = \sum_{j \vee k = i} {a_j b_k} $$

우선 아래와 같은 연산을 정의하자.

$$ \mathcal{FWHT}\{a\}_i=\sum_{j\vee i=i}\{a_j\}$$

그러면 아래와 같은 성질이 성립한다.

$$ \begin{split} \mathcal{FWHT}\{a\}_i \cdot \mathcal{FWHT}\{b\}_i & = \left( \sum_{j\vee i=i}{a_j}\right) \left( \sum_{k\vee i=i}{b_k} \right) \\ & = \sum_{j\vee i=i} \sum_{k\vee i=i} a_j b_k  \\ & = \sum_{(j\vee k)\vee i=i} a_j b_k \\ & = \sum_{l\vee i=i} \sum_{j\vee k=l} a_j b_k \\ & = \mathcal{FWHT}\{ a \circledast b \} _i \end{split}$$

or convolution이 element-wise multiplication으로 바뀐 것이다. 따라서 $\mathcal{FWHT}$와 그 역연산을 빠르게($O(NlogN)$) 구하면 된다.

이 $\mathcal{FWHT}$는 SOS dp와 같다. 그 계산 방법을 다시 한 번 살펴보자.

 

$A^0 = \mathcal{FWHT}\{a^0\}$과 $A^1 = \mathcal{FWHT}\{a^1\}$을 이미 구한 상황에서 $A=\mathcal{FWHT}\{a\}$를 구하고자 한다.

$\sum_{j \vee (0|i) = 0|i}{a_j} = \sum_{(0|j) \vee (0|i) = 0|i}{a_{0|j}}$이므로 $A_{0|i} = A^0_{i}$이다.

$\sum_{j \vee (1|i) = 1|i}{a_j} = \sum_{(0|j) \vee (0|i) = 0|i}{a_{0|j}} + \sum_{(1|j) \vee (1|i) = 1|i}{a_{1|j}}$이므로 $A_{1|i} = A^0_{i}+A^1_{i}$이다. 따라서 아래와 같이 쓸 수 있다.

$$ A = (A^0, A^0 + A^1)$$

이를 구현한 코드는 아래와 같다. (SOS dp에서의 코드와 사실상 같다.)

void fwht_or(vector<ll> &a) {
    ll n = a.size();
    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];
}

 

$A=(A', A'')$일 때 $A^0$와 $A^1$을 구하는 역연산은 아래와 같다.

$$\begin{split}  A^0 & = A' \\ A^1 & = A'' - A' \end{split}$$

이를 구현한 코드는 아래와 같다.

void ifwht_or(vector<ll> &a) { 
    ll n = a.size();
    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 < min(h, n - l - h); i++)
                a[l + h + i] -= a[l + i];
}

그리고 이 모든 연산이 bit의 순서와는 관련 없다는 성질을 이용하면, 위 코드에서 s를 n부터 줄이는 것이 아니라 1부터 늘려도 된다. 그러면 $\mathcal{FWHT}$와 그 역연산을 하나의 함수로 처리할 수 있다.

void fwht_or(vector<ll> &a, bool inv) {
    ll n = a.size();
    int dir = inv ? -1 : 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 < h; i++)
                a[l + h + i] += dir * a[l + i];
}

Bitwise and

$$\{a \circledast b \}_i = \sum_{j \wedge k = i} {a_j b_k} $$

and convolution도 or convolution과 비슷하게 해결할 수 있다.

$$\mathcal{FWHT} \{ a \}_i = \sum_{j \wedge i = i}{a_i}$$

위와 같이 놓으면 or에서와 마찬가지로 and convolution이 element-wise multiplicatoin으로 바뀐다.

$$\mathcal{FWHT} \{ a \circledast b \} _i = \mathcal{FWHT} \{ a \} _i \cdot \mathcal{FWHT} \{ b \} _i$$

$\mathcal{FWHT}$는 아래와 같은 성질을 만족한다.

$$A = (A^0 + A^1, A^1)$$

구현은 아래와 같다.

void fwht_and(vector<ll> &a, bool inv) {
    ll n = a.size();
    int dir = inv ? -1 : 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 < h; i++)
                a[l + h] += dir * a[l + h + i];
}

Bitwise xor

$$\{a \circledast b \}_i = \sum_{j \oplus k = i} {a_j b_k} $$

xor convolution은 조금 더 복잡하다.

우선, $\textrm{popcount}(x)$가 $x$에서 1인 bit의 개수를 의미한다 하자. 그리고 아래와 같은 연산을 정의하자.

$$ x \otimes y = \textrm{popcount}(x \wedge y) \: \textrm{mod} \: 2$$

$\mathcal{FWHT}$를 아래와 같이 정의하자.

$$\mathcal{FWHT} \{ a \} _i = \sum_{j \otimes i = 0} {a_j} - \sum_{j \otimes i = 1} {a_j}$$

그러면 아래와 같이 정리할 수 있다.

$$\mathcal{FWHT} \{ a \} _i \cdot \mathcal{FWHT} \{ b \}_i  = \left( \sum_{j \otimes i = 0} {a_j} - \sum_{j \otimes i = 1} {a_j} \right)\left( \sum_{k \otimes i = 0} {b_k} - \sum_{k \otimes i = 1} {b_k} \right) \\  = \sum_{j \otimes i = 0} \sum_{k \otimes i = 0} a_j b_k + \sum_{j \otimes i = 1} \sum_{k \otimes i = 1} a_j b_k - \sum_{j \otimes i = 1} \sum_{k \otimes i = 0} a_j b_k - \sum_{j \otimes i = 0} \sum_{k \otimes i = 1} a_j b_k$$

위 식을 정리하기 위해 $\oplus$와 $\otimes$의 성질을 이용해야 한다.

$$\begin{split} \textrm{popcount}(x \oplus y) & \equiv \textrm{popcount}(x) + \textrm{popcount}(y) - 2\textrm{popcount}(x \wedge y) \: & (\textrm{mod} \: 2) \\ & \equiv \textrm{popcount} (x) + \textrm{popcount}(y) \: & (\textrm{mod} \: 2) \end{split}$$

$$\textrm{popcount}((j \wedge i) \oplus (k \wedge i)) \equiv \textrm{popcount}(j \wedge i) + \textrm{popcount}(k \wedge i) \: (\textrm{mod} \: 2)$$

그리고 $ (j \wedge i) \oplus (k \wedge i) = (j \oplus k) \wedge i$이므로 아래와 같다.

$$\textrm{popcount}((j \oplus k) \wedge i) \equiv \textrm{popcount}(j \wedge i) + \textrm{popcount}(k \wedge i) \: (\textrm{mod} \: 2)$$

따라서 아래와 같다.

$$(j \oplus k) \otimes i \equiv j \otimes i + k \otimes i \: (\textrm{mod} \: 2)$$

이를 $\mathcal{FWHT}$ 식에 적용해주면

$$\begin{split} \mathcal{FWHT} \{ a \} _i \cdot \mathcal{FWHT} \{ b \} _i & = \sum_{(j \oplus k) \otimes i = 0} a_jb_k - \sum_{(j \oplus k) \otimes i = 1} a_j b_k \\ & = \mathcal{FWHT} \{ a \circledast b \}_i \end{split}$$

xor convolution이 element-wise multiplication으로 바뀐다.

 

이제 $\mathcal{FWHT}$를 빠르게 구하자.

$$\begin{split} A_{0|i} & = \sum_{j \otimes (0|i) = 0} a_j - \sum_{j \otimes (0|i) = 1} a_j \\ & = \sum_{(0|j) \otimes (0|i) = 0} a^0_j + \sum_{(1|j) \otimes (0|i) = 0} a^1_j - \sum_{(0|j) \otimes (0|i) = 1} a^0_j - \sum_{(1|j) \otimes (0|i) = 1} a^1_j \\ & = \sum_{j \otimes i = 0} a^0_j + \sum_{j \otimes i = 0} a^1_j - \sum_{j \otimes i = 1} a^0_j - \sum_{j \otimes i = 1} a^1_j \\ & = A^0_i + A^1_i \end{split}$$

$$\begin{split} A_{1|i} & = \sum_{j \otimes (1|i) = 0} a_j - \sum_{j \otimes (1|i) = 1} a_j \\ & = \sum_{(0|j) \otimes (1|i) = 0} a^0_j + \sum_{(1|j) \otimes (1|i) = 0} a^1_j - \sum_{(0|j) \otimes (1|i) = 1} a^0_j - \sum_{(1|j) \otimes (1|i) = 1} a^1_j \\ & = \sum_{j \otimes i = 0} a^0_j + \sum_{j \otimes i = 1} a^1_j - \sum_{j \otimes i = 1} a^0_j - \sum_{j \otimes i = 0} a^1_j \\ & = A^0_i - A^1_i \end{split}$$

따라서 아래와 같다.

$$ A = (A^0 + A^1, A^0 - A^1)$$

$A=(A', A'')$일 때 역연산은 아래와 같다.

$$\begin{split} A^0 & = \frac{A' + A''}{2} \\ A^1 & = \frac{A' - A''}{2} \end{split}$$

이를 구현하면 아래와 같다.

void fwht_xor(vector<ll> &a, bool inv) {
    ll n = a.size();
    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;
            }
        }
    }
}