算法学习笔记-FWT

引入

什么是 FWT (快速沃尔什变换) ? 请看这道题: P4717 【模板】 快速莫比乌斯/沃尔什变换 (FMT/FWT)

如果不会 FWT , 最朴素的想法就是两重循环:

for j in range(2 ** n):
    for k in range(2 ** n):
        C[j | k] = (C[j | k] + A[j] * B[k]) % mod

这个做法的时间复杂度是, 无法接受.

有没有什么办法能加速这一过程呢? 这就需要用到 FWT 了.

FWT 核心思想

FWT 的思想和 FFT 非常相似 (没事, 不懂 FFT 、卷积之类的东西也能理解 FWT).

1. 设计一种变换, 将原序列变成一个新序列.

2. 这种变换需要有这样的神奇性质: 我们所要求的(这里用代表位运算卷积) , 它的变换结果, 可以直接通过与逐位相乘 得到. 也就是说,.

3. 这样做是因为, 逐位相乘的复杂度只有, 比卷积快多了.

4. 此外, 我们还需要一个逆变换, 把变成, 毕竟才是我们的最终目标.

整个流程就是:

只要 FWT 和 IFWT 足够快 (剧透一下, 它们能做到), 那么就能大大加速原来的卷积运算.

FWT 公式推导

为了实现 FWT 与 IFWT, 大神们采用了分治思想. (不知道他们是怎么想到这个方法的, 我水平有限, 只会做结论的推导, 无法解释为什么要这样做.)

假设原序列的长度为, 且下标从开始, 我们把对半分成两半:

于是,.

假设我们可以递归地分别对和做 FWT, 那么关键问题就是, 如何从 和合并出?

1. 按 OR 卷积

目标

定义变换

对于 OR 卷积, 我们定义:

何意味

j | i = i 意味着二进制下是的子集, 也就是是所有满足 是的子集的的和.

神秘性质

我们上面说, FWT 需要满足神秘性质. 我们在此处定义的这个变换是否满足该性质?

, 展开后, 每一项是, 其中且.

显然, 如果且, 那么它们的并集也满足.

重写一下求和的下标:

哎! 括号里不就是的定义吗? 于是上式就变成:

因此, 这个定义是合法的.

推导 FWT(A) 的计算方法

回到最开始没有展开说明的关键问题: 如何用和组合出?

分类讨论一下~

1. 对于前半部分(下标从到):

因为, 既然, 那么肯定有.

于是.

2. 对于后半部分(下标从到):

我们令, 其中和上面一样是到.

此时可以来自, 也可以来自, 再分类讨论一下.

• 如果来自, 则. 为了让, 就需要. 这部分的和就是刚刚的

• 如果来自, 我们令, 于是, 所以.

两个部分求和一下,. 注意这里的是序列的逐位相加.

结论:

记,,, 则结论也可以写成:

接下来求逆变换 IFWT 的合并规则, 其实就是已知和, 要求和, 显然,.

结论:

其实, 这里的 FWT 就相当于做 SOS DP.

2. 按 AND 卷积

过程和 OR 卷积类似.

目标

定义变换

对于 AND 卷积, 我们定义:

何意味

j & i = i 意味着二进制下是的子集, 也就是是所有满足 是的子集的的和.

神秘性质

和 OR 卷积一样, 可以验证该定义具有的性质.

推导 FWT(A) 的计算方法

同样是分治.

1. 对于前半部分(下标从到): j 可以来自(记来自的为) 或. 不论是哪种情况, 都可有, 也就是说两部分都有贡献, 所以.

2. 对于后半部分(下标): j 只能来自, 因此.

结论:

也可以写成:

逆变换就是:

3. 按 XOR 卷积

这是最神秘的部分.

目标

定义变换

对于 XOR 卷积, 我们定义:

非常神秘的定义, 使我的大脑直接过载宕机, 推不动柿子了 QwQ .

但这坨东西就是满足的性质.

FWT(A) 的计算方法

直接上结论了.

简记 FWT 为 F , IFWT 为 I.

注意, 除以要按乘以的模逆元来处理.

FWT 代码实现

分治, 一般的写法是递归, 但效率没有循环高, 也不够优雅.

其实分治的本质就是按块长从小到大合并.

• 块长为(mid = 1): 操作,,
• 块长为(mid = 2): 操作,,
• 直到块长为(mid =).

接下来分别给出注释版和适当压行版代码实现.

使用的时候需要定义好 mod 和 inv2.

注释版

struct FWT {
    static void OR(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            // mid << 1 就是块长
            // i 就是块的起始位置
            // 我们取出了块 [i, i + mid * 2)
            // 然后枚举 j, 也就是块内的位置
            // a[i+j] 对应 A0, a[i+j+mid] 对应 A1
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    // FWT(A) = (FWT(A0), FWT(A0) + FWT(A1)), 只需要修改后半段
                    if(type == 1){
                        a[i + j + mid] = (a[i + j + mid] + a[i + j]) % mod;
                    }
                    // IFWT = (IFWT(A0), IFWT(A1) - IFWT(A0)), 只需要修改后半段
                    else{
                        a[i + j + mid] = (a[i + j + mid] - a[i + j] + mod) % mod;
                    }
                }
            }
        }
    }
    static void AND(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    // FWT(A) = (FWT(A0) + FWT(A1), FWT(A1)), 只需要修改前半段
                    if(type == 1){
                        a[i + j] = (a[i + j] + a[i + j + mid]) % mod;
                    }
                    // IFWT(A) = (IFWT(A0) - IFWT(A1), IFWT(A1))
                    else{
                        a[i + j] = (a[i + j] - a[i + j + mid] + mod) % mod;
                    }
                }
            }
        }
    }
    static void XOR(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    ll x = a[i + j];
                    ll y = a[i + j + mid];
                    // FWT(A) = (FWT(A0) + FWT(A1), FWT(A0) - FWT(A1))
                    a[i + j] = (x + y) % mod;
                    a[i + j + mid] = (x - y + mod) % mod;
                    // IFWT(A) 还需要除以 2
                    if(type == -1){
                        a[i + j] = (a[i + j] * inv2) % mod;
                        a[i + j + mid] = (a[i + j + mid] * inv2) % mod;
                    }
                }
            }
        }
    }
};

无注释, 适当压行版

struct FWT {
    static void OR(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    a[i + j + mid] = (a[i + j + mid] + a[i + j] * type + mod) % mod;
                }
            }
        }
    }
    static void AND(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    a[i + j] = (a[i + j] + a[i + j + mid] * type + mod) % mod;
                }
            }
        }
    }
    static void XOR(std::vector<ll>& a, int type){
        int n = a.size();
        for(int mid = 1; mid < n; mid <<= 1){
            for(int i = 0; i < n; i += (mid << 1)){
                for(int j = 0; j < mid; ++j){
                    ll x = a[i + j], y = a[i + j + mid];
                    a[i + j] = (x + y) % mod;
                    a[i + j + mid] = (x - y + mod) % mod;
                    if(type == -1){
                        a[i + j] = (a[i + j] * inv2) % mod;
                        a[i + j + mid] = (a[i + j + mid] * inv2) % mod;
                    }
                }
            }
        }
    }
};

应用

例1. Dominoes

Solution

总共只有种多米诺骨牌, 每个多米诺骨牌有 选/不选 两种状态, 于是总共就有种状态.

将多米诺骨牌的选择集合用二进制数表示, 用表示状态的胜负情况, 胜为, 否则为. 设询问所给的多米诺骨牌集合为, 则答案即为.

可以用 SOS DP 做, 或者直接做一个 FWT OR变换, 也就是, 则即为答案.

那么如何对特定的求呢? 可以把多米诺骨牌看成有向边, 多米诺骨牌集合就是这些边所构成的图. 显然, 获胜条件等价于该图为一条欧拉路径. 用并查集维护连通性, 并记录节点的即可. 欧拉路径的判断条件为: 图连通 且 奇数度数的节点数不超过.

总的计算量约为, 可轻松通过此题.

例2. P6097 【模板】子集卷积

Solution

如果只是, 那就是普通的 FWT OR 了. 但这里多了一个条件 i & j = 0, 也就是和不交.

当感觉棘手的时候, 不妨观察一下条件有什么性质.

• 如果的某位为, 则和的对应位都为.
• 如果的某位为, 则和的对应位必定是一个为, 一个为.

可以推出:.

也就是说, 如果我们想计算, 其中, 则我们只能使用满足的和.

这启发我们把序列按 popcnt 分组.

• 用储存, 前提是, 否则为.
• 用储存, 前提是, 否则为.
• 用储存, 前提是, 否则为.

现在, fa[p] () 可以看作一个长为的序列了, 如果对每一个, 都分别对 fa[p] 和 fb[p] 做 FWT_OR 变换, 也就是:

for p in range(n+1): FWT_OR(fa[p])
for p in range(n+1): FWT_OR(fb[p])

然后, 记,, 根据条件, 对于每一个, 只需要枚举, 就能求出了! 也就是:

对应的伪代码为:

for p in range(n+1): FWT_OR(fa[p])
for p in range(n+1): FWT_OR(fb[p])

# 卷积计算 fc[p][mask]
for mask in range(2 ** n):
    for p in range(n+1):
        fc[p][mask] = 0
        for q in range(p):
            fc[p][mask] += fa[q][mask] * fb[p-q][mask]
# 别忘了逆变换回去
for p in range(n+1): IFWT_OR(fc[p][mask])

最后,就是, 整个做法的时间复杂度为.

例3. UVA13277 Xor Path

Solution

记 d[u] 表示到的边权异或距离, 显然.

记 cnt[x] 表示到的距离为的路径数, 则对于每个, 答案即为:

对 cnt 做 FWT_XOR 即可.

记得的答案要减去. (排除条的路径) , 且最终答案还要除以, 因为题目要求的是无序二元组.

---

拓展阅读:

• OI-Wiki
• znpdco 的博客