快速沃尔什变换 FWT
给定两个数组 $A,B$，求出 $A,B$ 经过下列变换之后得到的 $C$ 数组，快速计算。
$$
C_i = \sum_{j \oplus k=i} A_jB_k
$$
容易发现，快速沃尔什变换就是对下标进行卷积，当 $\oplus$ 为加法的时候就是 FFT/NTT。
其中 $\oplus$ 为 $\operatorname{or},\operatorname{and},\operatorname{xor}$ 中的任意一个，它们分别被称作：或卷积，与卷积，异或卷积。
或卷积
求解或卷积，即为求解：
$$
C_i = \sum_{j | k=i} A_jB_k
$$
首先如果我们设 $A0_i = \sum_{j \in i} A_j,B0_i=\sum_{j \in i} B_j$，那么我们有：
$$
\begin{aligned}
C0i &= \sum_{j \in i} C_j \\
&= \sum_{k|l \in i} A_kB_l \\
&= \sum_{k\in i,l \in i} A_kB_l \\
&= \sum_{k\in i} A_k \sum_{l \in i} B_l \\
&= A0_i B0_i \\
\end{aligned}
$$
而 $A \to A0,B\to B0,C0\to C$ 都很好算，高维前缀和即可，所以时间复杂度就被优化成了 $O(n2^n)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| inline void gor(ll *a,ll *b,ll *c,ll n){
for(ll i=0;i<n;i++) fa[i]=a[i],fb[i]=b[i];
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) (((fa[k+i]+=fa[k])>=mod)&&(fa[k+i]-=mod)),(((fb[k+i]+=fb[k])>=mod)&&(fb[k+i]-=mod));
}
}
for(ll i=0;i<n;i++) c[i]=fa[i]*fb[i]%mod;
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) ((c[i+k]-=c[k])<0&&(c[i+k]+=mod));
}
}
for(ll i=0;i<n;i++) c[i]=(c[i]%mod+mod)%mod;
}
|
与卷积
求解与卷积，即为求解：
$$
C_i = \sum_{j \And k=i} A_jB_k
$$
与或卷积一样的，只不过全部 $j \in i$ 换成了 $i \in j$，公式如下：
首先如果我们设 $ A0_i = \sum_{i \in j} A_j , B0_i=\sum_{i \in j} B_j $，那么我们有：
$$
\begin{aligned}
C0i &= \sum_{i \in j} C_j \\
&= \sum_{i \in (k \And l)} A_kB_l \\
&= \sum_{i\in k,i \in l} A_kB_l \\
&= \sum_{i\in k} A_k \sum_{i \in l} B_l \\
&= A0_i B0_i \
\end{aligned}
$$
求解 C 是一样的，所以就是高维后缀和，时间复杂度 $O(n2^n)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
| inline void gand(ll *a,ll *b,ll *c,ll n){
for(ll i=0;i<n;i++) fa[i]=a[i],fb[i]=b[i];
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) (((fa[k]+=fa[k+i])>=mod)&&(fa[k]-=mod)),(((fb[k]+=fb[k+i])>=mod)&&(fb[k]-=mod));
}
}
for(ll i=0;i<n;i++) c[i]=fa[i]*fb[i]%mod;
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) ((c[k]-=c[i+k])<0&&(c[k]+=mod));
}
}
for(ll i=0;i<n;i++) c[i]=(c[i]%mod+mod)%mod;
}
|
异或卷积
求解异或卷积，即为求解：
$$
C_i = \sum_{j \operatorname{xor} k=i} A_jB_k
$$
这个推导有点复杂，所以我们考虑设代码中的 $f_k=a,f_{k+i}=b$，那么或卷积是 $b \gets a+b$，与卷积是 $a \gets a+b$，异或卷积是 $a \gets a+b,b \gets a-b$。
然后最后我们要逆运算回去，异或与与和或不同，因为我们两者加起来发现多了一倍，而不是用 $a \gets a-b$ 或者 $b \gets b-a$，因此每次用 $a \gets \frac{a+b}2,b \gets \frac{a-b}{2}$ 即可。
具体的推导过程可以用当 $n=2^1$ 的时候特殊情况推导，这个留作复习的时候思考吧。
代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
| inline void gxor(ll *a,ll *b,ll *c,ll n){
for(ll i=0;i<n;i++) fa[i]=a[i],fb[i]=b[i];
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++){
fa[k] += fa[k+i],fa[k+i] = fa[k]-fa[k+i]-fa[k+i];
fa[k] %= mod,fa[k+i] %= mod;
fb[k] += fb[k+i],fb[k+i] = fb[k]-fb[k+i]-fb[k+i];
fb[k] %= mod,fb[k+i] %= mod;
}
}
}
for(ll i=0;i<n;i++) c[i]=fa[i]*fb[i]%mod;
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++){
c[k] += c[k+i],c[k+i] = c[k]-c[k+i]-c[k+i];
c[k] = c[k]*inv2%mod,c[k+i] = c[k+i]*inv2%mod;
}
}
}
for(ll i=0;i<n;i++) c[i]=(c[i]%mod+mod)%mod;
}
|
FWT 的一些性质
上面的三个运算统称为位运算卷积，其实只有或卷积是 FWT，其他的有的叫做莫比乌斯变换等等。
并且由 $A \to A0$ 的过程叫做 FWT 变换，从 $A0\to A$ 的过程叫做 FWT 逆变换（IFWT），其实跟 FFT/NTT 很像。
运算律
我们可以发现，每次从 $ A $ 数组变化到 $ A0 $ 数组都可以看作是一个矩阵乘法线性变换的过程，比如 $ A0_{i}=\sum_{j \in i} A_j $ 可以写作 $ p_{j,i}=[j \in i] $（转移矩阵），因此 FWT 具有结合律和分配律而不具有交换律。
所以如果我们在某个算法中会用到 FWT 变换的话，大部分情况都可以在程序最开始先令 $ A \to A0 $，然后用 $A0$ 进行变换，最后用 $ A0 $ 还原 $ A $ 即可，但是要注意，算法里面必须进行的是线性变换才可以这么处理，例题可以参见 P5406 THUPC2019 找树。
或者我们如果遇到了要求 $ A_i = \sum_{n_0}\sum_{n_1}\sum_{n_2}\dots\sum_{n_k}[n_0|n_1|n_2|\dots|n_k=i]B_{n_0}B_{n_1}\dots B_{n_k} $ 这类题目的话，可以对 $B$ 直接进行 FWT 变换得到 $ B0 $，然后令 $ A0=B0^{k+1} $，最后通过 $ A0 \to A $ 即可，这样的话时间复杂度通常会少一个 $ \log $，是一个很大的优化。（证明可以通过上面的分配律或者把式子写出来展开即可）
分位考虑
如果我们要求对于二进制下的第 $pos$ 位采用了规定的 或/与/异或 的运算法则，那 FWT 还适用吗？
答案是肯定的，我们只需要在代码的第一层枚举中记录位数，然后位数是什么就运用哪种变换就可以了，比如下面的 FWT 和 IFWT 就是这样的。
顺带解释一下第一层 for 循环是枚举高维前缀和处理的位数，第二层是枚举没有处理的位数，第三层是枚举处理过的位数，可以想象一下 $ \text{second_num}+1\text{（第一层枚举的位）}+\text{third_num} $。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
| inline poly fwt(poly a){
poly fa = a;
for(ll i=1,pos=0;i<(1<<w);i*=2,pos++){
for(ll j=0;j<(1<<w);j+=2*i){
for(ll k=j;k<j+i;k++){
if(s[pos]=='|'){
((fa.a[k+i]+=fa.a[k])>=mod)&&(fa.a[k+i]-=mod);
}
if(s[pos]=='&'){
((fa.a[k]+=fa.a[k+i])>=mod)&&(fa.a[k]-=mod);
}
if(s[pos]=='^'){
((fa.a[k]+=fa.a[k+i])>=mod)&&(fa.a[k]-=mod),fa.a[k+i] = fa.a[k]-fa.a[k+i]*2;
}
}
}
}
for(ll i=0;i<(1<<w);i++) fa.a[i]=(fa.a[i]%mod+mod)%mod;
return fa;
}
inline poly ifwt(poly a){
poly c = a;
for(ll i=1,pos=0;i<(1<<w);i*=2,pos++){
for(ll j=0;j<(1<<w);j+=2*i){
for(ll k=j;k<j+i;k++){
if(s[pos]=='|'){
((c.a[k+i]-=c.a[k])<0)&&(c.a[k+i]+=mod);
}
if(s[pos]=='&'){
((c.a[k]-=c.a[k+i])<0)&&(c.a[k]+=mod);
}
if(s[pos]=='^'){
c.a[k] += c.a[k+i],c.a[k+i] = c.a[k]-c.a[k+i]*2;
c.a[k] = c.a[k]*inv2,c.a[k+i] = c.a[k+i]*inv2;
}
}
}
}
for(ll i=0;i<(1<<w);i++) c.a[i]=(c.a[i]%mod+mod)%mod;
return c;
}
|
例题依旧是上面那道 THUPC 的题目。
子集卷积及其相关操作
子集卷积
给出数组 $A,B$ 求数组 $C$，其中：
$$
C_i = \sum_{j|k=i,j \And k=0} A_jB_k
$$
首先，我们可以发现如果没有 $j \And k=0$ 这个限制就是 FWT 或卷积，但是有一个限制就是 $j \And k=0$，考虑转化这个限制变成 $\operatorname{popcount}(j)+\operatorname{popcount}(k)=\operatorname{popcount}(i)$，于是如果我们给 $A$ 扩展一维，给 $B$ 扩展一维，并令 $A_{i}=A0_{\operatorname{popcount(i)},i}$，$B$ 同理，那么我们得到：
$$
C0_{\operatorname{popcount}(i),i} = \sum_{j|k=i} A0_{\operatorname{popcount(j)},j}B0_{\operatorname{popcount(k)},k}
$$
那么 $j$ 和 $k$ 的 $\And$ 的限制就没有了，又因为 FWT 对加法具有分配律，所以我们可以先对 $A0$ 和 $B0$ 的 $pos$ 个（$pos$ 是二进制位数）维度做 FWT 变换，然后令：
$$
C1_{b,i} = \sum_{b=0}^a A1_{a,i}B1_{b-a,i}
$$
最后对 $C1$ 做 IFWT 就得到了 $C0$ 的值，进而可以得到 $C$ 的值。（其实应该是枚举了一个 $b$ 就算一次累加答案，但是因为 FWT 的分配律，所以可以枚举所有 $b$ 把结果加在一起之后再 IFWT）这个时候我们发现有些值虽然它二进制下 $1$ 的个数不等于当前维度，但是它的位置上依然有值，这是正常的（这个值没用），如果这个值对后面的计算有影响，清零即可。
总的时间复杂度是 $O(2^nn^2)$，代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
| inline void gor1(ll *a,ll n){
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) (((a[k+i]+=a[k])>=mod)&&(a[k+i]-=mod));
}
}
}
inline void gor2(ll *a,ll n){
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++) ((a[i+k]-=a[k])<0&&(a[i+k]+=mod));
}
}
for(ll i=0;i<n;i++) a[i]=(a[i]%mod+mod)%mod;
}
inline void times(ll *a,ll *b,ll *c,ll n){
ll limit = 0;
for(ll i=0;i<n;i++) f0[__builtin_popcount(i)][i]=a[i],f1[__builtin_popcount(i)][i]=b[i],limit=max(limit,(ll)__builtin_popcount(i));
for(ll i=0;i<=limit;i++) gor1(f0[i],n),gor1(f1[i],n);
for(ll i=0;i<=limit;i++) for(ll j=0;i+j<=limit;j++) for(ll k=0;k<=n;k++) f2[i+j][k]=(f2[i+j][k]+f0[i][k]*f1[j][k])%mod;
for(ll i=0;i<=limit;i++) gor2(f2[i],n);
for(ll i=0;i<n;i++) c[i]=f2[__builtin_popcount(i)][i];
for(ll i=0;i<=limit;i++) for(ll j=0;j<=n;j++) f1[i][j]=f2[i][j]=f0[i][j]=0;
}
|
子集卷积的逆运算
如题，这一个小节主要探讨的是子集卷积除法的一些问题。
给出数组 $A,C$，求数组 $B$，这三个数组满足下面的条件：
$$
C_i = \sum_{j|k=i,j \And k=0} A_jB_k
$$
根据上面子集卷积的定义，我们得到了 $A1,B1$ 之后，通过：
$$
C1_{b,i} = \sum_{b=0}^a A1_{a,i}B1_{b-a,i}
$$
得到了 $C1$，现在同理，我们得到了 $A1,C1$，也可以通过这个式子从 $0 \sim n$ 递推得到 $B1$ 的值，很简单，这里有一道例题，大概就是把题目的 dp 转化为知道上面式子的 $C$ 和 $A$ 求 $B$ 的问题进而求解，时间复杂度依然是 $O(n^22^n)$。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| inline void div(ll *a,ll *b,ll *c,ll n){
ll limit = 0;
for(ll i=0;i<n;i++) f0[__builtin_popcount(i)][i]=a[i],f1[__builtin_popcount(i)][i]=b[i],limit=max(limit,(ll)__builtin_popcount(i));
for(ll i=0;i<=limit;i++) gor1(f0[i],n),gor1(f1[i],n);
for(ll i=0;i<n;i++){
ll inv = qmi(f1[0][i],mod-2,mod);
f2[0][i]=f0[0][i]*inv%mod;
for(ll j=1;j<=limit;j++){
f2[j][i]=f0[j][i];
for(ll k=0;k<j;k++) f2[j][i]=(f2[j][i]-f2[k][i]*f1[j-k][i]%mod+mod)%mod;
f2[j][i]=f2[j][i]*inv%mod;
}
}
for(ll i=0;i<=limit;i++) gor2(f2[i],n);
for(ll i=0;i<n;i++) c[i]=f2[__builtin_popcount(i)][i];
for(ll i=0;i<=limit;i++) for(ll j=0;j<=n;j++) f1[i][j]=f2[i][j]=f0[i][j]=0;
}
|
子集卷积的其他运算
同理，还有开方等等其他操作，这里简单介绍一下开方，而 $\ln,\exp$ 的问题目前不需要掌握并且需要较高的数学知识因此此处不展开。
开方，也就是：
$$
C_i = \sum_{j|k=i,j \And k=0} A_jA_k
$$
知道 $C$ 求 $A$，还是先把 $C1$ 算出来，然后得到了：
$$
C1_{b,i} = \sum_{b=0}^a A1_{a,i}A1_{b-a,i}
$$
然后我们从 $0$ 开始，$0$ 需要通过二次剩余求出来，然后其他项就可以递推了，也很简单，时间复杂度依然是简单的 $O(n^22^n)$，二次剩余基本不占时间。
动态子集卷积
这个通常用于 dp 优化里面，也就是我们知道了 $G,Q$ 两个数组，要想通过下面的方式得到 $F$ 数组：
$$
F_S Q_S = \sum_{T\in S,T \ne 0} G_TF_{S/T}
$$
$T$ 一般不可能是空集，如果强制要求 $T$ 包含某一个元素的话依然可以做。
首先我们发现 $\operatorname{popcount}(S/T)<\operatorname{popcount}(S)$，于是我们可以第 $i$ 次处理出所有 $F_S$ 满足 $\operatorname{popcount}(S)=i$，因为 $G$ 和 $Q$ 是确定的，于是我们就可以通过所有 $\operatorname{popcount}(S)<i$ 的 $F_S$ 得到这一层的 $F_S$，同样的，我们把 $\operatorname{popcount}$ 放到维度里面，假如我们得到了 $\operatorname{popcount}(S)<i$ 的 $F_S$，我们考虑如何求出 $\operatorname{popcount}(T)=i$ 的 $F_T$（像普通子集卷积那样给 $F$ 增加一维）：
$$
F_{\operatorname{popcount}(T),T}Q_T=\sum_{i=1}^{\operatorname{popcount}(T)} \sum_{j|k=T} G_{i,j}F_{\operatorname{popcount}(T)-i,k}
$$
这样的话，求出来的答案一定是正确的，因为不合法的状态都会累加到不合法的 $F$ 上，完成整个过程之后把 $F_{i,j}(\operatorname{popcount}(j) \ne i)$ 置为 $0$ 即可。
容易发现，每次后面两个和式中的值都是确定的，因此我们只需要预处理出 $G$ 每一层 FWT 后的结果 $G0$ 和 $F$ 每一层 FWT 后的结果 $F0$ 即可。
然后就可以得到：
$$
F0_{\operatorname{popcount}(T),T}Q_T=\sum_{i=1}^{\operatorname{popcount}(T)} G0_{i,T}F0_{\operatorname{popcount}(T)-i,T}
$$
于是再 IFWT 回去即可，时间复杂度为 $O(n^22^n)$，例题是 P4221 WC2018 州区划分，参考代码如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
| #include<bits/stdc++.h>
#define ll int
#define mod 998244353
#define N 22
namespace IO{
inline char nc(){
static char buf[1000000],*p=buf,*q=buf;
return p==q&&(q=(p=buf)+fread(buf,1,1000000,stdin),p==q)?EOF:*p++;
}
inline ll read(){
ll res = 0,w = 1;
char c = nc();
while(c<'0'||c>'9')w=(c=='-'?-1:w),c=nc();
while(c<='9'&&c>='0')res=res*10+c-'0',c=nc();
return res*w;
}
char obuf[1<<21],*p34=obuf;
inline void pc(char c){
p34-obuf<=(1<<20)?(*p34++=c):(fwrite(obuf,p34-obuf,1,stdout),p34=obuf,*p34++=c);
}
inline void write(ll x){
if(x>9) write(x/10);
pc(x%10+'0');
}
}
using namespace IO;
using namespace std;
ll n,m,p,i,j,k,x,y,et,la[N],to[N<<5],ne[N<<5],w[N],du[N],fath[N],pos,f[N][1<<N-1],g[N][1<<N-1],cntt[1<<N-1],inv[1<<N-1],sum[1<<N-1];
inline ll gf(ll x){return x==fath[x]?x:fath[x]=gf(fath[x]);}
inline ll qmi(ll a,ll b,ll p){
ll res = 1%p,t = a;
while(b){
if(b&1) res=1ll*res*t%p;
t=1ll*t*t%p;
b>>=1;
}
return res;
}
inline void fwt(ll *a,ll n){
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++){
(a[k+i]+=a[k])>=mod&&(a[k+i]-=mod);
}
}
}
}
inline void ifwt(ll *a,ll n){
for(ll i=1;i<n;i*=2){
for(ll j=0;j<n;j+=2*i){
for(ll k=j;k<j+i;k++){
(a[k+i]-=a[k])<0&&(a[k+i]+=mod);
}
}
}
}
int main(){
n=read(),m=read(),p=read();
for(i=1;i<=m;i++){
x=read(),y=read();
et++,ne[et]=la[x],la[x]=et,to[et]=y;
et++,ne[et]=la[y],la[y]=et,to[et]=x;
}
for(i=1;i<=n;i++) w[i]=read();
for(i=1;i<(1<<n);i++){
ll cnt=1;
for(j=1;j<=n;j++) fath[j]=j,du[j]=0;
for(j=0;j<n;j++){
if((i>>j)&1){
pos=j,cntt[i]=cntt[i-(1<<j)]+1,sum[i]=sum[i-(1<<j)]+w[j+1];
for(k=la[j+1];k;k=ne[k]){
if((i>>(to[k]-1))&1) fath[gf(j+1)]=gf(to[k]),du[j+1]++;
}
}
}
for(j=0;j<n;j++) if(du[j+1]%2==1||(((i>>j)&1)&&gf(pos+1)!=gf(j+1))) cnt=0;
sum[i] %= mod,inv[i] = qmi(qmi(sum[i],p,mod),mod-2,mod);
if(!cnt) g[cntt[i]][i] = (g[cntt[i]][i]+qmi(sum[i],p,mod))%mod;
}
for(i=0;i<=n;i++) fwt(g[i],1<<n);
f[0][0]=1,fwt(f[0],1<<n);
for(i=1;i<=n;i++){
for(j=1;j<=i;j++) for(k=0;k<(1<<n);k++) f[i][k]=(f[i][k]+1ll*g[j][k]*f[i-j][k])%mod;
ifwt(f[i],1<<n);
for(j=0;j<(1<<n);j++){
if(cntt[j]==i) f[i][j]=1ll*f[i][j]*inv[j]%mod;
else f[i][j]=0;
}
if(i<n) fwt(f[i],1<<n);
}
cout<<f[n][(1<<n)-1]<<endl;
return 0;
}
|