天天看点

FWT

另一种卷积

多项式乘法的卷积是这样的:

Ci=∑j+k=iAjBk

然而还有种神奇的卷积( ⊕ 是位运算):

Ci=∑j⊕k=iAjBk

FWT可以快速求这种卷积。

快速沃尔什变换

因为FFT用了系数 → 点值,用点值快速求卷积,再点值 → 系数的方法,所以我们也想让沃尔什变换具有这样的特点,于是开始构造。

离散沃尔什变换

定义 DWT(A) 是向量 A 经过离散沃尔什变换得到的向量,令变换系数为 f(i,j) ,那么:

DWT(A)i=∑j=0n−1Ajf(i,j)DWT(A)iDWT(B)i=DWT(C)i

这样没法求 f(i,j) 啊,我们来推一下:

DWT(A)iDWT(B)i=DWT(C)i∑j=0n−1Ajf(i,j)∑k=0n−1Bkf(i,k)=∑t=0n−1Ctf(i,t)∑j=0n−1∑k=0n−1AjBkf(i,j)f(i,k)=∑t=0n−1∑j⊕k=tAjBkf(i,t)=∑t=0n−1∑j⊕k=tAjBkf(i,j⊕k)∑j=0n−1∑k=0n−1AjBkf(i,j)f(i,k)=∑j=0n−1∑k=0n−1AjBkf(i,j⊕k)

也就是说 f(i,j)f(i,k)=f(i,j⊕k) 。

然后就可以对不同的位运算求出不同的变换系数,直接上结论:

and:f(i,j)=[i and j=i]or:f(i,j)=[i and j=j]xor:f(i,j)=(−1)count(i and j)

其中 count(i) 表示 i 在二进制下 1 的个数。这些怎么得到的?翰爷表示:

并不知道是怎么来的,但是容易验证是对的。

这里就说明一下比较神奇的 xor ( and 和 or 很容易验证):

f(i,j)f(j,k)=(−1)count(i and j)+count(i and k)f(i,j xor k)=(−1)count(i and (j xor k))

count(i and j)+count(i and k) 求的是 i 和 j 共有的 1 的数量加上 i 和 k 共有的 1 的数量,我们再考虑下面那个式子表示什么。因为异或相同为 0 ,不同为 1 ,那么只有当 j 和 k 某位上不同时才有贡献,和上面的式子相比,消去了 j 和 k 某位上均为 1 的贡献,即 −2k ,不改变奇偶性,所以是相等的。

f(i,j) 有一个很重要性质:可以二进制拆分。用 ik 表示 i 二进制下的第 k+1 高位,则:

f(i,j)=f(i0,j0)f(i1,j1)f(i2,j2)⋯f(ilen−1,jlen−1)

接下来开始考虑离散沃尔什变换( n 补到 2k ):

DWT(A)i=∑j=0n−1Ajf(i,j)=∑j=0n2−1Ajf(i,j)+∑j=n2n−1Ajf(i,j)=∑j=0n2−1Ajf(i0,j0)f(i[1,len−1],j[1,len−1])+∑j=n2n−1Ajf(i0,j0)f(i[1,len−1],j[1,len−1])=f(i0,0)∑j=0n2−1Ajf(i[1,len−1],j[1,len−1])+f(i0,1)∑j=n2n−1Ajf(i[1,len−1],j[1,len−1])

喜闻乐见的规模减半了,于是( i∈[0,n2) ):

DWT(A)i=f(0,0)DWT(AL)i+f(0,1)DWT(AR)iDWT(A)i+n2=f(1,0)DWT(AL)i+f(1,1)DWT(AR)i

和FFT很像啊,只不过FFT是奇偶拆半,而FWT是左右拆半。正因为FWT是左右拆半,直接迭代就可以了,不需要和FFT一样二进制翻转。

逆沃尔什变换

根据离散沃尔什变换得到的递归式,我们发现只需要解一下二元一次方程就可以了。

模板

inline void FWT(int *a,int n,int f){
    for (int k=;k<n;k<<=)
        for (int i=;i<n;i+=(k<<))
            for (int j=;j<k;j++)
                if (f==){
                    int x=a[i+j],y=a[i+j+k];
                    //and:a[i+j]+=a[i+j+k];
                    //or :a[i+j+k]+=a[i+j];
                    //xor:a[i+j]=x+y;a[i+j+k]=x-y;
                } else{
                    int x=a[i+j],y=a[i+j+k];
                    //and:a[i+j]-=a[i+j+k];
                    //or :a[i+j+k]-=a[i+j];
                    //xor:a[i+j]=(x+y)/2;a[i+j+k]=(x-y)/2;
                }
}