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 ceil}$

初始情况下有 $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(aoplus 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( egin{array}{ccc} H{n-1} & H{n-1} \ H{n-1} & -H{n-1} end{array} ight)$$

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

#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;
}