FWT

Fast Walsh-Hadamard Transform

.pre

今天本來想看FFT的應用的...翻翻picks的部落格發現了好東西

Fast Walsh-Hadamard Transform

,感覺挺有趣的...就寫了一下...orz

.intro

Fast Walsh-Hadamard Transform(

FWT

)是用來解決一類卷積問題的.((當整篇文章看完後)如果你還是不知道什麼樣的卷積能做,戳這裡,最後一篇vfk大大的論文

集合幂級數與xxxxxx

)

準确地說,我們給出兩個長度都為\(2^m\)的序列\(A_0\dots A_{2^m-1}\),\(B_0\dots B_{2^m-1}\),顯然對于每一個下标我們都可以把它分解成一個唯一的長度為\(m\)的二進制數(開頭補0).我們需要求出一個它們的xor卷積,即一個序列\(C_0..C_{2^m-1}\)滿足:

\[C_k=\sum_{i\otimes j=k}A_iB_j\]

其中\(a\otimes b\)表示\(a\) xor \(b\).

FWT的基本思想就是對每一位分開讨論,恩,因為xor是按位做的,我們顯然可以按位讨論..

然後Picks不知道怎麼做了,恩.

然而窩也不知道怎麼做了,恩.

然後Picks有結論!背結論大法吼:

(等等窩先出個表示法)

由于是一個長度為\(2^m\)的序列,我們可以将此序列按數字大小二分,兩半在原序列裡是連續的且長度相同.顯然此時我們把這個序列按照最高位分開了.我們記這個過程為\(A=(A_0,A_1)\),\(A_0\)就是小的那一半.

我們定義對于一個序列做FWT為\(tf(A)\),然後\(tf(|A|=1,A)=A_0\).\(|A|\)表示\(A\)的長度.

(結論)

這時\(tf(A)=(tf(A_0)+tf(A_1),tf(A_0)-tf(A_1))\).(錯了去找Picks> <)

然後...好像就沒了?

當然FWT隻是對一個序列A進行的操作...其實這相當于一個FFT,要求出卷積還需要再對B最一次tf(b),把結果在同一個位置上的數相乘,然後再做IFFT(口胡,其實是IFWT).

IFWT也很簡單,\(itf(A)=(itf\left(\frac{A_0+A_1}{2}\right),itf\left(\frac{A_0-A_1}{2}\right))\).

(補充一個表示法)

對于兩個等長的序列,我們定義

\[A+B=C \rightarrow C_k=A_k+B_k\]

\[A-B=C \rightarrow C_k=A_k-B_k\]

\[A\times B=C \rightarrow C_k=A_k\times B_k\]

(總結)

那麼我們就會做FWT了,因為

\[A*B=itf(tf(A)\times tf(B))\]

.code

void fwt_xor(ll* x,int r){
    for(int i=2,p=1;i<=r;i<<=1,p<<=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k){
                ll kx=x[k],ky=x[k+p];
                x[k]=kx+ky,x[k+p]=kx-ky;
            }
        }
    }
}
void fwt_xor_inv(ll* x,int r){
    for(int i=r,p=r>>1;p;i>>=1,p>>=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k){
                x[k]=(x[k]+x[k+p])>>1;
                x[k+p]=x[k]-x[k+p];
            }
        }
    }
}                

代碼未經測試,慎用.

.ext

.op.and

\[tf(A)=(tf(A_0)+tf(A_1),tf(A_1))\]

\[itf(A)=(itf(A_0)-itf(A_1),itf(A_1))\]

void fwt_and(ll* x,int r){
    for(int i=2,p=1;i<=r;i<<=1,p<<=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k) x[k]+=x[k+p];
        }
    }
}
void fwt_and_inv(ll* x,int r){
    for(int i=2,p=1;i<=r;i<<=1,p<<=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k) x[k]-=x[k+p];
        }
    }
}                

.op.or

\[tf(A)=(tf(A_0),tf(A_1)+tf(A_0))\]

\[itf(A)=(itf(A_0),itf(A_1)-itf(A_0))\]

void fwt_or(ll* x,int r){
    for(int i=2,p=1;i<=r;i<<=1,p<<=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k) x[k+p]+=x[k];
        }
    }
}
void fwt_or_inv(ll* x,int r){
    for(int i=2,p=1;i<=r;i<<=1,p<<=1){
        int u=r-i;
        for(int j=0,ux=p;j<=u;j+=i,ux+=i){
            for(int k=j;k<ux;++k) x[k+p]-=x[k];
        }
    }
}                

.ps

關于FWT modulo prime

我們可以顯然地處理FWT modulo prime,隻需要取個模就好了(不是廢話麼= =)..

還有...關于如何求\(2^{-1}\equiv x\pmod{p}\)中的\(x\)...直接

x=(p+1)>>2

...顯然滿足性質..

FWT for xor 的模數需要與2互素.

.pps

我可能會出一道裸題..恩就是這個鬼卷積...放到某奇怪的OJ上供測試- -.

轉載于:https://www.cnblogs.com/tmzbot/p/4668158.html