FWT (Fast Walsh-Hadamard Transform) 快速沃爾什變換

要求 \(C(x)=\sum_{i\otimes j=x}A(i)B(j)\)，\(\otimes\) 是二進位制下的一種變換，如 \(\cap,\cup,\oplus\) 或者一些自定義運算。

類似 \(FFT\)，我們對 \(A(x),B(x)\) 做線性變換變為 \(A'(x),B'(x)\)，然後對應值相乘，即 \(C'(x)=A'(x)B'(x)\)，對 \(C'(x)\) 做逆變換後得到 \(C(x)\) 。

設在變換中，\(F'(x)=\sum_{i=0}^n g(x,i)F(i)\)：

\[\begin{align} C'(x)&=A'(x)B'(x) \\ \sum_{i=0}^n g(x,i)C(i)&=\sum_{j=0}^n g(x,j)A(j)\sum_{k=0}^ng(x,k)B(k) \\ \sum_{i=0}^ng(x,i)\sum_{j\otimes k=i}A(j)B(k)&=\sum_{j=0}^n\sum_{k=0}^n g(x,j)g(x,k)A(j)B(k) \\ \sum_{j=0}^n\sum_{k=0}^n A(j)B(k)g(x,j\otimes k)&=\sum_{j=0}^n\sum_{k=0}^n A(j)B(k)g(x,j)g(x,k)\\ g(x,i\otimes j)&=g(x,i)g(x,j) \end{align} \]

第 4 步到第 5 步可以互推，取特定的 \(A(j)=B(k)=1\)

，剩下為 \(0\) 即證。

於是若 \(g(x,i)\) 滿足:

\[g(x,i\otimes j)=g(x,i)g(x,j) \]

就可以作為 \(FWT\) 變換的一種方案，當然，由於我們在求值，\(g(x,i)\) 也需要是可以逆變換的。

這裡主要考慮按位的二進位制變換（詭異如 \(i\otimes j=\text{count}(i)\cdot \text{count}(j)\) 這種通常不太可做），那麼其實只要使得 \(0\) 和 \(1\) 處的關係滿足上式，然後令最終的對應 \(g(x,i)\) 是每一位的 \(g'(x,i)\) 之乘積即可。同時由於要可逆，變換系數應是非平凡的，於是乎對於常見的 \(\otimes\)

運算就很好構造了：

• \(\cup\)： 即 \(g(x,i)=[i\sub x]\)

• \(\cap\)： 即 \(g(x,i)=[i\supset x]\)

• \(\oplus\)：稍微麻煩一些，不過注意到

  \[\begin{cases} g(x,0)=g(x,0)^2=g(x,1)^2\\ g(x,1)=g(x,0)g(x,1) \end{cases} \]

  而且 \(x\) 分別為 \(0,1\) 時兩者的係數還要不一樣，

  於是：\(g(x,i)=(-1)^{|x\cap i|}\)

  其它的一些對應方案不一定可以，如 \(((1,-1),(-1,1)),\;g(x,i)=(-1)^{|x\oplus i|}\)

  ，這就沒法求逆。

 

現在我們已經會構造合適的線性變換 \(g(x,i)\) 了，考慮對 \(F(x)\) 求出 \(F'(x)\) 。這裡可以借鑑 \(FFT\) 的思路，從低位到高位逐步變換，可以發現，實際上就是 \(g(x,i)\) 每次多乘上當前位的 \(g'(x,i)\) 的係數。

設 \(a_0,a_1\) 為當前位是 \(0\) 和 \(1\) 的部分對應的兩個序列，那麼有：

\[F'(a)=\text{merge}(g'(0,0)F'(a_0)+g'(0,1)F'(a_1),\;g'(1,0)F'(a_0)+g'(1,1)F'(a_1)) \]

而逆變換則類似矩陣求逆：

\[a=\text{merge}\left(\frac{g'(1,1)a_0-g'(0,1)a_1}{g'(0,0)g'(1,1)-g'(0,1)g'(1,0)},\frac{g'(1,0)a_0-g'(0,0)a_1}{g'(1,0)g'(0,1)-g'(1,1)g'(0,0)}\right) \]

每一位的貢獻是獨立的，因此逆變換不用倒過來從高位到低位還原，和正變換從低到高也是可以的。

於是每次 \(O(n\log n)\) 求解即可。

 

三者的 code：

void OR(ll *a, int op){
    for(int i = 1; i < 1<<n; i <<= 1)
        for(int j = 0; j < 1<<n; j += i<<1)
            rep(k,0,i-1) (a[i+j+k] += a[j+k] * op) %= mod;
}
void AND(ll *a, int op){
    for(int i = 1; i < 1<<n; i <<= 1)
        for(int j = 0; j < 1<<n; j += i<<1)
            rep(k,0,i-1) (a[j+k] += a[i+j+k] * op) %= mod;
}
void XOR(ll *a, int op){
    for(int i = 1; i < 1<<n; i <<= 1)
        for(int j = 0; j < 1<<n; j += i<<1) rep(k,0,i-1){
            ll v1 = a[j+k], v2 = a[i+j+k];
            a[j+k] = (v1 + v2) % mod, a[i+j+k] = (v1 + mod - v2) % mod;
            if(op == -1) (a[j+k] *= inv2) %= mod, (a[i+j+k] *= inv2) %= mod;
        }
}

 

一個應用：

HDOJ 6966 I love sequences

給定序列 \(\{a_n\},\{b_n\}\)，求 \(\displaystyle \sum_{p=1}^n\sum_{k=1}^{+\infty} d_{p,k}\cdot c_{p}^k\)，其中 \(\displaystyle d_{p,k}=\sum_{k=i\oplus j}a_i\cdot b_j\;(1\leq i,j\leq \frac{n}{p})\)，\(\oplus\) 為三進位制下按位取 \(\gcd\) 。\(n\leq 2\times 10^5\)

還是按位考慮，\(i,j,\gcd(i,j)\) 是這麼個關係：

對於 \(i=0,1,2\) 都滿足 \(g(x,i)^2=g(x,i)\)，那麼取值只能是 \(0\) 或 \(1\) 。可以發現滿足條件的 \(g(x,*)\) 有三組：\((1,0,1),(1,1,1),(1,0,0)\)，於是將這三組任意分配給三行即可：

\[g=\begin{pmatrix} 1&0&1\\ 1&1&1\\ 1&0&0 \end{pmatrix} \]

逆變換就是對矩陣 \(g\) 求逆：

\[g^{-1}=\begin{pmatrix} 0&0&1\\ -1&1&0\\ 1&0&-1 \end{pmatrix} \]

求一個 \(d_p\) 是 \(O(len\log len)\) 的，\(\sum len = O(n\ln n)\)，所以時間複雜度是 \(O(n\log^2 n)\) 。

HDOJ 關掉前寫的程式碼，\(g\) 的順序稍有不同。

#include<iostream>
#include<cstdio>
#include<cstring>
#define rep(i,a,b) for(int i = (a); i <= (b); i++)
#define per(i,b,a) for(int i = (b); i >= (a); i--)
#define N 660000
#define mod 1000000007
#define ll long long
using namespace std;

inline int read(){
    int s = 0, w = 1;
    char ch = getchar();
    while(ch < '0' || ch > '9'){ if(ch == '-') w = -1; ch = getchar(); }
    while(ch >= '0' && ch <= '9') s = (s<<3)+(s<<1)+(ch^48), ch = getchar();
    return s*w;
}

ll A[N], B[N], C[N], D[N], tA[N], tB[N];
int n, m;

void Trans(ll &x, ll &y, ll &z){
    ll a = x, b = y, c = z;
    x = a, y = (a+b+c)%mod, z = (a+c)%mod;
}
void iTrans(ll &x, ll &y, ll &z){
    ll a = x, b = y, c = z;
    x = a, y = (b-c+mod)%mod, z = (c-a+mod)%mod;
}

void FWT(ll *a, int id){
    for(int i = 1; i < m; i *= 3)
        for(int j = 0; j < m; j += (i*3)) rep(k,0,i-1){
            if(id == 0) Trans(a[j+k], a[j+k+i], a[j+k+2*i]);
            else iTrans(a[j+k], a[j+k+i], a[j+k+2*i]);
        }
}

int main(){
    n = read();
    rep(i,1,n) tA[i] = A[i] = read();
    rep(i,1,n) tB[i] = B[i] = read();
    rep(i,1,n) C[i] = read();

    ll ans = 0;
    rep(p,1,n){
        int lim = n/p;
        rep(i,1,lim) A[i] = tA[i], B[i] = tB[i];
        m = 1;
        while(m <= lim) m *= 3;
        
        FWT(A, 0), FWT(B, 0);
        rep(i,0,m-1) D[i] = (A[i]*B[i])%mod;
        FWT(D, 1);

        ll pow = 1;
        rep(i,1,m) (pow *= C[p]) %= mod, (ans += pow*D[i]) %= mod;
        rep(i,0,m) A[i] = B[i] = D[i] = 0;
    }

    cout<<ans<<endl;
}

這是官方題解：

每個數列採用的 \(FWT\) 變換矩陣甚至可以不同，算是從一個更高的角度看待快速沃爾什變換了。