1.13 上午-FWT

前言

勿让将来,辜负曾经

劝退了好叭!

正文

知识点

引子

如果屏幕前的你接触过 FFT 或者 NTT 的话,可能会更好理解这一类加速多项式卷积运算的优化算法

好叭,也许你没学过,单着并不影响我们学习 FWT。这个神奇的小算法是用于优化位运算卷积的,可以将 \(O(n^2)\) 的暴力枚举做法转化为 \(O(n \log n)\) 的优秀做法

为了方便表述,我们记多项式 \(A,B,C\) 的第 \(i\) 次项的系数 分别为 \(A_i,B_i,C_i\)

什么是位运算卷积?相信各位都解过高次方程,对于两个多项式 \(A,B\),计算 \(A \times B\) 的过程被称作加法卷积。形式化地,有:

\[C_k =\sum_{k=i+j} A_i \times B_j \]

强调:上面这个东西是加法卷积,不是位运算卷积!(加法卷积是 FFT 和 NTT 用于优化的问题)

卷积形式

注意到这个求和号下方有一个神奇的小式子,叫做 \(k=i+j\),那么我们把加号改成位运算,就是位运算卷积,于是乎就出现了——

  1. 或运算卷积

\[C_k =\sum_{k=i \cup j} A_i \times B_j \]

  1. 与运算卷积

\[C_k =\sum_{k=i \cap j} A_i \times B_j \]

  1. 异或运算卷积

\[C_k =\sum_{k=i \oplus j} A_i \times B_j \]

接下来我们要尝试使用 FWT 对上述卷积进行复杂度优化

写在废话之前

如果没有接触过 FFT,你可以选择略过下面三个自然段

FFT 可以用于优化加法卷积,其原理完全可以类比到 FWT 上。一个多项式存在两种表示方法——系数表示法和点值表示法。

所谓的加法卷积就是系数表示法之间的某种运算,由于多项式长度都是 \(O(n)\) 的,所以朴素暴力的做法就是 \(O(n^2)\) 的。然而,描述一个多项式只需要 \(n\) 个点去表示即可(也就是点值表示法),时间复杂度 \(O(n)\)

所以,FFT 能够优化时间复杂度的原理是:先将多项式的系数表示法转化为点值表示法(DFT),可以借助单位根加速(FFT),然后点值表示法进行运算,最后再将点值表示法转回系数表示法(IFFT)

FWT 也可以效仿 FFT,给多项式做一些变换(FWT),使得他们之间的卷积运算不再消耗大量的时间复杂度,然后再变换回去(IFWT)。具体地,我们希望完成如下操作:

\[A \to FWT[A],B \to FWT[b] \newline FWT[C] = FWT[A] \times FWT[B] \newline FWT[C] \to C \]

针对于不同的位运算,FWT 变换有不同的构造方式。云落能力有限,无法给出构造的详细思维过程,只能给出构造的结论以及其正确性证明

FWT_or(FWT 加速或运算卷积)

构造方式:\(FWT[P]_i = \sum_{j \cup i = i} P_j\)

说人话就是,对于多项式 \(FWT[P]\) 的第 \(i\) 项,它等于原多项式 \(P\) 中所有 \(P_j\) 的和,其中 \(j\) 满足 j|i=i(写成代码形式兴许会好理解些?)

先不考虑如何 \(O(n \log n)\) 解决这个变换问题,先考虑该构造的正确性。简单推一下式子:

\[\begin{aligned} FWT[A]_i \times FWT[B]_i &= \Big( \sum_{j \cup i = i} A_j \Big) \times \Big( \sum_{k \cup i = i} B_k \Big) \newline &= \sum_{j \cup i = i} \sum_{k \cup i = i} A_j \times B_k \newline &= \sum_{(j \cup k) \cup i = i} A_j \times B_k \newline &= FWT[C]_i \nonumber \end{aligned} \]

卷积的部分搞定了,那么该怎么在优秀的复杂度中实现 \(P \to FWT[P]\) 捏?

考虑分治

位运算有一些非常好的性质,叫做按位运算(废话)。也就是说,我们可以尝试,顺次取出最高位,将一个规模较大的问题划分成规模较小的问题。然后噼里啪啦的就一顿分治就 \(O(n \log n)\) 力!

这么干讲不太好懂,还是湿讲比较好

钦定 \(FWT[P]\)\(2^n\) 项(长度不足可以用系数 \(0\) 补齐),那么将这 \(2^n\) 项分为 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\)\(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\) 两部分,进行一些分类讨论……

为了方便表述,我们将原序列 \(P\) 也效仿上述划分方式,并记前半部分为 \(P_0\),后半部分为 \(P_1\)

对于 \(FWT[P]\) 来说,考虑前半部分的贡献来源。由于或运算的神奇性质,显然不可能由 \(P_1\) 贡献,也就只能由 \(P_0\) 贡献。因此,\((FWT[P]_0,...,FWT[P]_{2^{n-1}-1}) = FWT[P_0]\)。同理,考虑后半部分的贡献。很显然无论这一项是什么,他都会给 \(FWT[P]\) 的后半部分造成贡献。因此,\(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1} = FWT[P_0] + FWT[P_1]\)

结合上述分析,我们发现:在或运算卷积下,\(FWT[P]\) 总是满足:

\[FWT[P]=\text{merge}(FWT[P_0],FWT[P_1]+FWT[P_0]) \]

其中 merge 表示多项式的拼接,(可以类比字符串的拼接哈!)

从代码实现上来说,递归分治的写法肯定没问题啦,就是常数有点大,所以建议写成循环的形式,外面套一层长度枚举即可

正变换 FWT 说的差不多了,接下来该说逆变换 IFWT 力!

现在相当于已知 \(FWT[P]\),求 \(IFWT[FWT[P]]\)(也就是 \(P\))。还是回到刚才多项式拆半的过程,注意到贡献给 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\) 的只有 \(P_0\),所以反过来也是一样的。同理,贡献给 \(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\) 的既有 \(P_0\) 又有 \(P_1\),那么求 \(P_1\) 就直接把 \(P_0\) 的部分刨去即可

形式化地,变换如下:

\[IFWT[P'] = \text{merge}(IFWT[{P'}_0],IFWT[{P'}_1]-IFWT[{P'}_0]) \]

从代码实现上来说,FWT 和 IFWT 的结构极为相似,甚至可以说只差一个正负号,所以可以传入一个参数 type 来表示进行的是正变换还是逆变换,是一种节省码量的小技巧

然后就没有然后了,进度条 \(25 \% / 100 \%\)

FWT_and(FWT 加速与运算卷积)

有了上面或运算的铺垫,与运算讲起来可能更舒服一些

构造方式:\(FWT[P]_i = \sum_{j \cap i = i} P_j\)

正确性证明:

\[\begin{aligned} FWT[A]_i \times FWT[B]_i &= \Big( \sum_{j \cap i = i} A_j \Big) \times \Big( \sum_{k \cap i = i} B_k \Big) \newline &= \sum_{j \cap i = i} \sum_{k \cap i = i} A_j \times B_k \newline &= \sum_{(j \cap k) \cap i = i} A_j \times B_k \newline &= FWT[C]_i \nonumber \end{aligned} \]

正变换:

\[FWT[P]=\text{merge}(FWT[P_0]+FWT[P_1],FWT[P_1]) \]

逆变换:

\[IFWT[P'] = \text{merge}(IFWT[{P'}_0]-IFWT[{P'}_1],IFWT[{P'}_1]) \]

与运算和或运算换汤不换药,自己手推吧(其实是云落敲不动了,好累 qaq)!

然后就没有然后了,进度条 \(50 \% / 100 \%\)

FWT_xor(FWT 加速异或运算卷积)

这个还是要好好说一下的!

为了方便表述,定义一种运算 \(\circ\),其运算法则可以被表述为 \(x \circ y = \operatorname{popcount}(x \cap y) \bmod 2\)

说人话就是,\(x,y\) 二进制按位与后的 “\(1\)” 的个数的奇偶性

该运算满足

  1. 交换律:\(x \circ y = y \circ x\)

  2. 结合律:\((x \circ y) \circ z = x \circ (y \circ z)\)

  3. 分配律:\(x \circ (y \oplus z) = (x \circ y) \oplus (x \circ z)\)

请注意分配律中间的符号表示异或运算

构造方式:\(FWT[P]_i = \sum_{i \circ j = 0} P_j - \sum_{i \circ j = 1} P_j\)

正确性证明也很抽象,诸位凑合看,要是对公式过敏的话就把这个构造方式背下来

\[\begin{aligned} FWT[A]_i \times FWT[B]_i &= \Big( \sum_{i \circ j = 0} A_j - \sum_{i \circ j = 1} A_j \Big) \times \Big( \sum_{i \circ k = 0} B_k - \sum_{i \circ j = 1} B_k \Big) \newline &= \Big( \sum_{i \circ j = 0} A_j \times \sum_{i \circ k = 0} B_k + \sum_{i \circ j = 1} A_j \times \sum_{i \circ j = 1} B_k \Big) - \Big( \sum_{i \circ j = 1} A_j \times \sum_{i \circ k = 0} B_k + \sum_{i \circ j = 0} A_j \times \sum_{i \circ j = 1} B_k \Big) \newline &= \sum_{i \circ (j \oplus k) = 0} A_j \times B_k - \sum_{i \circ (j \oplus k) = 1} A_j \times B_k \newline &= FWT[C]_i \nonumber \end{aligned} \]

对于异或运算的 FWT 变换,仍旧是使用分治——对多项式 \(P\) 进行拆半操作

考虑 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\) 的贡献来源,发现不论是 \(P_0\) 还是 \(P_1\) 都可以给它造成贡献(建议自己手模一下,建议打 “\(\circ\)” 运算的真值表)。同理,对于 \(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\),发现 \(P_0\) 可能会多余贡献,且多出来那部分恰好是 \(P_1\)(因为 \(0 \cap 1 = 0,1 \cap 1 = 1\)),所以要在 \(P_0\) 中刨去 \(P_1\) 给的那一部分贡献

综上 ,正变换递归式如下:

\[FWT[P] = merge(FWT[P_0] + FWT[P_1], FWT[P_0] - FWT[P_1]) \]

逆变换的式子写起来也不难。在若干年以前,即便是我们并不会使用二元一次方程组,但是我们依旧可以解决这个问题。

根据我们小学二年级就学过的…… ——毕导

已知两数之和,以及两数之差,求两个数分别是——?

钦定诸位都会上面这个问题,那么逆变换的形式可以写出来了,形如:

逆变换易得:

\[IFWT[P'] = merge(\frac{IFWT[{P'}_0] + IFWT[{P'}_1]}{2}, \frac{IFWT[{P'}_0] - IFWT[{P'}_1]}{2}) \]

如果你听了云落的歪门邪道码量优化小技巧,要注意逆变换的时候 type 不要传入 \(-1\),而应传入 \(\frac{1}{2}\)

然后就没有然后了,进度条 \(75 \% / 100 \%\)

子集卷积

相比于或运算卷积,多了一条限制,形如 \(i \cap j = 0\)。做一步转化,结合另一条限制 \(i \cup j = k\),显然有 \(|i| + |j| = |k|\),也就是转化为了集合大小上的关系。所以直接把原来的多项式系数按照集合大小进行分类,对于相同集合大小的多项式系数进行 FWT 变换,最后多项式 \(A,B\) 合并的时候暴力加法卷积即可

建议结合代码食用

然后就没有然后了,进度条 \(100 \% / 100 \%\)

我们终会重逢

一题一解

T1 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)(P4717)

链接

本来想直接放题解代码的,谁复制粘贴谁就习题棕名大礼包,想想算了,还是放自己的代码哈!

点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=(1<<17)+5,mod=998244353,inv=499122177;
int n,A[maxn],B[maxn];
int a[maxn],b[maxn],c[maxn];
inline void init(){
    memcpy(a,A,sizeof(A));
    memcpy(b,B,sizeof(B));
    memset(c,0,sizeof(c));
    return;
}
inline void mul(){
    for(int i=0;i<n;i++){
        c[i]=a[i]*b[i]%mod;
    }
    return;
}
inline void output(){
    for(int i=0;i<n;i++){
        cout<<c[i]<<' ';
    }
    cout<<endl;
    return;
}
inline void FWT_or(int *P,int type){
    for(int i=1;i<n;i<<=1){
        int p=i<<1;
        for(int j=0;j<n;j+=p){
            for(int k=0;k<i;k++){
                P[i+j+k]=(P[i+j+k]+(P[j+k]*type%mod+mod)%mod)%mod;
            }
        }
    }
    return;
}
inline void FWT_and(int *P,int type){
    for(int i=1;i<n;i<<=1){
        int p=i<<1;
        for(int j=0;j<n;j+=p){
            for(int k=0;k<i;k++){
                P[j+k]=(P[j+k]+(P[i+j+k]*type%mod+mod)%mod)%mod;
            }
        }
    }
    return;
}
inline void FWT_xor(int *P,int type){
    for(int i=1;i<n;i<<=1){
        int p=i<<1;
        for(int j=0;j<n;j+=p){
            for(int k=0;k<i;k++){
                int X=P[j+k],Y=P[i+j+k];
                P[j+k]=(X+Y)%mod;
                P[i+j+k]=(X-Y+mod)%mod;
                P[j+k]=(P[j+k]*type)%mod;
                P[i+j+k]=(P[i+j+k]*type)%mod;
            }
        }
    }
    return;
}
inline void OR(){
    init();
    FWT_or(a,1);
    FWT_or(b,1);
    mul();
    FWT_or(c,-1);
    output();
    return;
}
inline void AND(){
    init();
    FWT_and(a,1);
    FWT_and(b,1);
    mul();
    FWT_and(c,-1);
    output();
    return;
}
inline void XOR(){
    init();
    FWT_xor(a,1);
    FWT_xor(b,1);
    mul();
    FWT_xor(c,inv);
    output();
    return;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin>>n;
    n=(1<<n);
    for(int i=0;i<n;i++){
        cin>>A[i];
    }
    for(int i=0;i<n;i++){
        cin>>B[i];
    }
    OR();
    AND();
    XOR();
    return 0;
}

T2 【模板】子集卷积(P6097)

链接

板,不想多说,现在这个东西不卡常数咯!

温馨提示:模数是 \(10^9 + 9\) 而不是 \(10^9 +7\)

点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=22,maxn=(1<<20)+2,mod=1e9+9;
int n,a[MAXN][maxn],b[MAXN][maxn];
int len,c[MAXN][maxn];
inline int lowbit(int x){
    return x&(-x);
}
inline int cal(int x){
    int res=0;
    while(x){
        res++;
        x-=lowbit(x);
    }
    return res;
}
inline void FWT_or(int *P,int type){
    for(int i=1;i<n;i<<=1){
        int p=i<<1;
        for(int j=0;j<n;j+=p){
            for(int k=0;k<i;k++){
                P[i+j+k]=(P[i+j+k]+(P[j+k]*type%mod+mod)%mod)%mod;
            }
        }
    }
    return;
}
inline void Juan(){
    for(int i=0;i<=len;i++){
        for(int j=0;j<=i;j++){
            for(int k=0;k<n;k++){
                c[i][k]=(c[i][k]+a[j][k]*b[i-j][k]%mod)%mod;
            }
        }
    }
    return;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin>>n;
    len=n;
    n=(1<<n);
    for(int i=0;i<n;i++){
        cin>>a[cal(i)][i];
    }
    for(int i=0;i<n;i++){
        cin>>b[cal(i)][i];
    }
    for(int i=0;i<=len;i++){
        FWT_or(a[i],1);
        FWT_or(b[i],1);
    }
    Juan();
    for(int i=0;i<=len;i++){
        FWT_or(c[i],-1);
    }
    for(int i=0;i<n;i++){
        cout<<c[cal(i)][i]<<' ';
    }    
    return 0;
}

T3 卡牌(P8292)

链接

不大会 FWT 的做法……只会大力容斥

如果你们还对寿司晚宴有印象的话,显然有一个比较经典的 trick,形如根号分治(云落也不知道这么说是否准确)

对于 \(\ge 43\) 的质因数,最多出现一次。所以小质数状态压缩,大质数特殊处理

计算方案,考虑容斥

显然有 \(2^n\) 减去至少不含一个小质数的,加上至少不含两个小质数……显然是对于小质数来说,正确性和复杂度都有保证

现在把大质数的贡献加进来,注意到每次询问最多需要一个大质数。也就是对于每个大质数,按需求给答案 \(\times 2^{cnt_i}\) 或者 \(\times (2^{cnt_i} - 1)\)

剩下的就交给容斥即可

点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=1e6+5,maxv=2e3+5,maxp=13,mod=998244353;
int n,m,cnt[maxv];
vector<int> que;
bool isprime[maxv];
int prime[maxv],rev[maxn],tot;
vector<int> fac[maxn];
int p[maxv];
int g[(1<<maxp)+5],f[(1<<maxp)+5][305];
inline void euler(){
    memset(isprime,true,sizeof(isprime));
    isprime[1]=false;
    for(int i=2;i<maxv;i++){
        if(isprime[i]){
            prime[++tot]=i;
            rev[i]=tot;
		}
        for(int j=1;j<=tot&&i*prime[j]<maxv;j++){
            isprime[i*prime[j]]=false;
            if(i%prime[j]==0){
                break;
			}
        }
    }
    return;
}
inline int qpow(int a,int b){
    int res=1;
    while(b){
        if(b&1){
            res=res*a%mod;
        }
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin>>n;
    for(int i=1;i<=n;i++){
        int x;
        cin>>x;
        cnt[x]++;
    }
    euler();
    for(int i=2;i<maxv;i++){
        for(int j=1;j<=tot;j++){
            if(i%prime[j]==0){
                if(j<=13){
                    p[i]|=(1<<j-1);
                }
                fac[i].push_back(j);
            }
        }
    }
    for(int i=0;i<(1<<maxp);i++){
        for(int j=2;j<maxv;j++){
            if(p[j]&i){
                continue;
            }
            int x=fac[j][fac[j].size()-1];
            f[i][x]=(f[i][x]+cnt[j])%mod;
            g[i]=(g[i]+cnt[j])%mod;
        }
        g[i]=(g[i]+cnt[1])%mod;
    }
    cin>>m;
    while(m--){
        que.clear();
        int c;
        cin>>c;
        int state=0;
        for(int i=1;i<=c;i++){
            int x;
            cin>>x;
            que.push_back(x);
            if(rev[x]<=13){
                state|=(1<<rev[x]-1);
            }
        }
        sort(que.begin(),que.end());
        int res=0;
        for(int i=0;i<(1<<13);i++){
            if((i|state)!=state){
                continue;
            }
            int tol=g[i],ans=1;
            for(int j=0;j<c;j++){
                if(que[j]<=41){
                    continue;
                }
                int x=f[i][rev[que[j]]];
                ans=ans*(qpow(2,x)-1)%mod;
                tol=(tol-x+mod)%mod;
            }
            ans=(ans*qpow(2,tol))%mod;
            int sum=0;
            for(int j=0;j<=13;j++){
                sum=(sum+(i>>j&1))%mod;
            }
            if(sum&1){
                res=(res-ans+mod)%mod;
            }else{
                res=(res+ans)%mod;
            }
        }
        cout<<res<<endl;
    }
    return 0;
}

后记

T3 的 FWT 做法留坑后填(主要是云落只能理解 T3 第一篇 luogu 题解的做法,其它的太抽象了)

完结撒花!

posted @ 2025-03-09 18:20  sunxuhetai  阅读(34)  评论(0)    收藏  举报