Fast Walsh–Hadamard transform


考虑变换

$$\hat{A_x} = \sum_{i\ or\ x = x}{ A_i }$$

记 $S_{t}(A,x) = \sum_{c(i,t)\ or\ c(x,t)=c(x,t),\ i \le |A|}{A_i}$
则 $\hat{A} = S_{\lceil log_2n \rceil}$

初始情况下有 $S_0$ 被拆分为 $n$ 段, $S_0([A_i],i) = A_i$
考虑每次将相邻两段合并。

记 $B0 = S_t([A0,A1],x)$ 的前一半,记 $B1$ 为后一半。
则有
$$B0 = A0 \\ B1 = A0 + A1$$
从而考虑将序列长度补为 $2^t$,不停合并序列即可实现时间复杂度 $O(nt)$ 的变换方法,逆变换只需要按上述变换解方程即可。

考虑应用这个变换
试求$$C_x = \sum_{a\ or\ b = x} A_a B_b$$
考虑应用变换的性质
$$\hat{C_c} \\ = \sum{[x\ or\ c=c]C_x} \\ = \sum{[(a\ or\ b)\ or\ c = c]A_a B_b} \\= \sum{[(a\ or\ c = c)]A_a [(b\ or\ c = c)B_b ]} \\=\hat{A_c} \hat{B_c}$$

对于
$$C_x = \sum_{a\ and\ b = x}{A_a B_b}$$
我们考虑构造另外一种变换方法,$$[(a\ and\ b)\ and\ c = c] = [(a\ and\ c = c)][(b\ and\ c = c)]$$
新的变换方法同上,为$$B0 = A0 + A1 \\ B1 = A1$$
上述变换很容易推广到高维的情况,略。

那么考虑另外一种逻辑规则$$C_x = \sum_{a \oplus b = x} A_a B_b$$
类比上述两式我们应用本身的运算 $xor$,再考虑将其化为只含有0,1的逻辑函数。
尝试定义 $f(a,b) = a \oplus b$ 中一的个数取余2,从而有
$$([f(a,c)=0]-[f(a,c)=1])([f(b,c)=0] - [f(b,c)=1]) \\= [f(a \oplus b,c) = 0] - [f(a \oplus b,c) = 1]$$
我们将变换称为 $w$ 变换
从而 $w(C) = w(A)\cdot w(B)$
考虑如何求解 $w(A)$
假定,记$$B0 = B00 - B10 \\ B1 = B01 - B11$$
从而
$$B00 = A00 + A11, B01 = A10 + A01 \\ B10 = A10 + A01, B11 = A00 + A11 \\ B0 = A00+A11-A10-A01 \\ B1 = A10+A01-A00-A11$$
从而$$B0 = A0 - A1, B1 = A1 - A0$$
正确但是无法通过反解的方式进行逆变换。
可以发现关键问题在于我们要设法使得$B0, B1$的生成式不能线性相关,从而不可以使得$$f(a,b) = g(a \oplus b)$$
因为前者会导致两个递归式线性相关。
而使得$$f(a,b) = g(a\ or\ b) / g(a\ and\ b)$$
则可以得到两个并不线性相关的式子。
考虑修改 $f(a,b)$ 定义,经尝试得到$f(a,b) = a\ and\ b$ 中 1 的个数满足$$f(a\oplus b,c) = f(a,c) \oplus f(b,c)$$契合上述式子
这样我们会得到变换$$B0 = A0+A1 \\ B1 = A0-A1$$
从而逆变换为$$A0 = \frac{B0+B1}{2} \\ A1 = \frac{B0-B1}{2}$$
考虑和$fft$一样,统一正反变换,得到
$$B0 = \frac{A0+A1}{\sqrt 2} \\ B1 = \frac{A0-A1}{\sqrt 2}$$
反向变换相同。
我们考虑将原序列中各个元素对于变换后序列每一项的贡献,用矩阵写出来得到$Hadamard$ 矩阵。
$$H_n = \frac{1}{\sqrt 2}\left( \begin{array}{ccc} H{n-1} & H{n-1} \\ H{n-1} & -H{n-1} \end{array} \right)$$

后两个变换的简单测试程序:

#include <bits/stdc++.h>
using namespace std;
void W(int a[],int l,int r)
{
    if(l==r) return;
    int m = (l+r)>>1;
    W(a,l,m);
    W(a,m+1,r);
    for(int i=l;i<=m;i++)
    {
        a[i]-=a[i-l+1+m];
        a[i-l+1+m] = -a[i];
    }
}
void FWT(int a[],int l,int r)
{
    if(l==r) return;
    int m = (l+r)>>1;
    FWT(a,l,m);
    FWT(a,m+1,r);
    for(int i=l;i<=m;i++)
    {
        int x = a[i], y = a[i-l+1+m];
        a[i] = x+y;
        a[i-l+1+m] = x-y;
    }
}
void rFWT(int a[],int l,int r)
{
    if(l==r) return;
    int m = (l+r)>>1;
    for(int i=l;i<=m;i++)
    {
        int x = a[i], y = a[i-l+1+m];
        a[i] = (x+y)/2;
        a[i-l+1+m] = (x-y)/2;
    }
    rFWT(a,l,m);
    rFWT(a,m+1,r);
}
int sol(int now)
{
    int ans = 1;
    for(;now;now>>=1) if(now&1) ans=-ans;
    return ans; 
}
#define N 100010
int a[N],s[N];
int main()
{
    int t;
    cin >> t;
    for(int i=0;i<(1<<t);i++) scanf("%d",&a[i]);
    for(int i=0;i<(1<<t);i++)
    {
        s[i] = 0;
        for(int j=0;j<(1<<t);j++) s[i] += sol(i^j)*a[j];
    }
    W(a,0,(1<<t)-1);
    for(int i=0;i<(1<<t);i++) cout << s[i] - a[i] << endl;
    return 0;
}
View Code

 

posted @ 2018-07-10 10:21  lawyer'  阅读(207)  评论(0编辑  收藏  举报