上次的博客有点模糊的说...我把思路和算法实现说一说吧...
思路
关于快速沃尔什变换,为了方便起见,我们采用线性变换(非线性变换不会搞)。
那么,就会有一个变化前各数值在变换后各处的系数,即前一篇博文中的$f(i,j)$,表示线性变换中第$i$项到第$j$项的系数。
即
$$DWT(A)_i = sum_{j=0}^{n-1} A_j * f(i,j)$$
那么,我们既然要求$oplus$卷积在变换后等价于乘积,就有
$$DWT(A)_i * DWT(B)_i = DWT(C)_i$$
其中$C$是$A,B$的$oplus$卷积。
那么,将线性变换式及卷积定义代入,并令对应项系数相等,就可得到结论:
$$f(i, j) * f(i, k) = f(i, j oplus k)$$
那么,通过构造合适的$f(i,j)$,就可得出变换。
$f$函数的构造
我们以异或为例(因为这时最麻烦的)。
既然是位运算,那么我们就从位运算入手。
考虑位运算,我们只需要考虑一位的情况,对于多位运算需要将每位的结果相乘(注意是相乘,这很重要!)。
通过定义,我们得到
$$f(i, 0) * f(i, j) = f(i, 0 oplus j) = f(i, j)$$
$$f(i, 1) * f(i, 1) = f(i, 0)$$
又由于我们不能使所有系数都相同,那么只能取
$$egin{cases}
f(i, 0) &= 1\
f(0, j) &= 1\
f(1, 1) &= -1
end{cases}$$
那么此时
$$DWT(a_0, a_1) = (a_0 + a_1, a_0 - a_1)$$
可以验证其满足规则。
算法实现:
要实现$FWT$,我们仿照快速傅里叶变换,首先将$A$分成两半$A_0, A_1$(按二进制最高位,即直接取前一半和后一半)递归调用$FWT$,然后,
合并时,根据下面的式子($i_0,i_1$分别表示$i$的最高位和剩余位):
$$egin{aligned}
DWT(A)_i &= sum_{j=0}^{n-1} f(i, j)A_j\
&= sum_{j=0}^{n/2-1}f(i,j)A_j + sum_{j=n/2}^{n-1}f(i,j)A_j\
&= sum_{j=0}^{n/2-1}f(i_0,j_0)f(i_1,j_1)A_j + sum_{j=n/2}^{n-1}f(i_0,j_0)f(i_1,j_1)A_j\
&= sum_{j=0}^{n/2-1}f(i_0,0)f(i_1,j_1)A_j + sum_{j=n/2}^{n-1}f(i_0,1)f(i_1,j_1)A_j\
&= f(i_0, 0)sum_{j=0}^{n/2-1}f(i_1,j_1)A_j + f(i_1)sum_{j=n/2}^{n-1}f(i_0,1)f(i_1,j_1)A_j\
&= f(i_0, 0)*DWT(A_0)_{i_1} + f(i_0, 1)*DWT(A_1)_{i_1}end{aligned}$$
其中由第2行到第3行是因为我们对多位的$f$直接定义为各位的$f$相乘,即
$$f(i,j)=f(i_0,j_0)*f(i_1,j_1)$$
那么,算法实现就简单多了。
FWT(A, len)
1 divide A into (A0, A1)
2 FWT(A0, len/2)
3 FWT(A1, len/2)
4 for i=0 to len/2
5 A[i]=f00*A0[i]+f01*A1[i]
6 A[i+len/2]=f10*A0[i]+f11*A1[i]
关于逆变换嘛。。。你只需要告诉自己:我只是要处理刚被$FWT$变换过的数组,我只需要把它倒过来。
那么
IFWT(A, len)
1 divide A into (A0, A1)
2 for i=0 to len/2
3 solve the equation set, where the unknown numbers are A0[i] and A1[i]:
4 f00*A0[i]+f01*A1[i]=A[i]
5 f10*A0[i]+f11*A1[i]=A[i+len/2]
6 IFWT(A0, len/2)
7 IFWT(A1, len/2)
(当然解二元一次方程组只需要手算出来就好啦)
把三种运算(即与,或,异或)的f给出来:
$$egin{array}{c|cccc}
Ops & f_{00} & f_{01} & f_{10} & f_{11} \
hline
And & 1 & 1 & 0 & 1 \
Or & 1 & 0 & 1 & 1 \
Xor & 1 & 1 & 1 & -1
end{array}
$$
附三种运算变换代码:
1 void FWT_And(int *A, int len) { 2 if (len == 1) return; 3 int len2 = len >> 1; 4 FWT_And(A, len2); 5 FWT_And(A + len2, len2); 6 for (int i = 0; i < len2; ++i) 7 A[i] += A[i + len2]; 8 } 9 void IFWT_And(int *A, int len) { 10 if (len == 1) return; 11 int len2 = len >> 1; 12 for (int i = 0; i < len2; ++i) 13 A[i] -= A[i + len2]; 14 IFWT_And(A, len2); 15 IFWT_And(A + len2, len2); 16 } 17 18 void FWT_Or(int *A, int len) { 19 if (len == 1) return; 20 int len2 = len >> 1; 21 FWT_Or(A, len2); 22 FWT_Or(A + len2, len2); 23 for (int i = 0; i < len2; ++i) 24 A[i + len2] += A[i]; 25 } 26 void IFWT_Or(int *A, int len) { 27 if (len == 1) return; 28 int len2 = len >> 1; 29 for (int i = 0; i < len2; ++i) 30 A[i + len2] -= A[i]; 31 IFWT_Or(A, len2); 32 IFWT_Or(A + len2, len2); 33 } 34 35 void FWT_Xor(int *A, int len) { 36 if (len == 1) return; 37 int len2 = len >> 1; 38 FWT_Xor(A, len2); 39 FWT_Xor(A + len2, len2); 40 for (int i = 0; i < len2; ++i) { 41 int x = A[i], y = A[i + len2]; 42 A[i] = x + y; 43 A[i + len2] = x - y; 44 } 45 } 46 void IFWT_Xor(int *A, int len) { 47 if (len == 1) return; 48 int len2 = len >> 1; 49 for (int i = 0; i < len2; ++i) { 50 int x = A[i], y = A[i + len2]; 51 A[i] = (x + y) >> 1; 52 A[i + len2] = (x - y) >> 1; 53 } 54 IFWT_Xor(A, len2); 55 IFWT_Xor(A + len2, len2); 56 }