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 (\emptyset) = f (\emptyset)

$g(\varnothing)=f(\varnothing)$

g ({x_{1}}) = f ({x_{1}}) + f (\emptyset)

$g\Big(\{x_1\}\Big)=f\Big(\{x_1\}\Big)+f(\varnothing)$

g ({x_{2}}) = f ({x_{2}}) + f (\emptyset)

$g\Big(\{x_2\}\Big)=f\Big(\{x_2\}\Big)+f(\varnothing)$

g ({x_{1}, x_{2}}) = f ({x_{1}, x_{2}}) + f ({x_{1}}) + f ({x_{2}}) + f (\emptyset)

$g\Big(\{x_1,x_2\}\Big)=f\Big(\{x_1,x_2\}\Big)+f\Big(\{x_1\}\Big)+f\Big(\{x_2\}\Big)+f(\varnothing)$

发现其看上去一副可以容斥的样子。于是我们稍加思索，将最后一个式子改成了这样：

g ({x_{1}, x_{2}}) = f ({x_{1}, x_{2}}) + g ({x_{1}}) + g ({x_{2}}) - g (\emptyset)

$g\Big(\{x_1,x_2\}\Big)=f\Big(\{x_1,x_2\}\Big)+g\Big(\{x_1\}\Big)+g\Big(\{x_2\}\Big)-g(\varnothing)$

再考虑更高维的样式，就有 $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)$ 。

于是考虑

\begin{aligned} \hat{f} (x) \times \hat{g} (x) & = (\sum_{x ⊙ y = 0} f (y) - \sum_{x ⊙ y = 1} f (y)) \times (\sum_{x ⊙ y = 0} g (y) - \sum_{x ⊙ y = 1} g (y)) \\ = (\sum_{x ⊙ y = 1} \sum_{x ⊙ z = 1} f (y) g (z) + \sum_{x ⊙ y = 0} \sum_{x ⊙ z = 0} f (y) g (z)) - (\sum_{x ⊙ y = 1} \sum_{x ⊙ z = 0} f (y) g (z) + \sum_{x ⊙ y = 0} \sum_{x ⊙ z = 1} f (y) g (z)) \\ = \sum_{(x ⊙ y) xor (x ⊙ z) = 0} f (y) g (z) - \sum_{(x ⊙ y) xor (x ⊙ z) = 1} f (y) g (z) \\ = \sum_{x ⊙ (y xor z) = 0} f (y) g (z) - \sum_{x ⊙ (y xor z) = 1} f (y) g (z) \\ = \hat{h} (x) \end{aligned}

$\begin{aligned}\hat f(x)\times\hat g(x)&=\Big(\sum\limits_{x\odot y=0}f(y)-\sum\limits_{x\odot y=1}f(y)\Big)\times\Big(\sum\limits_{x\odot y=0}g(y)-\sum\limits_{x\odot y=1}g(y)\Big)\\&=\Big(\sum\limits_{x\odot y=1}\sum\limits_{x\odot z=1}f(y)g(z)+\sum\limits_{x\odot y=0}\sum\limits_{x\odot z=0}f(y)g(z)\Big)-\Big(\sum\limits_{x\odot y=1}\sum\limits_{x\odot z=0}f(y)g(z)+\sum\limits_{x\odot y=0}\sum\limits_{x\odot z=1}f(y)g(z)\Big)\\&=\sum\limits_{(x\odot y)\operatorname{xor}(x\odot z)=0}f(y)g(z)-\sum\limits_{(x\odot y)\operatorname{xor}(x\odot z)=1}f(y)g(z)\\&=\sum\limits_{x\odot (y\operatorname{xor}z)=0}f(y)g(z)-\sum\limits_{x\odot (y\operatorname{xor}z)=1}f(y)g(z)\\&=\hat h(x)\end{aligned}$

于是该 $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$ 。考虑解方程。

首先，第一条方程是众所周知的——

t_{0} + t_{1} + t_{2} + t_{3} = m

$t_0+t_1+t_2+t_3=m$

下面，考虑找出第二条方程。

因为 $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$ 。

联立得到

{\begin{cases} t_{0} + t_{1} + t_{2} + t_{3} = m \\ t_{0} + t_{1} - t_{2} - t_{3} = k_{1} \\ t_{0} - t_{1} + t_{2} - t_{3} = k_{2} \\ t_{0} - t_{1} - t_{2} + t_{3} = k_{3} \end{cases}

$\begin{cases}t_0+t_1+t_2+t_3=m\\t_0+t_1-t_2-t_3=k_1\\t_0-t_1+t_2-t_3=k_2\\t_0-t_1-t_2+t_3=k_3\end{cases}$

手工爆解得到

{\begin{cases} t_{0} = \frac{m + k_{1} + k_{2} + k_{3}}{4} \\ t_{1} = \frac{m + k_{1}}{2} - t_{0} \\ t_{2} = \frac{m + k_{2}}{2} - t_{0} \\ t_{3} = \frac{m + k_{3}}{2} - t_{0} \end{cases}

$\begin{cases}t_0=\dfrac{m+k_1+k_2+k_3}{4}\\t_1=\dfrac{m+k_1}{2}-t_0\\t_2=\dfrac{m+k_2}{2}-t_0\\t_3=\dfrac{m+k_3}{2}-t_0\end{cases}$

搞定。时间复杂度 $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;
}

posted @ 2021-03-31 15:45 Troverld 阅读(2322) 评论(4) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 阿里最新开源QwQ-32B，效果媲美deepseek-r1满血版，部署成本又又又降低了！
· 单线程的Redis速度为什么快？
· SQL Server 2025 AI相关能力初探
· AI编程工具终极对决：字节Trae VS Cursor，谁才是开发者新宠？
· 展开说说关于C#中ORM框架的用法！

Troverld

博客渲染出问题了那是博客的锅，本人什么都不会做的。