SOSDP（子集DP）、高维前缀和、FMT（快速莫比乌斯变换）、FWT（快速沃尔什变换）学习笔记
这里是SOSDP（子集DP）、高维前缀和、FMT（快速莫比乌斯变换）、FWT（快速沃尔什变换）学习笔记。
在接下来的讲解中，我们将会发现前三个东西是本质相同的，而最后一个东西是其扩展。
本博客在半年前就已经有计划去写了，只不过当时题目较少没能如愿。现在攒了几道题，开始动笔。
引入
我们首先来看一道最板子的例题：
全集 \(\mathbb{U}\) 中共有 \(n\) 个元素。对于全部 \(\mathbb{S}\subseteq\mathbb{U}\)，定义了一函数 \(f(\mathbb S)\)。现在，再定义另一函数 \(g(\mathbb S)=\sum\limits_{\mathbb T\subseteq\mathbb S}f(\mathbb T)\)。已知全部 \(f(\mathbb S)\) 的值，求全体 \(g(\mathbb S)\) 的值。
数据范围：\(n\leq 22\)。
暴力1：
\(2^n\) 地分开枚举 \(\mathbb{S}\) 和 \(\mathbb T\)，复杂度 \(O(4^n)\)。
暴力2：
使用子集枚举的trick，将复杂度优化到 \(O(3^n)\)。
暴力3：
我们不妨先从维数少一点的情形考虑，比如 \(n=1\)。此时共有两个集合，\(\varnothing\) 和 \(\{x_1\}\)。显然，\(g(\varnothing)=f(\varnothing),g\Big(\{x_1\}\Big)=f\Big(\{x_1\}\Big)+f(\varnothing)\)。
现在再考虑维数稍大一点的情形，比如 \(n=2\)。此时，有
发现其看上去一副可以容斥的样子。于是我们稍加思索，将最后一个式子改成了这样：
再考虑更高维的样式，就有 \(g(\mathbb S)=f(\mathbb S)+\sum\limits_{\mathbb T\subset\mathbb S}(-1)^{|\mathbb S|-|\mathbb{T}|+1}g(\mathbb T)\)。
虽然看上去比我们前两种的暴力更高级了，但是复杂度仍为 \(O(3^n)\)，甚至比暴力2常数更劣。
高维前缀和
上述容斥的思想启发我们，什么情形下我们还要用容斥？
二维前缀和，它与上述二维时情形十分类似。
于是我们发现，我们要求的东西本质上是高维前缀和。
SOSDP
我们发现，前面的算法中，我们没有省掉任何的枚举。这启发我们通过DP的方式来减少枚举量。
于是就有了这种美妙的东西——SOSDP，全称 Sum Over Subsets(SOS) DP，中文名”子集DP“。
顾名思义，其计算的是子集之和，也即我们问题中要求的东西。
我们不妨先定下一个枚举顺序——比如说，先处理所有含物品 \(x_1\) 的集合，再处理所有含物品 \(x_2\) 的集合……
于是我们可以找到如下的DP方式：枚举物品 \(x_i\)，再从小到大枚举所有的集合 \(\mathbb S\)。若 \(x_i\in\mathbb S\)，则 \(g\Big(\mathbb S\setminus\{x_i\}\Big)\rightarrow g(\mathbb S)\)，其中 \(\setminus\) 符号是子集删去符号。
显然，当枚举 \(x_1\) 时，缺失 \(x_1\) 的 \(\mathbb T\) 转移到了 \(\mathbb S\)；枚举 \(x_2\) 时，缺失 \(x_1,x_2\) 中任意多个的 \(\mathbb{T}\) 转移到了 \(\mathbb{S}\)；枚举 \(x_3\) 时……
于是我们发现该DP方式枚举了所有 \(\mathbb{T}\subset\mathbb S\)。再加上 \(f(\mathbb{S})\) 本身，即得到 \(f(\mathbb{S})\)。
可以发现，该算法复杂度仅为 \(O(n2^n)\)。通过SOSDP，我们得以实现了高维前缀和。
以下是高维前缀和模板的代码（虽然就很短一点）：
void SUM(int *a,int n){for(int i=0;i<n;i++)for(int j=0;j<(1<<n);j++)if(j&(1<<i))a[j]+=a[j^(1<<i)];}
高维差分
众所周知，前缀和的逆运算是差分。
虽然高维差分的定义比较困难，但是换一种说法就能很轻易地理解——我们能否找到一种方式，使得已知 \(g\) 数组，求出 \(f\) 数组？
明显，因为我们已经知道 \(g\) 的递推式，故只需将循环全部反向、加号变减号，即可做到高维差分。
但是，如果我们仔细看看的话，就会发现因为一个 \(\mathbb S\) 中某个元素最多出现一次，所以此二循环的枚举顺序无论正反是均可的。于是，从高维前缀和到高维差分，唯一的变化就是加号变减号。
高维后缀和
高维后缀和在下面也有着诸多应用。其定义是 \(h(\mathbb S)=\sum\limits_{\mathbb S\subseteq\mathbb T}f(\mathbb T)\)。但是，通过将所有东西取补集，其就变成了高维前缀和，也就可以类似地处理了。同理也可以处理高维后缀差分。
快速莫比乌斯变换
一定要将其与关于 \(\mu\) 函数的“莫比乌斯反演”分清！
我们考虑已知两个集合上函数 \(f\) 和 \(g\)，要求出 \(h(\mathbb S)=\sum\limits_{\mathbb I\cup\mathbb J=\mathbb S}f(\mathbb{I})g(\mathbb J)\)。因为若把 \(\mathbb{I,J}\) 均看作 \(2^n\) 维的状态的话，其等价于 \(I\text{ or }J=S\)。故这种操作也可以被叫做 \(f\text{ or }g\)。
现在，考虑对 \(f\) 进行高维前缀和操作，得到一个新的函数 \(\hat{f}\)。考虑计算 \(\hat{h}(\mathbb{S})=\hat{f}(\mathbb{S})\times\hat{g}(\mathbb{S})\)。假如我们把 \(\hat{f}(\mathbb S)\) 与 \(\hat{g}(\mathbb S)\) 的定义代入的话，会发现得到的 \(\hat h(\mathbb S)\) 的实际意义也是 \(h\) 的高维前缀和。于是对其作高维差分即可得到 \(h\) 本身。
我们发现，这种 \(\hat f\)，实际上跟FFT中的“DFT”操作特别类似——同样是对一个函数作运算得到另一个函数，然后用运算完的操作进行按位求积，最后再进行反向操作将其变换为原函数。而整套流程下来，实现了卷积的效果，因此这套流程被叫做“\(\text{or}\) 卷积”。在原函数和新函数间建立双射的操作，被称作“快速莫比乌斯变换Fast Mobius Transform(FMT)”。
总结一下。通过对函数进行莫变（高维前缀和），可以得到新函数。对新函数进行按位求积，可以得到答案的新函数。对其进行莫变的反操作（也被称作快速莫比乌斯反演FMI，但是因为与普通的莫反太过相似故我们不用这个称呼，直接叫其“反莫变”）（实际上就是高维差分），得到答案。整一套流程相当于 or 卷积。
同理，我们也可以得到 and 卷积的操作——只不过使用的是高维后缀和和高维后缀差分罢了。
or 卷积与 and 卷积的复杂度都是 \(O(n2^n)\)。但是，如果我们重新令 \(n\) 表示状态总数（序列长度）的话，就会发现复杂度是 \(n\log n\)，同FFT复杂度一致。
到这里我们先打住吧，来看几道例题。
I.CF165E Compatible Numbers
实际上已经在我的任务计划里待了半年没能来写题解了。
我们开一个数组 \(f\)，并且令全集（所有位上全一的数）为 \(U\)。当我们有一个数 \(x\) 时，将 \(f\) 的第 \(x\) 位设作 \(1\)。然后，对其作高维前缀和（但实际上不需要高维前缀和，理论上应该叫“高维前缀或”——拿所有子集的东西取或），之后判断对于每个 \(x\) 其 \(f\) 数组上的位置 \((U\operatorname{xor}x)\) 是否有值即可。
代码：
#include<bits/stdc++.h>
using namespace std;
const int lim=1<<22;
int n,a[1001000],f[lim];
int main(){
scanf("%d",&n),memset(f,-1,sizeof(f));
for(int i=1;i<=n;i++)scanf("%d",&a[i]),f[a[i]]=a[i];
for(int j=0;j<22;j++)for(int i=1;i<lim;i++)if((i&(1<<j))&&f[i^(1<<j)]!=-1)f[i]=f[i^(1<<j)];
for(int i=1;i<=n;i++)printf("%d ",f[(lim-1)^a[i]]);
return 0;
}
II.CF383E Vowels
我们考虑对于所有集合求出正确单词的数量。在一个数组上标注全部单词后，对其作高维前缀和即可。
代码：
#include<bits/stdc++.h>
using namespace std;
const int lim=(1<<24);
int n,f[1<<24],res;
char s[5];
int main(){
scanf("%d",&n);
for(int i=0;i<n;i++)scanf("%s",s),f[(1<<(s[0]-'a'))|(1<<(s[1]-'a'))|(1<<(s[2]-'a'))]++;
for(int j=0;j<24;j++)for(int i=1;i<lim;i++)if(i&(1<<j))f[i]+=f[i^(1<<j)];
for(int i=0;i<lim;i++)res^=(n-f[i])*(n-f[i]);
printf("%d\n",res);
return 0;
}
III.CF449D Jzzhu and Numbers
考虑求出 \(f(x)\) 表示 and 结果有 \(x\) 作为子集的方案数，则对 \(f(x)\) 作高维反向差分即可得到 and 结果恰为 \(x\) 的方案数。则答案即为差分后的 \(f(0)\)。考虑如何求出原本的 \(f(x)\)。显然，\(f(x)\) 需要参与 and 的所有东西都有 \(x\) 作为子集，故我们做一个高维后缀和即可求出有 \(x\) 作为子集的元素数量，再把它做一个以 \(2\) 为底的幂即可得到总方案数，也即 \(f(x)\)。
但需要注意的是，要排除一个数都不选的方案。
代码：
#include<bits/stdc++.h>
using namespace std;
const int lim=1<<20;
const int mod=1e9+7;
int f[lim],n,pov[1001000];
int main(){
scanf("%d",&n),pov[0]=1;
for(int i=1,x;i<=n;i++)scanf("%d",&x),f[x]++,pov[i]=(pov[i-1]<<1)%mod;
for(int j=0;j<20;j++)for(int i=lim-1;i>=0;i--)if(!(i&(1<<j)))(f[i]+=f[i|(1<<j)])%=mod;
for(int i=0;i<lim;i++)f[i]=pov[f[i]]-1;
for(int j=0;j<20;j++)for(int i=lim-1;i>=0;i--)if(!(i&(1<<j)))(f[i]=f[i]+mod-f[i|(1<<j)])%=mod;
printf("%d\n",f[0]);
return 0;
}
IV.CF582E Boolean Function
考虑建出表达式树，并令 \(f(x,y)\) 表示若树上节点 \(x\) 在全部 \(n\) 个函数中的取值是状态 \(y\) 的方案数。于是拼接两个儿子的答案时就必须使用 and 和 or 卷积进行处理。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
char s[510];
int ch[510][2],S,mat[510],rt,f[510][1<<16],n,g[5],h[4],all,A[1<<16],B[1<<16];
stack<int>stk;
int build(int l,int r){
if(l==r)return l;
int x=mat[l]+1;
ch[x][0]=build(l+1,mat[l]-1);
ch[x][1]=build(mat[r]+1,r-1);
return x;
}
void ORFMT(int *a,int tp){for(int i=0;i<n;i++)for(int j=0;j<all;j++)if(j&(1<<i))(a[j]+=(mod+tp*a[j^(1<<i)])%mod)%=mod;}
void OR(int *a,int *b,int *c){
for(int i=0;i<all;i++)A[i]=a[i],B[i]=b[i];
ORFMT(A,1),ORFMT(B,1);
for(int i=0;i<all;i++)A[i]=1ll*A[i]*B[i]%mod;
ORFMT(A,-1);
for(int i=0;i<all;i++)(c[i]+=A[i])%=mod;
}
void ANDFMT(int *a,int tp){for(int i=0;i<n;i++)for(int j=all-1;j>=0;j--)if(j&(1<<i))(a[j^(1<<i)]+=(mod+tp*a[j])%mod)%=mod;}
void AND(int *a,int *b,int *c){
for(int i=0;i<all;i++)A[i]=a[i],B[i]=b[i];
ANDFMT(A,1),ANDFMT(B,1);
for(int i=0;i<all;i++)A[i]=1ll*A[i]*B[i]%mod;
ANDFMT(A,-1);
for(int i=0;i<all;i++)(c[i]+=A[i])%=mod;
}
void dfs(int x){
// printf("%d\n",x);
if(!ch[x][0]&&!ch[x][1]){
if('a'<=s[x]&&s[x]<='d')f[x][h[s[x]-'a']]=1;
if('A'<=s[x]&&s[x]<='D')f[x][g[s[x]-'A']]=1;
if(s[x]=='?')for(int i=0;i<4;i++)f[x][g[i]]++,f[x][h[i]]++;
// printf("%d:",x);for(int i=0;i<all;i++)printf("%d ",f[x][i]);puts("");
return;
}
dfs(ch[x][0]),dfs(ch[x][1]);
if(s[x]=='&'||s[x]=='?')AND(f[ch[x][0]],f[ch[x][1]],f[x]);
if(s[x]=='|'||s[x]=='?')OR(f[ch[x][0]],f[ch[x][1]],f[x]);
// for(int i=0;i<all;i++)printf("%d ",f[x][i]);puts("");
}
int main(){
scanf("%s",s+1),S=strlen(s+1);
for(int i=1;i<=S;i++){
if(s[i]=='(')stk.push(i);
if(s[i]==')')mat[stk.top()]=i,mat[i]=stk.top(),stk.pop();
}
rt=build(1,S);
scanf("%d",&n),all=1<<n;
for(int i=0,x;i<n;i++)for(int j=0;j<5;j++)scanf("%d",&x),g[j]|=x<<i;
for(int i=0;i<4;i++)h[i]=(~g[i])&(all-1);
dfs(rt);
printf("%d\n",f[rt][g[4]]);
return 0;
}
V.CF1326F2 Wise Men (Hard Version)
Observation 1. 一个状态的答案与且仅与其中连续的
1的段的长度构成的可重集有关。举例来说，应有 \(ans(\overline{100110111})=ans(\overline{011101011})\)。
不太严谨的证明可以用因为状态涵盖了所有阶乘，所以每一组分段方法必定一一对应着另一组。
于是，我们可以考虑枚举每一组可重集计算答案。明显，所有可重集数量是 \(18\) 的整数划分，共 \(1582\) 种。如果是 \(1582\times2^{18}\) 的复杂度，四亿出头次计算还是能卡过去的。
现在考虑如何统计一组可重集的解。如果建立一一对应关系的话，不仅要求 \(1\) 的段数及长度吻合，还要求两段 \(1\) 之间必须有至少一个 \(0\) 进行分割，而这是不简单的。但是，假如我们只要求 \(1\) 的段数和长度相同，并不要求 \(0\) 的分割的话，就只需要在求出答案后进行高维差分即可（有一些 \(1\) 被当成 \(0\) 来看了）。
于是我们现在考虑为可重集中每一段的长度对应排列中一段链。假设可重集中长度序列为 \({d_1,d_2,\dots,d_m}\)，而我们找到的一堆链的集合分别是 \(S_1,S_2,\dots,S_m\)，则需要满足 \(d_i=|S_i|,S_i\cap S_j=\varnothing\)。
于是我们考虑设 \(g(S)\) 表示集合 \(S\) 可以连成的链的数量。其可以通过设 \(f(S,x)\) 表示集合为 \(S\) 且终点为 \(x\) 的链的数量在 \(O(2^nn^2)\) 的时间内DP出来。于是一组可重集的答案就是 \(\sum\limits_{\{S\}}\Big[d_i=|S_i|\Big]\Big[S_i\cap S_j=\varnothing\Big]\prod g(S_i)\)。
\(S_i\cap S_j=\varnothing\) 这个限制很像我们接下来学到的子集卷积中的限制（可以往下翻去先学掉）。于是我们扩展 \(g\) 的定义，设 \(g(i,S)\) 表示 \(S\) 集合的链数，并强制只有 \(|S|=i\) 的状态才是合法的。于是只要计算 \(h=\prod\limits_{i=1}^mg(d_m)\) （因为我们将 \(g\) 扩张至二维，所以这里是卷积操作），则该可重集对应的所有状态的答案均是 \(h(\text{fullset})\)。
在计算多个式子的卷积时，可以一次性先全 FMT 掉，按位全部点乘在一起后，再 IFMT。并且，因为我们只要求 \(h(\text{fullset})\) 处的值，最后连 IFMT 都不需要，直接容斥一波带走。因为我们可以在预处理中就 FMT 过，所以按照上述算法，计算可重集的答案只需要 \(O(\text{size}2^n)\) 即可，其中 \(\text{size}\) 是可重集大小。并且，若我们是在搜索的过程中处理的话，复杂度还能被进一步优化到 \(1582\times 2^n\)，按照我们上述结论，是可以通过的。
则总复杂度 \(O\Big(2^n(1582+n^2)\Big)\)。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int n,lim,tp;
ll f[18][1<<18],g[19][1<<18];
char s[20][20];
void ORFMT(ll*a,int tp){for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))a[j]+=tp*a[j^(1<<i)];}
void ANDFMT(ll*a,int tp){}
ll h[19][1<<18],res[1<<18];
map<vector<int>,vector<int> >mp;
vector<int>v;
void dfs(int lef,int bon){
if(!lef){
ll ret=0;for(int i=0;i<lim;i++)if(__builtin_popcount((lim-1)^i)&1)ret-=h[tp][i];else ret+=h[tp][i];
for(auto i:mp[v])res[i]+=ret;return;
}
for(int i=1;i<=min(lef,bon);i++){
v.push_back(i);
for(int j=0;j<lim;j++)h[tp+1][j]=h[tp][j]*g[i][j];
tp++,dfs(lef-i,i),tp--;
v.pop_back();
}
}
int main(){
scanf("%d",&n),lim=1<<n;
for(int i=0;i<n;i++)scanf("%s",s[i]);
for(int i=0;i<n;i++)f[i][1<<i]=1;
for(int j=0;j<lim;j++)for(int i=0;i<n;i++)if(j&(1<<i))for(int k=0;k<n;k++)if(!(j&(1<<k))&&s[i][k]=='1')f[k][j|(1<<k)]+=f[i][j];
// for(int i=0;i<n;i++){for(int j=0;j<lim;j++)printf("%d ",f[i][j]);puts("");}
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)g[__builtin_popcount(j)][j]+=f[i][j];
// for(int i=0;i<lim;i++)printf("%d ",g[__builtin_popcount(i)][i]);puts("");
for(int i=0;i<=n;i++)ORFMT(g[i],1);
for(int i=0;i<(1<<(n-1));i++){
vector<int>v;
int len=1;
for(int j=0;j<n-1;j++)if(i&(1<<j))len++;else v.push_back(len),len=1;
v.push_back(len);
// printf("%d:",i);for(auto j:v)printf("%d ",j);puts("");
sort(v.rbegin(),v.rend()),mp[v].push_back(i);
}
h[0][0]=1,ORFMT(h[0],1),dfs(n,n);
// for(int i=0;i<(1<<(n-1));i++)printf("%d ",res[i]);puts("");
for(int i=0;i<n-1;i++)for(int j=0;j<(1<<(n-1));j++)if(j&(1<<i))res[j^(1<<i)]-=res[j];
for(int i=0;i<(1<<(n-1));i++)printf("%lld ",res[i]);
return 0;
}
快速沃尔什变换
在三种常见位运算符（OR，AND，XOR）中，FMT 仅能解决前两种，而对 XOR 卷积却根本无从下手。
于是，这就是快速沃尔什变换（FWT）出场的时候了。
FWT 同 FMT 一样，也可实现 OR 和 AND 卷积。不同的是，FMT 是从DP方向推出的结论，而 FWT 是从变换方向推出的结论。殊途同归，它们在 OR 和 AND 卷积中得到了一样的结果，不同的是 FWT 有着更像 FFT 的结构。但正因如此，OR 和 AND 的 FWT 要比 FMT 难写很多，且不方便理解。所以，这两种情形推荐直接使用 FMT，而把 FWT 留到 XOR 卷积上使用。
具体而言，我们希望得到 \(h=f\operatorname{xor}g\)。构造变换 \(f\rightarrow\hat f\)，我们希望有 \(\hat h(x)=\hat f(x)\times\hat g(x)\)。
现在，我们给出一个定义：设 \(i\odot j=\operatorname{popcount}(i\operatorname{and}j)\bmod2\)，其中 \(\text{popcount}\) 表示计算一个二进制数中 \(1\) 的数量 。
则应有 \((i\odot j)\operatorname{xor}(i\odot k)=i\odot(j\operatorname{xor}k)\)。具体可以按位考虑然后发现每一位在模 \(2\) 意义下全部相等得证。
有了这个性质，我们便设 \(\hat f(x)=\sum\limits_{x\odot y=0}f(y)-\sum\limits_{x\odot y=1}f(y)\)。
于是考虑
于是该 \(f\rightarrow\hat f\) 是合法的 XOR 变换。考虑如何代码实现此变换及其逆变换。
考虑分治实现。在分治的第 \(i\) 层，我们只考虑所有数的最低 \(i\) 位。考虑两个前 \(i-1\) 位全部相同，唯独第 \(i\) 位一个是 \(0\) 一个是 \(1\) 的数。不妨设 \(0\) 的那个数是 \(x\)，\(1\) 的那个数是 \(y\)。则，对于 \(x\)，其与所有数的异或值不变，有 \(f'(x)=f(x)+f(y)\)；对于 \(y\)，因为多了一个 \(1\)，异或值翻转，故有 \(f'(y)=f(x)-f(y)\)。
然后，其逆运算则有 \(f'(x)=\dfrac{f(x)+f(y)}{2},f'(y)=\dfrac{f(x)-f(y)}{2}\)。
同 FFT 类似地实现即可。
时间复杂度 \(O(n2^n)\)，或者说 \(n\log n\)，取决于从哪个角度来看。
I.【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
就是模板啦。
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int inv2=499122177;
int n,a[1<<17],b[1<<17],A[1<<17],B[1<<17],c[1<<17],lim;
void FMTOR(int *f,int tp){
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(f[j]+=(mod+tp*f[j^(1<<i)])%mod)%=mod;
}
void FMTAND(int *f,int tp){
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(f[j^(1<<i)]+=(mod+tp*f[j])%mod)%=mod;
}
void FWTXOR(int *f,int tp){
for(int md=1;md<lim;md<<=1)for(int stp=md<<1,pos=0;pos<lim;pos+=stp)for(int i=0;i<md;i++){
int x=f[pos+i],y=f[pos+md+i];
f[pos+i]=(x+y)%mod;
f[pos+md+i]=(x-y+mod)%mod;
if(tp==-1)f[pos+i]=1ll*f[pos+i]*inv2%mod,f[pos+md+i]=1ll*f[pos+md+i]*inv2%mod;
}
}
void OR(int *f,int *g,int *h){
for(int i=0;i<lim;i++)A[i]=f[i],B[i]=g[i];
FMTOR(A,1),FMTOR(B,1);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
FMTOR(A,-1);
for(int i=0;i<lim;i++)c[i]=A[i];
}
void AND(int *f,int *g,int *h){
for(int i=0;i<lim;i++)A[i]=f[i],B[i]=g[i];
FMTAND(A,1),FMTAND(B,1);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
FMTAND(A,-1);
for(int i=0;i<lim;i++)c[i]=A[i];
}
void XOR(int *f,int *g,int *h){
for(int i=0;i<lim;i++)A[i]=f[i],B[i]=g[i];
FWTXOR(A,1),FWTXOR(B,1);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
FWTXOR(A,-1);
for(int i=0;i<lim;i++)c[i]=A[i];
}
int main(){
scanf("%d",&n),lim=1<<n;
for(int i=0;i<lim;i++)scanf("%d",&a[i]);
for(int i=0;i<lim;i++)scanf("%d",&b[i]);
OR(a,b,c);for(int i=0;i<lim;i++)printf("%d ",c[i]);puts("");
AND(a,b,c);for(int i=0;i<lim;i++)printf("%d ",c[i]);puts("");
XOR(a,b,c);for(int i=0;i<lim;i++)printf("%d ",c[i]);puts("");
return 0;
}
II.CF1119H Triple
开始涉及到 FWT 的本质了。
首先，先思考最暴力的做法：令 \(f_i\) 表示第 \(i\) 对三元组所对应的函数，其中 \(f_i(a_i)=X,f_i(b_i)=Y,f_i(c_i)=Z\)，其余位置均为 \(0\)（这里 \(X,Y,Z\) 大写是为了避免接下来讲解中引起歧义）。重复计算 \(m\) 个 \(f_i\) 的 \(\operatorname{xor}\) 卷积，最终得到的就是答案。复杂度为 \(mn2^n\)（这里认为共有 \(m\) 对三元组，值域是 \(2^n\)），显然吃不消。
考虑 \(f\rightarrow\hat f\) 的变换的实质。若我们把 \(f,\hat f\) 各自看作一个 \(n\) 维向量，则该变换相当于拿 \(f\) 左乘了一个变换矩阵得到 \(\hat f\)。不妨设该矩阵是 \(A\)。则，依据我们对 FWT 的理解，有 \(A_{i,j}=(-1)^{i\odot j}\)。
当我们计算所有 \(f_i\) 的 \(\text{xor}\) 卷积时，可以一次性对所有 \(f\) 作 FWT，将所有 \(\hat f\) 按位作积得到答案函数 \(\hat S\)，然后对 \(\hat S\) 快速沃尔什反演得到真正的 \(S\)。现在，因为 \(f\) 矩阵中有大量的空白，所以在计算 \(Af\) 的时候也不需要全部的东西，仅需枚举 \(a_i,b_i,c_i\) 这三行即可由 \(f\) 得到 \(\hat f\)，复杂度 \(O(2^n)\)。因为我们要对 \(m\) 对三元组各作一次，最终还要作一次常规沃反得到 \(S\)，所以总复杂度 \(O\Big((m+n)2^n\Big)\)，仍然不行。
没办法，我们考虑列出一些式子，从式子出手。
首先，我们的目标是快速求出 \(\hat S\)，因为从 \(\hat S\) 本身即可快速求 \(S\)。我们有 \(\hat S(x)=\prod\limits_{i=1}^m\hat f_i(x)\)。
考虑 \(\hat f(x)\) 的实质，就是 \(XA_{a,x}+YA_{b,x}+ZA_{c,x}\)，也即 \((-1)^{a\odot x}X+(-1)^{b\odot x}Y+(-1)^{c\odot x}Z\)。因为 \(X,Y,Z\) 是固定的，而它们前面的系数一共只有 \(8\) 种可能，所以只需要快速求出 \(8\) 种可能各有多少组，然后快速幂一下即可。
但我们有一种方法可以把 \(8\) 种消作 \(4\) 种——对每一对三元组 \(\{a,b,c\}\)，统一异或上 \(a\) 得到 \(\{0,a\operatorname{xor} b,a\operatorname{xor} c\}\)。则最终只需要在答案的下标上再异或掉一个 \(\operatorname{xor}\limits_{i=1}^ma_i\) 即可。设 \(u=a\operatorname{xor}b,v=a\operatorname{xor}c\)，则上式可以被表示成 \(X+(-1)^{u\odot x}Y+(-1)^{v\odot x}Z\)。
设仅针对 \(\hat S(x)\) 来说，四种情况（\(X+Y+Z,X+Y-Z,X-Y+Z,X-Y-Z\)）的次数分别是 \(t_0,t_1,t_2,t_3\)。考虑解方程。
首先，第一条方程是众所周知的——
下面，考虑找出第二条方程。
因为 \(Y\) 的系数和 \(Z\) 的系数是独立的（一个仅与 \(u\) 有关，一个仅与 \(v\) 有关），故先来考虑 \(Y\) 的系数——\((-1)^{u\odot x}\)。根据我们的观察，这可以被看作是仅有 \(g(u)=1\) 的 \(g\) 作 FWT 的矩阵的转移系数。故我们此处就设出上述的 \(g\)，然后计算 \(\sum\limits_{i=1}^m\hat g(x)\)，就能得到 \(\sum\limits_{i=1}^m(-1)^{u_i\odot x}\)。
因为在 \(t_0,t_1\) 中，\(Y\) 的系数为正的上述东西，而其余两个中则为负的上述东西，所以应有上述值等于 \(t_0+t_1-t_2-t_3\)。很明显这可以作为我们的第二个方程。
于是现在的问题转换为如何求出 \(\sum\limits_{i=1}^m\hat g(x)\)。因为我们前面已经说过 FWT 的实质是向量左乘矩阵，所以其满足分配律，也即只需求出 \((\widehat{\Sigma g})(x)\) 即可。显然，\(\Sigma g\) 是极好求出的。于是我们可以令一个 \(k_1\) 表示该值。
同理，我们也可以仅考虑 \(Z\) 的系数，然后类似地得到一个 \(t_0-t_1+t_2-t_3=k_2\)。
但是现在我们仅有三个方程，还差一个。稍微思考一下，认为这个方程应该由 \(Y,Z\) 共同给出。
比如说，令仅有 \(u\operatorname{xor}v\) 的位置为 \(1\) 的函数，对其 FWT，发现其系数就等于 \((-1)^{u\odot x}(-1)^{v\odot x}\)（这里仍然使用了我们 \((i\odot j)\operatorname{xor}(i\odot k)=i\odot(j\operatorname{xor}k)\) 的式子）。
于是我们得到了第四个方程——\(t_0-t_1-t_2+t_3=k_3\)。
联立得到
手工爆解得到
搞定。时间复杂度 \(O(m+n2^n)\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int inv2=499122177;
const int inv4=748683265;
int m,n,X,Y,Z,u[100100],v[100100],A,lim,k1[1<<17],k2[1<<17],k3[1<<17],S[1<<17];
int ksm(int x,int y){
int z=1;
for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;
return z;
}
void FWT(int *f,int tp){
for(int md=1;md<lim;md<<=1)for(int stp=md<<1,pos=0;pos<lim;pos+=stp)for(int i=0;i<md;i++){
int x=f[pos+i],y=f[pos+md+i];
f[pos+i]=(x+y)%mod;
f[pos+md+i]=(x+mod-y)%mod;
if(tp==-1)f[pos+i]=1ll*f[pos+i]*inv2%mod,f[pos+md+i]=1ll*f[pos+md+i]*inv2%mod;
}
}
int main(){
scanf("%d%d%d%d%d",&m,&n,&X,&Y,&Z),lim=1<<n;
for(int i=1,a,b,c;i<=m;i++)scanf("%d%d%d",&a,&b,&c),u[i]=a^b,v[i]=a^c,A^=a;
for(int i=1;i<=m;i++)k1[u[i]]++;FWT(k1,1);
for(int i=1;i<=m;i++)k2[v[i]]++;FWT(k2,1);
for(int i=1;i<=m;i++)k3[u[i]^v[i]]++;FWT(k3,1);
for(int i=0;i<lim;i++){
int t0=(0ll+m+k1[i]+k2[i]+k3[i])*inv4%mod;
int t1=(1ll*(m+k1[i])*inv2%mod+mod-t0)%mod;
int t2=(1ll*(m+k2[i])*inv2%mod+mod-t0)%mod;
int t3=(1ll*(m+k3[i])*inv2%mod+mod-t0)%mod;
S[i]=1ll*ksm((0ll+X+Y+Z)%mod,t0)*ksm((0ll+X+Y-Z+mod)%mod,t1)%mod*ksm((0ll+X-Y+Z+mod)%mod,t2)%mod*ksm((0ll+X-Y-Z+mod+mod)%mod,t3)%mod;
}
FWT(S,-1);
for(int i=0;i<lim;i++)printf("%d ",S[i^A]);
return 0;
}
子集卷积
定义两个函数的子集卷积是 \(h(i)=\sum\limits_{j\operatorname{or}k=i,j\operatorname{and}k=0}f(j)\times g(k)\)。
首先，\(j\operatorname{or}k\) 的限制直接 FMT 即可。但是 \(j\operatorname{and}k\) 的限制又要怎么办呢？
设 \(|i|\) 表示 \(i\) 二进制表达中 \(1\) 的数量。则，\(j\operatorname{and}k=0\)，等价于 \(|j|+|k|=|j\operatorname{or}k|\)。
于是我们设 \(f_p(x)\) 表示仅保留所有 \(|x|=p\) 的 \(x\) 后所得到的函数。则，应有 \(h_p=\sum\limits_{q=0}^pf_q\operatorname{or}g_{p-q}\)。
显然，对于每个 \(f_p\) 各作一遍 FMT，然后再 \(O(n^2)\) 地暴力进行卷积拼在一块，最后再对得到的所有 \(\hat h_p\) 作 FMI 变回去即可。
时间复杂度 \(O(n^22^n)\)，远优于暴力的 \(O(3^n)\)。
I.【模板】子集卷积
模板即可。
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+9;
#define bit __builtin_popcount
int n,a[21][1<<20],b[21][1<<20],c[21][1<<20],lim;
void FMT(int *f,int tp){for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(f[j]+=(mod+tp*f[j^(1<<i)])%mod)%=mod;}
int main(){
scanf("%d",&n),lim=1<<n;
for(int i=0;i<lim;i++)scanf("%d",&a[bit(i)][i]);
for(int i=0;i<lim;i++)scanf("%d",&b[bit(i)][i]);
for(int i=0;i<=n;i++)FMT(a[i],1),FMT(b[i],1);
for(int i=0;i<=n;i++)for(int j=0;j<=i;j++)for(int k=0;k<lim;k++)(c[i][k]+=1ll*a[j][k]*b[i-j][k]%mod)%=mod;
for(int i=0;i<=n;i++)FMT(c[i],-1);
for(int i=0;i<lim;i++)printf("%d ",c[bit(i)][i]);
return 0;
}
II.CF914G Sum the Fibonacci
经 典 多 合 一
\(a,b\) 这一坨用子集卷积解决，\(d,e\) 这一坨用 FWT 解决，然后三个并一块用 FMT-AND 解决。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int inv2=5e8+4;
#define bit __builtin_popcount
int m,a[18][1<<17],b[1<<17],c[1<<17],d[1<<17],e[1<<17],n=17,lim=1<<17,res;
void FMTOR(int *f,int tp){
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(f[j]+=(mod+tp*f[j^(1<<i)])%mod)%=mod;
}
void FMTAND(int *f,int tp){
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(f[j^(1<<i)]+=(mod+tp*f[j])%mod)%=mod;
}
void FWTXOR(int *f,int tp){
for(int md=1;md<lim;md<<=1)for(int stp=md<<1,pos=0;pos<lim;pos+=stp)for(int i=0;i<md;i++){
int x=f[pos+i],y=f[pos+md+i];
f[pos+i]=(x+y)%mod;
f[pos+md+i]=(x-y+mod)%mod;
if(tp==-1)f[pos+i]=1ll*f[pos+i]*inv2%mod,f[pos+md+i]=1ll*f[pos+md+i]*inv2%mod;
}
}
int main(){
scanf("%d",&m);
for(int i=1,x;i<=m;i++)scanf("%d",&x),a[bit(x)][x]++,c[x]++,d[x]++;
for(int i=0;i<=n;i++)FMTOR(a[i],1);
for(int i=n;i>=0;i--)for(int j=0;j<lim;j++){
int x=0;
for(int k=0;k<=i;k++)(x+=1ll*a[k][j]*a[i-k][j]%mod)%=mod;
a[i][j]=x;
}
for(int i=0;i<=n;i++)FMTOR(a[i],-1);
for(int i=0;i<lim;i++)b[i]=a[bit(i)][i];
FWTXOR(d,1);
for(int i=0;i<lim;i++)d[i]=1ll*d[i]*d[i]%mod;
FWTXOR(d,-1);
e[0]=0,e[1]=1;for(int i=2;i<lim;i++)e[i]=(e[i-1]+e[i-2])%mod;
for(int i=0;i<lim;i++)b[i]=1ll*b[i]*e[i]%mod,c[i]=1ll*c[i]*e[i]%mod,d[i]=1ll*d[i]*e[i]%mod;
FMTAND(b,1),FMTAND(c,1),FMTAND(d,1);
for(int i=0;i<lim;i++)e[i]=1ll*b[i]*c[i]%mod*d[i]%mod;
FMTAND(e,-1);
for(int i=0;i<n;i++)(res+=e[1<<i])%=mod;
printf("%d\n",res);
return 0;
}
