1

FWT/快速沃尔什变换 入门指南 - _maze

 1 year ago
source link: https://www.cnblogs.com/closureshop/p/how-to-use-FWT.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

来学点好玩的。


我们也许学过,FFT 可以解决一类卷积:

Ci=∑k+j=iAiBj

现在我们稍微变一下式子:

Ci=∑i=k&jAkBj
Ci=∑i=k∣jAkBj
Ci=∑i=k⊕jAkBj
上面那个圆圆的东西是异或

FWT 或

也就是这个式子:

Ci=∑i=k∣jAkBj

FWT 的定义

模仿 FFT,对于 A 我们想要在可接受时间内得到一个 FWT(A),使得

FWT(C)i=FWT(A)i×FWT(B)i

i 代表第 i 个位置的数。

这样我们就可以 O(n) 暴力乘起来了。

因此我们构造 FWT(A)i=∑j|i=iAj。

现在我们来证明这个构造的正确性。

FWT(A)i×FWT(B)i
=∑j|i=iAj×∑k|i=iBk
=∑j|i=i∑k|i=iAjBk

因为可以从 j|i=i 与 k|i=i 中得出 (j|k)|i=i,所以我们还可以消去另一个式子

=∑(j|k)|i=iAjBk

根据 FWT 的定义,这个式子就是 FWT(C)i。证毕。

另一种理解

涉及到了或运算,我们不妨把数全都变成二进制。于是我们想到一种经典转换:将二进制中的 0 与 1 转化为集合中一个数选还是不选。那么或操作代表什么呢?

bingo!两个集合的并集!

于是我们可以把上面的式子改写一个形式:

Ci=∑i=k∣jAkBj
Ci=∑i=k⋃jAkBj

注意看,这时 i,j,k 都是集合,只不过我们将其用二进制表示.

这样我们可以改写 FWT 的定义。

重新来一遍,FWT(A)i=∑k⊆iAk。

因此我们为 FWT 找到了新的定义,他代表着集合 i 的子集之和。我们有了更自然的推导:

∑j⊆iAj×∑k⊆iBk
=∑j,k⊆iAj×Bk
=∑x⊆i∑j⋃k=xAj×Bk
=∑x⊆iCx
=FWT(c)x

换句话说,我们对子集做了个前缀和操作(发现了吗?子集的前缀和进行或卷积与普通前缀和进行加法卷积具有相似性),并用前缀和相乘代替了原来的 O(n2) 相乘。

我们把原序列 A 按下标最高位是 0 还是 1 分成两部分 A0 与 A1 分治求解。显然,前半部分(最高位为 0 的部分)就是 FWT(A0),所以我们考虑后半部分的答案。

后半部分最高位为 1,因此此时“子集”这一概念不仅包含分治处理的他子集,还包括把最大值变为 0 后的,序列 A0 中同一位置的子集。要将 A0 中的同一位置加到当前答案上。

写成数学形式就是:

FWT(A)=merge(FWT(A0),FWT(A0)+FWT(A1))

上面的 merge 代表拼接,就是字面意思。

于是我们就能写出分治递归代码了!但为了常数着想,我们试着把递归这一步骤去掉。

去掉的部分并不难写,我们按照层数从小到大递归,不难发现第 i 层(从 0 开始编号,最底层为 0)就是枚举第 i 位是 0 还是 1,并且乱填其他数进行转移。

dp1xde4a.png

代码也简单:

void OR(mint *f){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)// k 为当前的层,o 仅用于穷举左边进行转移
        for(int i = 0; i < n; i += o)// 穷举左边
            for(int j = 0; j < k; j++){ // 穷举右边
                f[i + j + k] = f[i + j + k] + f[i + j];
            }
}

如何转回来

再看一眼转移的式子

FWT(A)=merge(FWT(A0),FWT(A0)+FWT(A1))

思考只有两个数的情况。此时 1 位置是不会变的,2 位置加上了 1 位置的贡献,要减去。

我们发现更大的情况也是一样的,只要依次把前面的贡献减去就好。

void IOR(mint *f){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)// k 为当前的层,o 仅用于穷举左边进行转移
        for(int i = 0; i < n; i += o)// 穷举左边
            for(int j = 0; j < k; j++){ // 穷举右边
                f[i + j + k] = f[i + j + k] - f[i + j];
            }
}

这两份代码显然是可以合并的。因此我们得到了 FWT 或 的全过程。

void OR(mint *f, int type){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j++){
                f[i + j + k] = f[i + j + k] + (f[i + j] * mint(type));
            }
}

FWT 与

和 或 差不多,只是要从 1 转移到 0。

可以发现,实际上我们用子集后缀和优化了运算。

void AND(mint *f, int type){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j++)
                f[i + j] = f[i + j] + (f[i + j + k] * mint(type));
}

FWT 异或

Ci=∑i=k⊕jAkBj

很遗憾,我并没有发现这个东西的集合意义,如果有大佬知道可以告诉我。。。

思考 FWT 的作用,我们想要把 AkBj 变成 AiBi 的形式,以此来简化运算。

我们考虑这样 n 个 n 维向量 b,b(i) 只有下标 i 处是 1,其他位置都是 0。

现在我们把 FWT 后的 A,B 看作系数,此时显然 A1b(1),A2b(2),...,Anb(n)=A1,A2,A3,A4...,An

显然,异或卷积对于乘法有分配律。

设异或卷积为 ∗,则

(∑i=1nAib(i))∗(∑i=1nBjb(j))
=(∑i=1n∑j=1nAiBj(b(i)∗b(j))

发现后面的东西可以简单表示,即 bi∗bi=bi,bi∗bj=0(i≠j)。

那么整个式子就是我们寻找的形式:

∑i=1nAiBi

而我们要做的事情无非是求出 FWT 之前的 bi。

太长不看版:异或 FWT 与原序列线性相关

既然这样,我们设 FWT(A)x=∑i=1ng(x,i)Ai。

那么因为 FWT(C)x=FWT(A)x×FWT(b)x

所以 ∑k=1ng(x,k)Ck=∑i=1ng(x,i)Ai×∑j=1ng(x,j)Bj。

整理一下可以得出:

∑k=1ng(x,k)Ck=∑i=1n∑j=1ng(x,i)g(x,j)×AiBj

把 Ck 用 A,B 表示可得:

∑k=1ng(x,k)∑k=i⊕jAiBj=∑i=1n∑j=1ng(x,i)g(x,j)×AiBj

更改求和顺序,我们枚举 i,j 可得:

∑i=1n∑j=1ng(x,i⊕j)AiBj=∑i=1n∑j=1ng(x,i)g(x,j)×AiBj

于是我们发现了 g 的关系:

g(x,i⊕j)=g(x,i)g(x,j)

现在问题来了,与 i,j 相关的什么东西,使异或之后的值等于原来两值的乘积?

于是我们可以想到有人托梦给我奇偶性。

具体的,我们发现异或前后 1 的个数奇偶性不变。原因如下:

按每一位依次考虑。如果第 i 位异或后为 1,那么原来必定有且仅有一个 1。个数不变

如果为 0,要么是两个 0,此时 1 的个数不变,要么是两个 1,此时 1 的个数减 2,奇偶性仍不变。

所以我们定义 g(x,i)=(−1)|i⋂x|。那么上式就等价于:

(−1)|(i⊕j)⋂x|=(−1)|i⋂x|(−1)|j⋂x|

根据上面的推论,左右两边奇偶性不变,与 后无非是减去两个相同的数,奇偶性还是不变。

于是我们得出 FWT 的转移式:

FWT(A)x=∑i=1n(−1)|i⋂x|Ai

考虑模仿前两个 FWT 的形式,讨论最高位 i 为 0 和为 1 两种情况。

原来最高位为 0,FWT 后的前 2i−1 个数最高位还是 0。由于 1&0=0,所以后 2i−1 个数的贡献为正。前半部分答案为 FWT(A0)+FWT(A1)。

FWT 后的后 2i−1 个数最高位变成了 1,此时 A0 的贡献还是正(因为 1&0=0)。但是此时后半部分加了 1,于是贡献要取反。后半部分答案为 FWT(A0)−FWT(A1)。

所以我们得出:

FWT(A)=merge(FWT(A0)+FWT(A1),FWT(A0)−FWT(A1))

当然,参照 或 FWT,我们可以写出不依赖递归的程序:

void XOR(mint *f){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)//具体意义参考 FWT 或
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j ++){
                mint x = f[i + j], y = f[i + j + k];
                f[i + j] = x + y;
                f[i + j + k] = x - y;
            }
}

实际上就是把贡献减去

IFWT(A)=merge(IFWT(A0)+IFWT(A1)2,IFWT(A0)−IFWT(A1))2

显然这两个东西是可以合并的。于是我们可以得出模板的完整代码:

#include<bits/stdc++.h>
using namespace std;

#define forp(i, a, b) for(int i = (a);i <= (b);i ++)
#define forc(i, a, b) for(int i = (a);i >= (b);i --)

const int maxn = 6e5 + 5;
const int mod = 998244353;

int read(){
    int u;cin >> u;return u;
}

class mint{
    private : int v;
    public:
        mint(){}
        int operator()(void)const{
            return v;
        }
        mint (const int &u){ 
            v = u % mod; 
        }
        mint operator+(const mint &a) const{ 
            int x = a.v + v;
            if(x >= mod) return mint(x - mod);
            if(x < 0) return mint(x + mod);
            return x;
        }
        mint operator-(const mint& a)const{
			return v < a.v ? v - a.v + mod : v - a.v;
		}
        mint operator*(const mint &a) const{
            return mint((1ll * a.v * v) % mod);
        }
};

mint qpow(mint u, int v){
    mint ans = mint(1);
    while(v){
        if(v & 1) ans = ans * u;
        u = u * u;
        v >>= 1;
    }
    return ans;
}
mint inv2 = qpow(2, mod - 2);

int n;
mint A[maxn], B[maxn], C[maxn];
mint g[maxn];

void OR(mint *f, int type){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j++){
                f[i + j + k] = f[i + j + k] + (f[i + j] * mint(type));
            }
}

void AND(mint *f, int type){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j++)
                f[i + j] = f[i + j] + (f[i + j + k] * mint(type));
}

void XOR(mint *f, int type){
    for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
        for(int i = 0; i < n; i += o)
            for(int j = 0; j < k; j ++){
                mint x = f[i + j], y = f[i + j + k];
                f[i + j] = x + y;
                f[i + j + k] = x - y;
                if(type == -1){
                    f[i + j] = f[i + j] * inv2;
                    f[i + j + k] = f[i + j + k] * inv2;
                }
            }
}

signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);

    n = (1 << read());
    forp(i, 0, n - 1) A[i] = mint(read());
    forp(i, 0, n - 1) B[i] = mint(read());

    OR(A, 1);OR(B, 1);
    forp(i, 0, n - 1) C[i] = (A[i] * B[i]);
    OR(C, -1);
    forp(i, 0, n - 1) cout << C[i]() << ' ';
    cout << endl;
    OR(A, -1);OR(B, -1);

    AND(A, 1);AND(B, 1);
    forp(i, 0, n - 1) C[i] = (A[i] * B[i]);
    AND(C, -1);
    forp(i, 0, n - 1) cout << C[i]() << ' ';
    cout << endl;
    AND(A, -1);AND(B, -1);

    XOR(A, 1);XOR(B, 1);
    forp(i, 0, n - 1) C[i] = (A[i] * B[i]);
    XOR(C, -1);
    forp(i, 0, n - 1) cout << C[i]() << ' ';
    cout << endl;
    XOR(A, -1);XOR(B, -1);
    return 0;
}

大概就是这样。。。应用也许会另外开坑吧。


感谢 xht 的博客,从这篇博客里我学到了 FWT 的基础知识。

感谢同校大佬 yllcm 为本人解释符号与定义。

感谢万能的U群群友 Untitled_unrevised 解释 FWT 的目的。

感谢万能的U群群友 rqy 学姐与 樱初音斗橡皮 解释为什么异或 FWT 是线性变换顺便发现了我在线性代数方面的巨大缺口


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK