SOSDP（子集DP）、高维前缀和、FMT（快速莫比乌斯变换）、FWT（快速沃尔什变换）学习笔记
这里是集合幂级数学习笔记。
本博客在半年前就已经有计划去写了，只不过当时题目较少没能如愿。现在攒了几道题，开始动笔。
O.定义
定义 \mathsf 字体的函数为集合幂级数。例：\(\mathsf F(x),\mathsf G(x)\)。
定义 \mathbb 字体的变量表示集合。例：\(\mathbb{S,T}\)。
定义常规字体的变量为常规变量（废话）或者形式幂级数（多项式）。例：\(x,y,f(x),h(z)\)。
一些符号：
\(\setminus\) 符号为集合删除符号。
\(\bb{\comp_U(V)}\) 符号为补集符号，表示 \(\bb{U\setminus V}\)。在 \(\bb U\) 为全集时也可省略。
I.集合幂级数
定义集合幂级数为形为 \(\mathsf F(x)=\sum\limits_{\mathbb S}f(\mathbb S)x^{\mathbb S}\) 的东西。
什么？你问一个变量是怎么对集合取幂的？
……应该结合形式幂级数来看，这只是形式啊，形式啦形式。
在明确全集 \(\mathbb U\) 的前提下，明显每个 \(\mathbb S\) 可以双射为一个二进制数 \(s\)，其每一位表示 \(\mathbb U\) 中的某个元素是否在 \(\mathbb S\) 中。
故集合幂级数也可以表示为 \(\sum\limits_{0\leq s<u}f(s)x^s\)。
当然，如果全集无限，显然这种表示法就不好处理了，还是得回到原始定义。
同时我们发现，有 \(\mathbb{I\cup/\cap/\otimes J=K}\Leftrightarrow i\operatorname{or/and/xor} j=k\)。其中 \(\otimes\) 被定义为类似于两个集合取 XOR 之类的操作，也即集合“对称差”。
II.快速莫比乌斯变换/集合并卷积/集合交卷积
我们首先来看一道最板子的例题：
\(|\mathbb U|=n\)。对于全部 \(\mathbb{S}\subseteq\mathbb{U}\)，定义了一函数 \(f(\mathbb S)\)。现在，再定义另一函数 \(g(\mathbb S)=\sum\limits_{\mathbb T\subseteq\mathbb S}f(\mathbb T)\)。现建立集合幂级数 \(\mathsf{F}(x),\mathsf{G}(x)\)。明显 \(\mathsf F(x)\) 已知，求 \(\mathsf G(x)\)。
数据范围：\(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常数更劣。
上述容斥的思想启发我们，什么情形下我们还要用容斥？
二维前缀和，它与上述二维时情形十分类似。
于是我们发现，我们要求的东西本质上是高维前缀和。
我们发现，前面的算法中，我们没有省掉任何的枚举。这启发我们通过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)];}
众所周知，前缀和的逆运算是差分。
虽然高维差分的定义比较困难，但是换一种说法就能很轻易地理解——我们能否找到一种方式，使得已知形式幂级数 \(\mathsf G\)，形式幂级数 \(\mathsf F\)？
明显，因为我们已经知道 \(\mathsf G\) 的递推式，故只需将循环全部反向、加号变减号，即可做到高维差分。
但是，如果我们仔细看看的话，就会发现因为一个 \(\mathbb S\) 中某个元素最多出现一次，所以此二循环的枚举顺序无论正反是均可的。于是，从高维前缀和到高维差分，唯一的变化就是加号变减号。
同理也有高维后缀和。其定义是 \(h(\mathbb S)=\sum\limits_{\mathbb S\subseteq\mathbb T}f(\mathbb T)\)。但是，通过将所有东西取补集，其就变成了高维前缀和，也就可以类似地处理了。同理也可以处理高维后缀差分。
有了之前的铺垫，理解快速莫比乌斯变换也就不难了。
注意是快速莫比乌斯变换，而非莫比乌斯反演。二者一定不能混淆。
我们考虑已知 \(\mathsf{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\) 维的状态（即全集 \(\mathbb U\) 中每个元素在或不在集合中）的话，其等价于 \(I\text{ or }J=S\)（注意这里是普通字体，因为这里并非集合，而是集合对应的常规变量）。故这种操作也可以被叫做 \(\mathsf F\operatorname{or}\mathsf G\)。
现在，考虑对 \(\mathsf F\) 进行高维前缀和操作，得到一个新的函数 \(\hat{\mathsf F}(x)=\sum\limits_{\mathbb S}\hat f(\mathbb S)x^\mathbb S\)。考虑令 \(\hat{h}(\mathbb{S})=\hat{f}(\mathbb{S})\times\hat{g}(\mathbb{S})\)（注意这里是点积，就和常规 FFT 之类类似）。假如我们把 \(\hat{f}(\mathbb S)\) 与 \(\hat{g}(\mathbb S)\) 的定义代入的话，会发现得到的 \(\hat h(\mathbb S)\) 就是 \(\hat{\mathsf H}(x)\) 的 \(x^{\mathbb S}\) 项系数。于是对其作高维差分即可得到 \(\mathsf H\) 本身。
我们发现，这种 \(\hat{\mathsf H}\) 操作，实际上跟FFT中的“DFT”操作特别类似——同样是对一个函数作运算得到另一个函数，然后用运算完的操作进行点积，最后再进行反向操作将其变换为原函数。而整套流程下来，实现了卷积的效果，因此这套流程被叫做“or 卷积”。在原函数和新函数间建立双射的操作，被称作“快速莫比乌斯变换Fast Mobius Transform(FMT)”。
总结一下。通过对函数进行 FMT（高维前缀和），可以得到新函数。对新函数进行按位求积，可以得到答案的新函数。对其进行 FMI（高维差分），得到答案。整一套流程相当于 or 卷积。
同理，我们也可以得到 and 卷积的操作——只不过使用的是高维后缀和和高维后缀差分罢了。
or 卷积与 and 卷积的复杂度都是 \(O(n2^n)\)。但是，如果我们重新令新 \(n\) 表示原本的 \(2^n\) 的话，就会发现复杂度是 \(n\log n\)，同FFT复杂度一致。
因为 or 卷积的实质意义是 \(\mathbb{I\cup J=S}\)，而 and 卷积的实质意义是 \(\mathbb{I\cap J=S}\)，所以二者也分别被叫做“集合并卷积”和“集合交卷积”。
到这里我们先打住吧，来看几道例题。
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;
}
VI.开黑团体
题意：给定 \(m\) 个 \(n\) 位二进制数，求有多少个 \(n\) 位二进制数可以表示成其中若干个数的或。
数据范围：\(n\leq28,m\leq\min(2^n,5\times10^5)\)。
空限 2G，时限 4s。
如何判定一个二进制数能否表示成若干个数的或？
一种方式是把每个数 \(a\) 看成集合幂级数 \(1+x^a\) 然后用 FMT 卷积卷在一块。但是这样不仅能够判定，还能够计数，太过泛用的结果就是无法优化。
令集合幂级数 \(\mathsf F\)，其中对于每个给出的数 \(a\)，均有 \(\mathsf F(a)=a\)。对其作高维前缀或，得到 \(\sf G\)。若有 \(\mathsf G(b)=b\)，则明显 \(b\) 可以被得出（因为其所有子集中数或起来可以得到其本身）。故我们若可以得到 \(\sf G\)，则可以在 \(O(2^n)\) 的时间内实现判定。
现在问题在于 \(\sf G\) 的求出。一种想法是按位处理，对于每位分开做 FMT，然后用 bitset 压位，复杂度 \(\dfrac{2^nn^2}\omega\)。
事实上无法通过。考虑用某种方法把所有位一并求出。
我们首先分析按位处理的情形：当处理 \(i\) 时，考虑全体第 \(i\) 位为 \(1\) 的数 \(a\)，把 \(\mathsf F(a)\) 赋成 \(1\)，然后做高位前缀或。
然后分析高位前缀或的情形：枚举 \(j\)，然后枚举 \(j\) 位上为 \(1\) 的数 \(k\)，若 \(\mathsf F(k)=0\) 且 \(\mathsf F(k-2^j)=1\)，把 \(\mathsf F(k)\) 赋成 \(1\)。
然后发现你若要求出第 \(i\) 位为 \(1\) 的情形，只需在枚举 \(j\) 的时候忽略 \(i\) 即可。用分治处理这个情形，复杂度 \(O(\dfrac{2^nn\log n}\omega)\)，事实上与 \(m\) 无关。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
int n,m,lim;
char buf[1<<23],*p1=buf,*p2=buf;
#ifndef Troverld
#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*(p1++))
#else
#define gc() getchar()
#endif
template<typename T>void read(T&x){
x=0;
char c=gc();
while(c>'9'||c<'0')c=gc();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=gc();
}
ull a[20][1<<22],b[1<<22],msk[6];
int tp,res;
void pr(ull x,int y=0){for(int i=0;i<64;i++)if((x>>i)&1)printf("%d ",(y<<6)|i);puts("");}
void solve(int l,int r){
// pr(a[tp][0]);
if(l==r){
if(l>=6){for(int i=0;i<(1<<22);i++)if(i&(1<<l-6))b[i]&=a[tp][i];}
else for(int i=0;i<(1<<22);i++)b[i]&=a[tp][i]|msk[l];
return;
}
int mid=(l+r)>>1;++tp;
for(int i=0;i<(1<<22);i++)a[tp][i]=a[tp-1][i];
for(int j=mid+1;j<=r;j++){
if(j>=6){for(int i=0;i<(1<<22);i++)if(i&(1<<j-6))a[tp][i]|=a[tp][i^(1<<j-6)];}
else for(int i=0;i<(1<<22);i++)a[tp][i]|=(a[tp][i]&msk[j])<<(1<<j);
}
solve(l,mid);
for(int i=0;i<(1<<22);i++)a[tp][i]=a[tp-1][i];
for(int j=l;j<=mid;j++){
if(j>=6){for(int i=0;i<(1<<22);i++)if(i&(1<<j-6))a[tp][i]|=a[tp][i^(1<<j-6)];}
else for(int i=0;i<(1<<22);i++)a[tp][i]|=(a[tp][i]&msk[j])<<(1<<j);
}
solve(mid+1,r);
tp--;
}
int main(){
freopen("team.in","r",stdin);
freopen("team.out","w",stdout);
for(int i=0;i<6;i++)for(int j=0;j<64;j++)if(!(j&1<<i))msk[i]|=1ull<<j;
// for(int i=0;i<6;i++)pr(msk[i]);
read(n),read(m);
memset(b,0xff,sizeof(b));
for(int i=1,x;i<=m;i++)read(x),a[0][x>>6]|=1ull<<(x&63);
solve(0,27);
for(int i=0;i<(1<<22);i++)res+=__builtin_popcountll(b[i]);
printf("%d\n",res);
return 0;
}
VII.[2022互测#3]整数
有一个长度为 \(n\) 的整数数组 \(a\)，且 \(a_i\in[0,R_i]\)。对于第 \(x\) 个
bit，令 \(\omega(x)=\sum\limits_{i=0}^{n-1}(\left\lfloor\dfrac{a_i}{2^x}\right\rfloor\bmod 2)2^i\)。求满足全体 \(\omega(x)\in\bb S\) 的方案数，其中 \(\bb S\) 是给定的集合。数据范围：\(n\leq18,R_i<2^{60}\)。保证 \(0\in\bb S\)。
显然的想法是从高位往低位数位 DP。记 \(f(\bb T)\) 为当前压上界的集合是 \(\bb T\) 的方案数。令 \(\bb R\) 为这位上为 \(1\) 的位置集合。
一种想法是确定这位上哪些位置贴上界，哪些位置不贴上界。进而把问题转化为一堆集合的交并补，带上一个“\(\bb S\) 中满足某些位置是 \(1\)，某些位置是 \(0\) 的元素计数”的系数。
然后这玩意统计起来很麻烦，进而无法简单搞。
这个做法的本质是确定了贴上界且为 \(1\) 的位中有多少个填 \(1\)，需要的信息是 \(O(2^{|\bb R|})\) 的。同时，我们还要统计不贴上界且不为 \(1\) 的位中有哪些属于 \(\bb S\)，需要的信息是 \(2^{n-|\bb R|}\) 的。总共需要 \(n\) 个 bit 才能描述一种决策。
事实上，我们可以反过来，思考哪些用 \(n\) 个 bit 描述的方案可以得到决策。事实上，我们只需用 \(n\) 个 bit 描述一个 \(\bb S\) 中的元素即可。
考虑我们如果上界状态是 \(\bb T\) 且本位上填的元素是 \(\bb{O\in S}\) 时的场合。我们需要满足：
- \(\bb{(O\cap T)\subseteq R}\)。
- 转移到 \(\bb{T\setminus(R\setminus O)}=\bb{(T\setminus R)\cup(O\cap T)}\)。
我们实际上仅在意 \(\bb{O\cap T}\)。于是，我们 \(2^{n-|\bb R|}\) 枚举 \(\bb{T\setminus R}\)，设为 \(\bb U\)。若一个 \(\bb O\) 能对 \(\bb U\) 产生贡献，必要条件是 \(\bb{O\cap U}=\varnothing\)。于是我们枚举 \(\bb O\)，\(\bb S\) 能贡献给 \(\complement(\bb{O\setminus T})\)。这个直接 FMT 搞一下即可。因为要对每个 \(\bb U\) 记录其有多少个合法的 \(\bb{O\cap R}\)，所以要做 \(2^{|\bb R|}\) 遍 \(2^{n-|\bb R|}\) 的 FMT。
这样，对于 \(\bb U\)，我们可以把所有合法 \(\bb O\) 摊到 \(\bb{O\cap R}\) 上。然后枚举 \(\bb{O\cap R}=\bb V\)，枚举 \(\bb{T\cap R=W}\)，则二者直接转移到 \(\bb{V\cap W}\)。
那么直接搞一个 FMT 卷积即可。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
const int mod=998244353;
void ADD(int&x,int y){if((x+=y)>=mod)x-=mod;}
int n,a[64];
char s[1<<18|1];
ll f[1<<18],g[1<<18];
ll h[1<<18],res;
int main(){
scanf("%d",&n);
for(int i=0;i<n;i++){
ull x;scanf("%lld",&x);
for(int j=0;j<60;j++)a[j]|=((x>>j)&1)<<i;
}
scanf("%s",s);
f[(1<<n)-1]=1;
for(int _=59;_>=0;_--){
// printf("%d(%d):\n",_,a[_]);
// for(int i=0;i<(1<<n);i++)printf("%d ",f[i]);puts("");
int U=(~a[_])&((1<<n)-1);
for(int i=0;i<(1<<n);i++)if(s[i]=='1')h[(i&a[_])|((~i)&U)]++;
for(int i=0;i<n;i++)if(U&(1<<i))
for(int j=0;j<(1<<n);j++)
if(j&(1<<i))
h[j^(1<<i)]+=h[j];
// for(int i=0;i<(1<<n);i++)printf("%d ",h[i]);puts("");
for(int j=0;j<n;j++)if(a[_]&(1<<j))
for(int k=0;k<(1<<n);k++)
if(k&(1<<j))
f[k^(1<<j)]+=f[k],
h[k^(1<<j)]+=h[k];
for(int k=0;k<(1<<n);k++)g[k]=(f[k]%mod)*(h[k]%mod)%mod;
for(int j=0;j<n;j++)if(a[_]&(1<<j))
for(int k=0;k<(1<<n);k++)
if(k&(1<<j))g[k^(1<<j)]-=g[k];
for(int i=0;i<(1<<n);i++){
f[i]=g[i]%mod,g[i]=h[i]=0;
if(f[i]<0)f[i]+=mod;
}
}
for(int i=0;i<(1<<n);i++)res+=f[i];
printf("%d\n",res%mod);
return 0;
}
III.快速沃尔什变换
在三种常见位运算符（OR，AND，XOR）中，FMT 仅能解决前两种，而对 XOR 卷积却根本无从下手。
于是，这就是快速沃尔什变换（FWT）出场的时候了。
FWT 同 FMT 一样，也可实现 OR 和 AND 卷积。不同的是，FMT 是从DP方向推出的结论，而 FWT 是从变换方向推出的结论。殊途同归，它们在 OR 和 AND 卷积中得到了一样的结果，不同的是 FWT 有着更像 FFT 的结构。但正因如此，OR 和 AND 的 FWT 要比 FMT 难写很多，且不方便理解。所以，这两种情形推荐直接使用 FMT，而把 FWT 留到 XOR 卷积上使用。
具体而言，我们希望得到 \(\mathsf{H=F\operatorname{xor}G}\)。构造变换 \(\mathsf{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)\)。
于是考虑
于是该 \(\mathsf{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}\)。（为什么是这个呢？等一会看到 \(K\) 进制 FWT 的时候你就知道啦）
同 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;
}
III.[BZOJ#3696]化合物
在 LCA 处计算贡献。然后发现就是对每个深度的点数构成的序列作异或卷积。那就 FWT 优化做到 \(O(nh\log h)\)。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int lim=512;
void FWT(ll*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++){
ll x=f[pos+i],y=f[pos+md+i];
f[pos+i]=x+y,f[pos+md+i]=x-y;
}
if(tp==-1)for(int i=0;i<lim;i++)f[i]/=lim;
}
int n,dep[100100],cnt[100100][lim];
ll res[lim],A[lim],B[lim];
vector<int>v[100100];
void dfs(int x){
cnt[x][dep[x]]=1;
for(auto y:v[x]){
dep[y]=dep[x]+1,dfs(y);
memset(A,0,sizeof(A)),memset(B,0,sizeof(B));
for(int i=dep[x];i<lim;i++)A[i-dep[x]]=cnt[x][i],B[i-dep[x]]=cnt[y][i],cnt[x][i]+=cnt[y][i];
// for(int i=0;i<4;i++)printf("%d ",A[i]);puts("");
// for(int i=0;i<4;i++)printf("%d ",B[i]);puts("");
FWT(A,1),FWT(B,1);
for(int i=0;i<lim;i++)res[i]+=A[i]*B[i];
}
}
int main(){
scanf("%d",&n);
for(int i=2,x;i<=n;i++)scanf("%d",&x),v[x].push_back(i);
dfs(1);
FWT(res,-1);
int k=lim-1;while(!res[k])k--;
for(int i=0;i<=k;i++)printf("%lld\n",res[i]);
return 0;
}
IV.异或
题意：给定长度为 \(2^n\) 的数组 \(a,b\)。对于每个 \(i\)，求出
\[\sum\limits_{\mathbb S\subseteq\{0,1,\dots,2^n-1\},|\mathbb S|=i}b_{\oplus_{x\in\mathbb S}x}\prod\limits_{x\in\mathbb S}a_x \]对 \(10^9+7\) 取模。
数据范围：\(1\leq n\leq11\)。
把元素写成二元生成函数的形式，其中 \(x\) 维是异或和，\(y\) 维是数量。
则我们要对 \(F_i=1+a_ix^iy\) 在 \(x\) 维作 FWT、\(y\) 维作卷积。
从两方面下手。FWT 完后还要 IFWT，所以我们对 \(b\) 作 IFWT、每个 \(F\) 分开在 \(x\) 上作 FWT，得到
到这个时候就可以扔掉 \(x\) 这一维。以下将上式中 \(y\) 改作 \(x\)。
令 \(G_i(x)=\prod\limits_{j=0}1+(-1)^{i\odot j}a_jx\)。这个式子可以类 FWT 地计算：FWT 中是 \(f_i+f_j\) 与 \(f_i-f_j\)，这里使用 \(G_i(x)G_j(x)\) 与 \(G_i(x)G_y(-x)\) 即可。
卷积用插值代替会让计算简便且迅速（别忘了，模数不是 NTT 模数）。
时间复杂度 \(O(n4^n)\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int inv2=5e8+4;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,lim,a[1<<11],b[1<<11];
void IFWT(){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=(md<<1))for(int i=0;i<md;i++){
int x=b[pos+i],y=b[pos+md+i];
b[pos+i]=1ll*(x+y)*inv2%mod;
b[pos+md+i]=1ll*(x+mod-y)*inv2%mod;
}
}
void FUNC(int*f,int*g){
int F[(1<<11)+1],G[(1<<11)+1];
for(int i=0;i<=lim;i++)F[i]=1ll*f[i]*g[i]%mod,G[i]=1ll*f[i]*g[lim-i]%mod;
for(int i=0;i<=lim;i++)f[i]=F[i],g[i]=G[i];
}
int f[1<<11][(1<<11)+1];
void FWT(){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=(md<<1))
for(int i=0;i<md;i++)FUNC(f[pos+i],f[pos+md+i]);
}
int g[(1<<11)+1];
int coe[(1<<11)+2],eoc[(1<<11)+2],h[(1<<11)+1];
void mul(int k){
for(int i=lim+1;i>=0;i--){
coe[i]=1ll*coe[i]*k%mod;
if(i)(coe[i]+=coe[i-1])%=mod;
}
// for(int i=0;i<=lim+1;i++)printf("%d ",coe[i]);puts("");
}
void div(int k){
if(!k){for(int i=0;i<=lim;i++)eoc[i]=coe[i+1];return;}
int inv=ksm(k);
for(int i=0;i<=lim;i++){
eoc[i]=coe[i];
if(i)(eoc[i]+=mod-eoc[i-1])%=mod;
eoc[i]=1ll*eoc[i]*inv%mod;
}
// for(int i=0;i<=lim;i++)printf("%d ",eoc[i]);puts("");
}
int main(){
freopen("xor.in","r",stdin);
freopen("xor.out","w",stdout);
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]);
IFWT();
for(int i=0;i<lim;i++)for(int j=0;j<=lim;j++)f[i][j]=(1ll*a[i]*(j+mod-(lim>>1))+1)%mod;
FWT();
for(int j=0;j<=lim;j++)for(int i=0;i<lim;i++)g[j]=(1ll*f[i][j]*b[i]+g[j])%mod;
// for(int i=0;i<=lim;i++)printf("%d ",g[i]);puts("");
coe[0]=1;
for(int i=0;i<=lim;i++)mul((mod+(lim>>1)-i)%mod);
for(int i=0;i<=lim;i++){
int lam=g[i];div((mod+(lim>>1)-i)%mod);
for(int j=0;j<=lim;j++)if(i!=j)lam=1ll*lam*ksm((i-j+mod)%mod)%mod;
for(int j=0;j<=lim;j++)h[j]=(1ll*lam*eoc[j]+h[j])%mod;
}
for(int i=1;i<=lim;i++)printf("%d ",h[i]);
return 0;
}
V.游戏
题意：给定一张 \(n\) 个点的图，图中的每条边上放了至少一堆石子，且石子数量 \(a\) 非零。
现在随机求出这张图的一棵生成树，然后选中的边上的每堆石子都有 \(\dfrac 34\) 的概率被清除、\(\dfrac14\) 的概率保留。
求用剩余石子堆玩 NIM 游戏先手必胜的概率。
数据范围：\(1\leq n\leq70,0\leq a<2^m\)，其中 \(m\) 是给定的常数，且 \(1\leq m\leq10\)。每条边上石子堆数不超过 \(1000\)。
NIM 游戏的异或等价性应该不用说了吧。定义概率向量 \(\vec p_i\) 表示第 \(i\) 条边上石子的异或和取到每个值的概率。则我们需要计算的就是
中除零外所有项之和，再除以生成树总数。
把 \(\vec p_i\) FWT 后，异或卷积就变成点积。点积对于每一位来说是独立的，于是我们对于每位分开作矩阵树定理再合并即可得到最终答案向量。
现在问题在于如何求出 \(\hat p_i\)，即 FWT 后的 \(\vec p_i\)。
\(\vec p_i\) 本身等于其中所有石子堆 \(a\) 的概率向量的异或卷积。而 \(a\) 对应的概率向量是 \(\dfrac34+\dfrac14x^a\)。对其 FWT 得到 \(\sum\limits_{i}\Big(\dfrac 12\Big)^{i\odot a}x^i\)，即仅在 \(i\odot a=1\) 的位上有系数 \(\dfrac12\)。
暴力对于每个 \(a\) 在所有系数非 \(1\) 的位上点积是不行的，因为其极限是 \(1000n^22^m\) 级别的。
但是这个式子是可以类 FWT 计算的。我们维护每个位置乘上的 \(\dfrac12\) 的系数，然后就和 FWT 类似了。（事实上，FWT 就可以被理解为一种计算与 \(i\odot j\) 有关的式子的次序）
时间复杂度 \(O(n^2m2^m+n^32^m)\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int inv2=499122177;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,m,lim,invlim;
int RED(int x){return x>=mod?x-mod:x;}
int M[1200][80][80];
int solve(int t){
int tp=1;
for(int i=1,j;i<n;i++){
j=i;
while(j<n&&!M[t][j][i])j++;
if(j>=n)return 0;
if(j!=i){swap(M[t][i],M[t][j]);tp=-tp;}
int inv=RED(mod-ksm(M[t][i][i]));
for(j=i+1;j<n;j++){
int lam=1ll*M[t][j][i]*inv%mod;
for(int k=i;k<n;k++)M[t][j][k]=(1ll*M[t][i][k]*lam+M[t][j][k])%mod;
}
}
tp=RED(tp+mod);
for(int i=1;i<n;i++)tp=1ll*tp*M[t][i][i]%mod;
return tp;
}
void read(int&x){
x=0;
char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
}
int A[1200],B[1200];
void INS(int x){if(__builtin_popcount(x)&1)B[x]++;else A[x]++;}
void FWT(){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=(md<<1))for(int i=0;i<md;i++){
int a1=A[pos+i],b1=B[pos+i],a2=A[pos+md+i],b2=B[pos+md+i];
A[pos+i]=RED(a1+a2);
B[pos+i]=RED(b1+b2);
A[pos+md+i]=RED(a1+b2);
B[pos+md+i]=RED(b1+a2);
}
}
int main(){
freopen("gemo.in","r",stdin);
freopen("gemo.out","w",stdout);
read(n),read(m),invlim=ksm(lim=1<<m);
for(int i=0;i<n;i++)for(int j=i+1,x,y;j<n;j++){
read(x);
if(!x)continue;
while(x--)read(y),INS(y);
FWT();
for(int k=0;k<lim;k++)M[k][i][j]=ksm(inv2,B[k]),A[k]=B[k]=0;
M[lim][i][j]=1;
for(int k=0;k<=lim;k++){
M[k][i][i]=RED(M[k][i][i]+M[k][i][j]);
M[k][j][j]=RED(M[k][j][j]+M[k][i][j]);
M[k][i][j]=M[k][j][i]=RED(mod-M[k][i][j]);
}
}
// for(int i=0;i<n;i++){for(int j=0;j<n;j++)M[i][j].print();puts("");}
int res=0;
for(int i=0;i<lim;i++)res=RED(res+solve(i));
res=1ll*res*invlim%mod;
res=RED(mod-res);
res=1ll*res*ksm(solve(lim))%mod;
printf("%d\n",(res+1)%mod);
return 0;
}
VI.[AGC034F] RNG and XOR
从 \(0\) 开始求出每个值首次出现的时间是没有希望的。事实上，因为是异或，\(0\to x\) 的期望时间等于 \(x\to0\) 的期望时间，于是我们就可以设 \(e_x\) 表示 \(x\to0\) 的期望时间。
事实上，这是期望题中常见的计算某值到达终点的期望时长的套路。
令 \(P\) 表示概率集合幂级数，于是有转移式 \(e_x=\begin{cases}0&(x=0)\\1+\sum\limits_ip_ie_{x\operatorname{xor}i}&(x\neq0)\end{cases}\)。
换句话说，\(E=I+EP+c_0\)，其中 \(c_0\) 是修正 \(e_0\) 的系数。
我们需要解出 \(c_0\)。怎么办呢？方法之一就是对左右同时求一个系数和。记 \(F\) 的系数和为 \(S(F)\)，则上式变成 \(S(E)=S(I)+S(E)S(P)+c_0\)。
关键的一点来了：因为 \(S(P)=1,S(I)=2^n\)，因此左右两半简单消掉了，得到 \(c_0=-2^n\)。
另：\(c_0\) 的值也可以手动对比 \(E\) 和 \(I+EP\) 的差量得出。
然后就可以把上式变成 \(E(1-P)=I-2^n\)。
FWT 求逆仅在非零项生效。注意依照 FWT 的定义，\(\hat P\) 的系数除了 \(\hat P_0=\sum P_i\) 以外，其它项的系数均并非所有项之和。又注意到所有项都是整数，所以除了 \(\hat P_0=1\) 外，其它项小于一。进而，\(1-P\) 在 FWT 后除了常数项外，其它项都非零，故可直接求逆得出。令 \(F\) 表示经历上述操作后，除了 \(F_0\neq\hat E_0\) 外其它项均相等的多项式。
在人脑中对 \(\hat E\) 作一遍 IFWT，有 \(E_0=\dfrac{\sum\hat E_i}{2^n}\)。又有 \(E_0=0\)，故 \(\hat E_0\) 可以直接由上式得出。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,p[1<<18],q[1<<18],lim,invlim,S;
void FWT(int*a,int tp){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=(md<<1))for(int i=0;i<md;i++){
int x=a[pos+i],y=a[pos+md+i];
a[pos+i]=(x+y)%mod;
a[pos+md+i]=(x+mod-y)%mod;
}
if(tp==-1)for(int i=0;i<lim;i++)a[i]=1ll*a[i]*invlim%mod;
}
int main(){
scanf("%d",&n),invlim=ksm(lim=1<<n);
for(int i=0;i<lim;i++)scanf("%d",&p[i]),S+=p[i];
S=ksm(S);
for(int i=0;i<lim;i++)p[i]=(mod-1ll*p[i]*S%mod)%mod,q[i]=1;
(++p[0])%=mod;
(q[0]+=mod-lim)%=mod;
// for(int i=0;i<lim;i++)printf("%d ",p[i]);puts("");
// for(int i=0;i<lim;i++)printf("%d ",q[i]);puts("");
FWT(p,1),FWT(q,1);
// for(int i=0;i<lim;i++)printf("%d ",p[i]);puts("");
// for(int i=0;i<lim;i++)printf("%d ",q[i]);puts("");
S=0;for(int i=1;i<lim;i++)p[i]=1ll*q[i]*ksm(p[i])%mod,(S+=p[i])%=mod;
p[0]=(mod-S)%mod;
FWT(p,-1);
for(int i=0;i<lim;i++)printf("%d\n",p[i]);
return 0;
}
VII.[ZJOI2019]开关
和上题一样，我们首先得到 \(E(1-P)=I-2^n\) 的式子。
转成 \(\widehat{E(1-P)}=\widehat{I-2^n}\)，进而变成 \(\hat E\cdot\widehat{1-P}=\widehat{I-2^n}\)。
记 \(F=1-P,G=I-2^n\)，则 \(\hat E\cdot\hat F=\hat G\)。
\(\hat F=\hat 1-\hat P=I-\hat P\)，\(\hat G=\hat I-\hat{2^n}=\sum\limits_{\mathbb S}x^{\mathbb S}2^{n-|\mathbb S|}\sum\limits_{i=0}^{|\mathbb S|}\dbinom{|\mathbb S|}{i}(-1)^{i}-2^nI\)。
依据二项式定理，上式变成 \(\hat G=\sum\limits_{\mathbb S}x^\mathbb S2^{n-|\mathbb S|}0^{|\mathbb S|}-2^nI=2^n-2^nI\)。
令 \(H=\hat E\)，则 \(\forall\mathbb S\neq\varnothing,H_\mathbb S=\dfrac{-2^n}{1+\sum\limits_{i\in\mathbb S}p_i-\sum\limits_{i\notin\mathbb S}p_i}\)。
考虑 \(H_{\varnothing}\)。有 \(H_{\varnothing}=-\sum\limits_{\mathbb S\neq\varnothing}H_\mathbb S\)。
考虑对 \(H\) 作 IFWT 得到 \(E\)。众所周知，IFWT 就等于 FWT 的结果除以 \(2^n\)，故
到这一步为止是大多数人都能推到的。但是然后呢？然后这个式子就有办法算了吗？
我苦恼了很久，最终看了题解，然后才注意到：\(\sum p_i\leq5\times10^4\)。
这不就可以了吗？直接随便 DP 一下即可。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,p[110],f[50100][2],s[110],res;
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%d",&s[i]);
for(int i=1;i<=n;i++)scanf("%d",&p[i]);
f[0][0]=1;
int S=0;
for(int i=1;i<=n;i++){
for(int j=S;j>=0;j--)for(int k=0;k<2;k++)(f[j+p[i]][k^s[i]]+=f[j][k])%=mod;
S+=p[i];
}
for(int j=1;j<=S;j++)res=(1ll*ksm(j)*f[j][1]+res)%mod;
printf("%d\n",1ll*res*S%mod);
return 0;
}
VIII.线性基
题意：给定若干个 \([0,2^m)\) 范围中的整数，其中 \(i\) 有 \(a_i\) 个。注意到就算两个数的 值相同，它们也被看做是 不同的数。
你可以选择至多 \(b_i\) 个元素 \(i\) 加入一个可重集；定义一个可重集的权值为其中所有元素的异或和，求每个值对应的可重集个数。
数据范围：\(m\leq20,0\leq b_i\leq a_i\leq10^5\)。
首先元素 \(i\) 在选择奇数个时有 \(i\) 的贡献，选择偶数个时则有 \(0\) 的贡献。
这个贡献是杨辉三角行上前缀和的形式，故没有封闭形式，直接分块求出即可。
于是我们得到 \(f_i(x)=A+Bx^i\) 的形式，然后只要求出全体 \(f\) 的异或卷积即可。
首先先提取一下，变成 \(\prod A\times \prod(1+\dfrac BAx^i)\) 的形式。然后对这个东西用魔改后的 FWT 处理即可。
复杂度 \(m2^m\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int BBB=100;
const int N=100000;
const int mod=998244353;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,lim,a[1<<20],b[1<<20],c[1<<20],fac[100100],inv[100100],pro=1;
void read(int&x){
x=0;
char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
}
int C(int x,int y){return x<y?0:1ll*fac[x]*inv[y]%mod*inv[x-y]%mod;}
int f[100100];
vector<int>v[100100];
int calc(int _a,int _b){
if(!_a)return 1;
int ret=C(_a-1,_b);
return _b&1?mod-ret:ret;
}
int d[1<<20],e[1<<20];
void WMUL(){
for(int i=0;i<lim;i++)e[i]=1;
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=md<<1)for(int i=0;i<md;i++){
int dx=d[pos+i],dy=d[pos+md+i];
int ex=e[pos+i],ey=e[pos+md+i];
d[pos+i]=1ll*dx*dy%mod;
e[pos+i]=1ll*ex*ey%mod;
d[pos+md+i]=1ll*dx*ey%mod;
e[pos+md+i]=1ll*ex*dy%mod;
}
// for(int i=0;i<lim;i++)printf("%d ",e[i]);
}
void IFWT(){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=md<<1)for(int i=0;i<md;i++){
int x=e[pos+i],y=e[pos+md+i];
e[pos+i]=(x+y)%mod;
e[pos+md+i]=(x+mod-y)%mod;
}
int invlim=ksm(lim);
for(int i=0;i<lim;i++)e[i]=1ll*e[i]*invlim%mod;
}
void print(int x){if(x>=10)print(x/10);putchar('0'+x%10);}
int main(){
freopen("base.in","r",stdin);
freopen("base.out","w",stdout);
fac[0]=1;for(int i=1;i<=N;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[N]=ksm(fac[N]);for(int i=N;i;i--)inv[i-1]=1ll*inv[i]*i%mod;
read(n),lim=1<<n;
for(int i=0;i<lim;i++)read(a[i]),v[a[i]].push_back(i);
for(int i=0;i<lim;i++)read(b[i]);
for(int i=0;i<=N;i+=BBB){
f[0]=1;for(int j=1;j<=i;j++)f[j]=(f[j-1]+C(i,j))%mod;
for(int j=i;j<i+BBB;j++)for(auto k:v[j]){
c[k]=f[min(b[k],i)];
for(int t=i;t<j;t++)c[k]=(2ll*c[k]+mod-C(t,b[k]))%mod;
}
}
for(int i=0;i<lim;i++)pro=1ll*pro*c[i]%mod,d[i]=1ll*calc(a[i],b[i])*ksm(c[i])%mod;
// for(int i=0;i<lim;i++)printf("[%d,%d]",c[i],d[i]);
WMUL();
IFWT();
for(int i=0;i<lim;i++)print(1ll*e[i]*pro%mod),putchar(' ');
// for(int i=0;i<lim;i++)printf("%d ",c[i]);puts("");
return 0;
}
IX.[UOJ#310][UNR #2]黎明前的巧克力
两个集合异或和相同？那么它们的并集的异或和就会恰为零。
于是我们找到所有异或和为零的子集 \(S\)，则其贡献即为 \(2^{|S|}\)。
于是我们可以把一个大小为 \(a\) 的物品的贡献写成集合幂级数 \(1+2x^a\)，则我们要计算所有 \(1+2x^a\) 的异或卷积的常数项系数。
可以使用魔改 FWT，但是……
我们有更好的做法。
考虑 \(1+2x^a\) 的性质：其在一些位上为 \(3\)，一些位上为 \(-1\)。
那么我们如果计算 \(\sum 1+2x^a\) 的 FWT，由 FWT 的线性性，每位上的结果恰为所有 \(1+2x^a\) FWT 后该位上结果的和。
知道了和，又知道这为上总共有 \(n\) 个函数贡献给它，那么就可以通过 解方程 解出这位上有多少个 \(-1\) 和多少个 \(3\)，然后简单快速幂即可。
复杂度对数。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
const int LG=20,lim=1<<LG;
int n,a[lim];
void FWT(int tp){
for(int md=1;md<lim;md<<=1)for(int pos=0;pos<lim;pos+=md<<1)for(int i=0;i<md;i++){
int x=a[pos+i],y=a[pos+md+i];
a[pos+i]=(x+y)%mod,a[pos+md+i]=(x+mod-y)%mod;
}
if(tp==-1){int invlim=ksm(lim);for(int i=0;i<lim;i++)a[i]=1ll*a[i]*invlim%mod;}
}
int main(){
scanf("%d",&n);
for(int i=1,x;i<=n;i++)scanf("%d",&x),a[0]++,a[x]+=2;
FWT(1);
for(int i=0;i<lim;i++){
int _1=((n*3+mod-a[i])%mod)>>2,_3=n-_1;
a[i]=ksm(3,_3);if(_1&1)a[i]=(mod-a[i])%mod;
}
FWT(-1);
printf("%d\n",(a[0]+mod-1)%mod);
return 0;
}
X.Virtual Self
这题和上题不一样：它没法简单地魔改 FWT 得到。
具体而言，它限制了恰选择 \(m\) 个多项式；这意味着我们如果要直接求卷积的话，就势必要维护 \(m\) 一维的大小，是不好的。
因此，本题只能运用解方程的方法，计算 \(x^a\) 关于每项的系数，然后处理从中选取恰 \(m\) 个的方案数。这关于每项是独立的，因为我们要计算的其实是 \(\sum\limits_{|S|=m}\prod\limits_{i\in S}x^{a_i}\)，而 \(x^{a_i}\) 在 FWT 后具有线性性，无论是相加还是求和就只针对于单点，所以可以解方程。
\(x^a\) 的系数是 \(\pm1\)。那么我们解出有 \(t\) 个 \(1\) 和 \(n-t\) 个 \(-1\)，则这位上的结果就是 \([x^m](1+x)^t(1-x)^{n-t}\)。
枚举分配给两式的系数，我们有结果为
这个式子是很棒的，因为它就化成了卷积的形式，直接卷积处理即可。复杂度对数。
有个神必防 AK test case，懒得搞了，不管了。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace Poly{
const int N=1<<21;
const int mod=998244353;
const int G=3;
int rev[N],povG[N],invG[N];
int ksm(int x,int y=mod-2){int z=1;for(;y;x=1ll*x*x%mod,y>>=1)if(y&1)z=1ll*z*x%mod;return z;}
void init(){for(int md=1;md<N;md<<=1)povG[md]=ksm(G,(mod-1)/(md<<1)),invG[md]=ksm(povG[md]);}
int lim,invlim;
int RED(int x){return x>=mod?x-mod:x;}
void ADJ(int&x){if(x>=mod)x-=mod;}
void NTT(int *a,int tp){
for(int i=0;i<lim;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int md=1;md<lim;md<<=1)
for(int pos=0,rt=(tp==-1?invG[md]:povG[md]);pos<lim;pos+=md<<1)
for(int i=0,w=1;i<md;i++,w=1ll*w*rt%mod){
int x=a[pos+i],y=1ll*w*a[pos+md+i]%mod;
a[pos+i]=RED(x+y),a[pos+md+i]=RED(x+mod-y);
}
if(tp==-1)for(int i=0;i<lim;i++)a[i]=1ll*a[i]*invlim%mod;
}
int A[N],B[N];
void mul(int*a,int*b,int*c,int LG){//using: Array A and B
invlim=ksm(lim=1<<LG);
for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(LG-1));
for(int i=0;i<lim;i++)A[i]=B[i]=0;
for(int i=0;i<(lim>>1);i++)A[i]=a[i],B[i]=b[i];
NTT(A,1),NTT(B,1);
for(int i=0;i<lim;i++)A[i]=1ll*A[i]*B[i]%mod;
NTT(A,-1);
for(int i=0;i<lim;i++)c[i]=A[i];
}
}
using namespace Poly;
int fac[N],inv[N],LG,n,m,a[N],INV[N];
void FWT(){
for(int md=1;md<(1<<LG);md<<=1)for(int pos=0;pos<(1<<LG);pos+=md<<1)for(int i=0;i<md;i++){
int x=a[pos+i],y=a[pos+md+i];
a[pos+i]=x+y,a[pos+md+i]=x-y;
}
}
void IFWT(){
for(int md=1;md<(1<<LG);md<<=1)for(int pos=0;pos<(1<<LG);pos+=md<<1)for(int i=0;i<md;i++){
int x=a[pos+i],y=a[pos+md+i];
a[pos+i]=(x+y)%mod,a[pos+md+i]=(x+mod-y)%mod;
}
int invlim=ksm(1<<LG);
for(int i=0;i<(1<<LG);i++)a[i]=1ll*a[i]*invlim%mod;
}
int f[N],g[N],bnm[N],res;
int main(){
scanf("%d%d",&m,&LG),init();
for(int i=0;i<(1<<LG);i++)scanf("%d",&a[i]),n+=a[i];
FWT();
fac[0]=1;for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[n]=ksm(fac[n]);for(int i=n;i;i--)inv[i-1]=1ll*inv[i]*i%mod;
INV[1]=1;for(int i=2;i<=n;i++)INV[i]=1ll*INV[mod%i]*(mod-mod/i)%mod;
for(int i=0;i<=m;i++){
f[i]=1ll*fac[m]*inv[i]%mod*inv[m-i]%mod;
if((m-i)&1)f[i]=mod-f[i];
}
for(int i=0,pro=1;i<=n;i++)g[i]=pro,pro=1ll*pro*(n-m-i)%mod*INV[i+1]%mod;
for(int i=0,pro=1;i<=n;i++)bnm[i]=pro,pro=1ll*pro*(n-i)%mod*INV[i+1]%mod;
// for(int i=0;i<=n;i++)printf("%d ",f[i]);puts("");
// for(int i=0;i<=n;i++)printf("%d ",g[i]);puts("");
// for(int i=0;i<=n;i++)printf("%d ",bnm[i]);puts("");
int K=0;while((1<<K)<=n)K++;mul(f,g,f,K+1);
for(int i=0;i<=n;i++)f[i]=1ll*f[i]*bnm[m]%mod*ksm(bnm[i])%mod;
// for(int i=0;i<=n;i++)printf("%d ",f[i]);puts("");
// for(int i=0;i<(1<<LG);i++)printf("%d ",a[i]);puts("");
for(int i=0;i<(1<<LG);i++)a[i]=f[(a[i]+n)>>1];
// for(int i=0;i<(1<<LG);i++)printf("%d ",a[i]);puts("");
IFWT();
// for(int i=0;i<(1<<LG);i++)printf("%d ",a[i]);puts("");
for(int i=0;i<(1<<LG);i++)res^=1ll*(i+1)*a[i]%mod;
printf("%d\n",res);
return 0;
}
IV.FWT 思想的应用
FWT 作为一种较为 simple（？）的数学模型，其思想可以被用来出题，也可以被用来解题。
XI.[UOJ#328][UTR #3]量子破碎
分析一下我们的操作：基于概率的询问和对数组的奇怪更新。
因为每次询问都会重构数组，所以我们往往只有两种解决方案：通过操作把所有值全都集中在目标处（也即 \(x\operatorname{xor}y\) 处），然后通过一次操作来百分百确定之；或者，通过多次操作缩小可能范围，最终得以确定结果。
那么这个奇怪更新让我们联想到 FWT。因为我们要保证 \(\sum a^2=1\)，所以将 FWT 矩阵乘上常数后变成 \(\begin{bmatrix}\dfrac1{\sqrt2}&\dfrac1{\sqrt2}\\\dfrac1{\sqrt2}&-\dfrac1{\sqrt2}\end{bmatrix}\)，然后带进去玩一玩，我们得到的数组是正常 FWT 乘上 \(\dfrac1{2^{n/2}}\) 的结果。
因为我们只有 \(a_x=a_y=\dfrac1{\sqrt2}\)，所以用上述奇怪 FWT 更新后就得到了 \(a_i=\dfrac1{\sqrt2\times2^{n/2}}\Big((-1)^{|x\operatorname{and}i|}+(-1)^{|y\operatorname{and}i|}\Big)\)。
这个式的值仅可能有三种，即 \(\pm C\) 和 \(0\)（其中 \(C\) 是某常数，太长省略不写），而 \(\pm C\) 平方后的结果相同，因此我们会等概率返回满足 \((-1)^{|x\operatorname{and}i|}=(-1)^{|y\operatorname{and}i|}\) 的 \(i\) 之一。
因为 \((|x\operatorname{and}i|\operatorname{and}1)\operatorname{xor}(|y\operatorname{and}i|\operatorname{and}1)=(|(x\operatorname{xor}y)\operatorname{and}i|\operatorname{and}1)=0\)，所以这么搞一遍，我们期望能 ban 掉一半的备选方案。
每轮操作期望干掉一半的数，期望 \(n\) 轮出结果。一轮操作需要 \(n+1\) 次，总次数期望 \(n(n+1)\)，实际使用的时候可以跑满 \(23\) 次，这样挂掉的概率就极低了。
代码：
#include"quantumbreak.h"
#include<bits/stdc++.h>
using namespace std;
const double C=1.0/sqrt(2);
bool ban[1<<16];
int query_xor(int n,int){
double a[2][2]={{C,C},{C,-C}};
for(int _=0;_<400/(n+1);_++){
for(int i=0;i<n;i++)manipulate(a,i);
int x=query();
// printf("HAHAHAI:%d\n",x);
for(int i=0;i<(1<<n);i++)if(__builtin_popcount(i&x)&1)ban[i]=true;
}
for(int i=1;i<(1<<n);i++)if(!ban[i])return i;
}
V.\(K\) 进制FWT
这部分内容因为偏多项式而非位运算，所以经过本人的考虑已经转到了FFT理性瞎扯中进行讲解（连带着IFWT中除以二的原因）。
VI.子集卷积
定义两个集合幂级数的子集卷积是 \(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|\)。
于是我们设 \(\mathsf F_p(x)\) 表示仅保留所有 \(|x|=p\) 的 \(x\) 后所得到的集合幂级数。则，应有 \(\mathsf H_p=\sum\limits_{q=0}^p\mathsf F_q\operatorname{or}\mathsf G_{p-q}\)。
显然，对于每个 \(\mathsf F_p\) 各作一遍 FMT，然后再 \(O(n^2)\) 地暴力进行卷积拼在一块，最后再对得到的所有 \(\hat{\mathsf 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;
}
III.[WC2018]州区划分
首先所有合法的集合（要么不存在欧拉回路，要么不连通）可以搜出。
然后DP就行了。设 \(S(\mathbb S)\) 表示 \(\mathbb S\) 集合中所有元素之和，则有 \(g'(\mathbb S)=\dfrac 1{S(\mathbb S)}\sum\limits_{\mathbb{I\cup J=S,I\cap J=\varnothing}}g(\mathbb I)S(\mathbb J)\)，其中 \(g'\) 是下一轮的DP数组，\(g\) 则是上一轮。
发现这么转移一定是从 \(|\mathbb S|\) 小的状态向 \(|\mathbb S|\) 大的状态转移。于是就用卷积依次转移即可。
时间复杂度 \(O(n^22^n)\)。实现可以很不精细。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
#define bit(x) __builtin_popcount(x)
int n,lim;
void FMT(int*a,int tp){for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(a[j]+=(mod+tp*a[j^(1<<i)])%mod)%=mod;}
int m,p,f[22][1<<21],g[22][1<<21],X[410],Y[410],dsu[21],cnt[21],a[21],s[1<<21],t[1<<21];
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
void merge(int x,int y){x=find(x),y=find(y);if(x!=y)dsu[x]=y;}
int main(){
scanf("%d%d%d",&n,&m,&p),lim=1<<n;
for(int i=1;i<=m;i++)scanf("%d%d",&X[i],&Y[i]),X[i]--,Y[i]--;
for(int i=0;i<n;i++)scanf("%d",&a[i]);
for(int i=0;i<lim;i++){
for(int j=0;j<n;j++)if(i&(1<<j))cnt[j]=0,dsu[j]=j,s[i]+=a[j];
s[i]=ksm(s[i],p),t[i]=ksm(s[i]);
f[bit(i)][i]=s[i];
for(int j=1;j<=m;j++)if(((1<<X[j])&i)&&((1<<Y[j])&i))cnt[X[j]]^=1,cnt[Y[j]]^=1,merge(X[j],Y[j]);
bool ok=false;
for(int j=0;j<n&&!ok;j++)if(i&(1<<j))ok|=cnt[j],find(j);
if(ok)continue;
for(int j=0;j<n&&!ok;j++)if(i&(1<<j))for(int k=j+1;k<n;k++)if((i&(1<<k))&&dsu[j]!=dsu[k]){ok=true;break;}
if(ok)continue;
f[bit(i)][i]=0;
}
for(int i=0;i<=n;i++)FMT(f[i],1);
memcpy(g,f,sizeof(f));
for(int i=0;i<=n;i++){
FMT(g[i],-1);
// for(int k=0;k<lim;k++)printf("%d ",g[i][k]);puts("");
for(int k=0;k<lim;k++)g[i][k]=1ll*g[i][k]*t[k]%mod;
FMT(g[i],1);
for(int j=1;i+j<=n;j++)for(int k=0;k<lim;k++)(g[i+j][k]+=1ll*g[i][k]*f[j][k]%mod)%=mod;
}
FMT(g[n],-1);
printf("%d\n",g[n][lim-1]);
return 0;
}
IV.CF1034E Little C Loves 3 III
乍一看是子集卷积板子；但是事实上本题时限太小导致子集卷积压根过不去。
怎么办呢？我们考虑原本的集合幂级数是 \(\sf F\)，其在 \(x^\mathbb S\) 前的系数是 \(f(\mathbb S)\)。
子集卷积的时候把它换成了占位幂级数 \(f(\mathbb S)y^{|\mathbb S|}\)。
但是事实上，这个 \(y\) 只要取一个足够大的数，就可以保证不同的系数间不会干涉。
那么具体应该取多大的 \(y\) 呢？事实上，在无取模的意义下只要取任一比结果大的数即可，而在取模的意义下更简单，只需取模数即可。本题中的模数只有 \(4\)，则我们用 \(f(\mathbb S)4^{|\mathbb S|}\) 的占位幂级数即可。卷积结束后只需除掉 \(4^{|\mathbb S|}\) 即可。
代码：
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<21)+5;
typedef unsigned long long ull;
int n,lim;
char s[N],t[N];
void FMT(ull*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)];
}
ull f[N],g[N];
int main(){
scanf("%d%s%s",&n,s,t),lim=1<<n;
for(int i=0;i<lim;i++)
f[i]=(ull)(s[i]-'0')<<(__builtin_popcount(i)<<1),g[i]=(ull)(t[i]-'0')<<(__builtin_popcount(i)<<1);
FMT(f,1),FMT(g,1);
for(int i=0;i<lim;i++)f[i]*=g[i];
FMT(f,-1);
for(int i=0;i<lim;i++)printf("%d",(f[i]>>(__builtin_popcount(i)<<1))&3);
return 0;
}
V.[TopCoder13347]AquaparkPuzzle
令 \(n=|c|\)，则 \(n\leq11\)。那么我们显然考虑状压……吗？
\(3^{11}\) 能状压，但是转移是不是需要矩阵快速幂呢？
因此我们认为，不能直接上裸的状压，还得稍微分析分析。
于是我们自然想到，可以容斥。
那么我们用总方案数，减去钦定一些位置为零或一的方案数。
这需要我们对于每个 每位是零，一或无限制 的状态计算方案数。
零的位是简单的，直接弃掉即可；无限制的位置也好办；关键在于一的位置。
一个暴力的想法是，\(f_{i,j,k,t}\) 表示：
- \(i\) 集合为无限制的集合。
- \(j\) 集合为限制一的集合。
- \(k\) 集合为限制一的集合中已经被满足的部分。
- \(t\) 表示一共用了 \(t\) 次机会。
总方案数分析得到是 \(n4^n\)，处在能接受的边缘。
考虑转移。一步转移是枚举 \(j\) 中剔除 \(k\) 后的子集，配上 \(i\) 中的一些东西。
处理出 \(g_{i,v}\) 表示 \(i\) 集合中用不超过 \(v\) 能找到的集合数。则暴力转移可以做到 \(n5^n\)。
考虑做一些优化。首先在 \(i\) 固定的时候，可以处理出 \(h_j\) 表示用一次操作覆盖 \(j\) 中所有点，可行的方案数。
那么事实上我们就要计算 \(h_{j,k}\)，\(j\) 维表示出现与否，\(k\) 维表示用了几次操作，然后合并是 \(j\) 维子集卷积，\(k\) 维做加法。
暴力做 \(k\) 轮子集卷积，可以在 \(O(n^32^n)\) 的时间内解决问题。外层还要枚举 \(i\)，总复杂度 \(n^33^n\)。
另外再扯一嘴容斥。因为一和零的情形，我们钦定它是一和零，它就是一和零，所以我们最后可以把每维分成两类：不知道合不合法的位，以及钦定不合法的位。一个状态对最终答案的系数，是 \(-1\) 的其中钦定不合法位数次幂，可以简单计算。
代码：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
void ADD(int&x,const int&y){if((x+=y)>=mod)x-=mod;}
void ADD(int&x,const ll&y){x=(x+y)%mod;}
class AquaparkPuzzle{
private:
int n,lim,s[1<<11],g[1010],h[12][1<<11],f[12][1<<11],C[12],pc[1<<11],res;
void FMT(int*a,int U,int tp){
for(int i=0;i<n;i++)if(U&(1<<i))for(int u=U;u;u=(u-1)&U)if(u&(1<<i))ADD(a[u],(mod+tp*a[u^(1<<i)])%mod);
}
public:
int countSchedules(int K,int m,vector<int>c){
n=c.size(),lim=1<<n;
for(int i=0;i<lim;i++)for(int j=0;j<n;j++)if(i&(1<<j))s[i]+=c[j],pc[i]++;
for(int _=0;_<lim;_++){//the set of free positions.
// puts("-------------------------------------------------");
// printf("%d:\n",_);
for(int i=_;;i=(i-1)&_){if(s[i]<=m)g[s[i]]++;if(!i)break;}
for(int i=1;i<=m;i++)g[i]+=g[i-1];
C[0]=ksm(g[m],K);
for(int i=0;i<n-pc[_];i++)
C[i+1]=1ll*C[i]*ksm(g[m])%mod*(K-i)%mod*ksm(i+1)%mod;
// printf("C:");for(int i=0;i<=n-pc[_];i++)printf("%d ",C[i]);puts("");
// printf("G:");for(int i=0;i<=m;i++)printf("%d ",g[i]);puts("");
int U=(lim-1)^_;
for(int u=U;u;u=(u-1)&U){if(s[u]<=m)h[pc[u]][u]=g[m-s[u]];else h[pc[u]][u]=0;f[0][u]=1;}
h[0][0]=0,f[0][0]=1;
for(int i=0;i<=n-pc[_];i++)FMT(h[i],U,1);
int now=C[0];
for(int i=0;i<n-pc[_];i++){
for(int u=n-pc[_];u>=i;u--){
for(int v=1;u+v<=n-pc[_];v++)for(int w=U;;w=(w-1)&U){
f[u+v][w]=(1ll*f[u][w]*h[v][w]+f[u+v][w])%mod;
if(!w)break;
}
for(int w=U;;w=(w-1)&U){f[u][w]=0;if(!w)break;}
}
for(int u=i+1;u<=n-pc[_];u++){
FMT(f[u],U,-1);
for(int w=U;;w=(w-1)&U){
if(pc[w]==u)ADD(now,1ll*f[u][w]*C[i+1]);else f[u][w]=0;
if(!w)break;
}
FMT(f[u],U,1);
}
}
for(int i=0;i<=n-pc[_];i++)for(int w=U;;w=(w-1)&U){f[i][w]=h[i][w]=0;if(!w)break;}
for(int i=0;i<=m;i++)g[i]=0;
if((n-pc[_])&1)ADD(res,mod-now);else ADD(res,now);
// printf("%d\n",now);
}
return res;
}
}my;
VI.大度一点
给定一张无向简单图，求有多少边的子集，满足仅保留这些边时，所有边的度数都不小于 \(2\)。
有 \(q\) 次询问，每次给定一个点集，求在该点集的导出子图中做上述计数的结果。
答案对 \(998244353\) 取模。
数据范围：\(n\leq19,q\leq10^5\)。
解法 1.容斥。
钦定若干点的度数小于等于一。记 \(f(\mathbb{S,T})\) 表示在 \(\mathbb S\) 的导出子图中，钦定 \(\mathbb T\) 中节点的度数小于等于一的方案数。则 \(\bb S\) 的答案即为 \(F(\bb S)=\sum\limits_{\bb{T\subseteq S}}(-1)^{|\bb T|}f(\bb{S,T})\)。
现在考虑 \(g\) 应如何计算。被钦定度数小于等于一的节点共可分成三种：
- 度数为零的节点。
- 度数为一，且连向另一个被钦定度数为一的节点的节点。
- 度数为一，且连向未被钦定节点的节点。
我们令 \(g(\bb{S,T})\) 为仅考虑第三种钦定节点的场合。由 \(g\) 到 \(f\) 的过程中，我们可以进行的操作是：
- 选择一个点 \(x\notin\bb S\) 并令其为孤立点，转移到 \((\bb S\cup\{x\},\bb T\cup\{x\})\)，系数为 \(-1\)。
- 选择一条边 \((x,y)\) 且 \(x,y\notin\bb S\) 并令其为孤立边，转移到 \((\bb S\cup\{x,y\},\bb T\cup\{x,y\})\)，系数为 \(1\)。
然后我们发现两种操作与 \(\bb T\) 无关，仅与 \(\bb S\) 有关。于是我们可以先由 \(g\) 求出 \(G\)，然后对 \(G\) 做类似 SOS-DP 的流程，即每次加入一个点或一条边。直接大力实现，复杂度 \(O(n^22^n)\)。
下面我们只需考虑如何求出 \(g\) 了。则，\(g(\bb{S,T})=2^{\bb{|E(S\setminus T)|}}\prod\limits_{x\in\bb T}|\bb E(\bb{S\setminus T},x)|\)。
然后我们发现，在 \(\bb{S\setminus T}\) 固定时，往 \(\bb T\) 中添加一个元素仅仅会让 \(g\) 乘上一位。于是我们用 \(O(2^n)\) 计算子集和的经典技巧，即可简单 \(O(3^n)\) 求出全体 \(g\)，进而求出 \(G\)，进而求出 \(F\)。
复杂度 \(O(3^n+n^22^n)\)。虽然 \(3^{19}\) 看上去很糟糕，但是正解是 \(O(n^32^n)\) 的，实际复杂度甚至还不如 \(2^n\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
void ADD(int&x,int y){if((x+=y)>=mod)x-=mod;}
int n,m,q;
int g[19];
int pc[1<<19],lb[1<<19];
int f[1<<19];
int bin[410];
int h[1<<19];
int d[19];
int main(){
freopen("big.in","r",stdin);
freopen("big.out","w",stdout);
scanf("%d%d",&n,&m);
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),x--,y--,g[x]|=1<<y,g[y]|=1<<x;
for(int i=1;i<(1<<n);i++)pc[i]=pc[i>>1]+(i&1),lb[i]=(i&1?0:lb[i>>1]+1);
bin[0]=1;for(int i=1;i<=n*n;i++)bin[i]=(bin[i-1]<<1)%mod;
for(int i=0;i<(1<<n);i++){
int E=0;
for(int j=0;j<n;j++)
if(i&(1<<j))E+=pc[g[j]&i];
else d[j]=pc[g[j]&i];
E>>=1;
h[0]=bin[E];
for(int U=((1<<n)-1-i),u=U;;u=(u-1)&U){
int x=U-u;
if(x)h[x]=1ll*h[x^(x&-x)]*d[lb[x]]%mod;
// printf("%d,%d|%d\n",i,x,h[x]);
if(pc[x]&1)ADD(f[i|x],mod-h[x]);else ADD(f[i|x],h[x]);
if(!u)break;
}
}
// for(int i=0;i<(1<<n);i++)printf("%d ",f[i]);puts("");
for(int j=0;j<n;j++)for(int i=0;i<(1<<n);i++)if(i&(1<<j))
ADD(f[i],mod-f[i^(1<<j)]);
// for(int i=0;i<(1<<n);i++)printf("%d ",f[i]);puts("");
for(int j=0;j<n;j++)for(int k=0;k<j;k++)if(g[j]&(1<<k))
for(int i=0;i<(1<<n);i++)if((i&(1<<j))&&(i&(1<<k)))
ADD(f[i],f[i^(1<<j)^(1<<k)]);
// for(int i=0;i<(1<<n);i++)printf("%d ",f[i]);puts("");
scanf("%d",&q);
for(int i=1;i<=q;i++){
char s[30];scanf("%s",s);
int x=0;for(int j=0;j<n;j++)if(s[j]=='1')x|=1<<j;
printf("%d\n",f[x]);
}
return 0;
}
VII.集合幂级数的高端操作
本章节特别鸣谢command block神仙的博客。
感觉自己在数数方面一直只能追随cmd神仙的脚步……
在这之前，我们首先定义本章中所有卷积全部为子集卷积，所有变换全部为 or 卷积的 FMT。
接着，我们换一种思路解释子集卷积：
原本我们是将集合幂级数 \(\mathsf F\) 拆做了多个集合幂级数 \(\mathsf F_p\)。但是现在，我们修改定义为将 \(\mathsf F\) 中的每个 \(f(\mathbb S)\)，换作占位形式幂级数 \(f(\mathbb S)x^{|\mathbb S|}\)——虽然该形式幂级数仅有 \(x^{|\mathbb S|}\) 一项非零。明显一样地进行 FMT、点积、FWI，然后仍取 \(x^{|\mathbb S|}\) 一项的系数即可。
- 单位元。
定义集合幂级数关于子集卷积的单位元 \(\epsilon=x^{\varnothing}\)。明显任意级数卷上单位元不变。
- 逆元。
定义 \(\mathsf F\) 的逆元 \(\mathsf F^{-1}\) 满足二者卷积为 \(\epsilon\)。
现在考虑如何求解逆元。考虑对扩展为占位形式幂级数的集合幂级数作 FMT，然后对每一位上的形式幂级数分别求逆，然后再 FMI 回来即可。
求逆也可 \(n^2\) 求。具体就依照求逆的式子算即可。
- \(\exp\)。
存在 \(\exp\) 的级数 \(\mathsf F\) 须保证 \(f(\varnothing)=0\)。
这样之后就定义 \(\exp\mathsf F=\sum\limits_{k=0}\dfrac{\mathsf F^k}{k!}\)。
当然按照定义肯定不好算。不过我们仍然可以 FMT 后每一位分开 \(\exp\)。
考虑如何 \(n^2\) \(\exp\) 。
若 \(e^{f(x)}=g(x)\)，
那么 \(\dfrac{\mathrm de^{f(x)}}{\mathrm dx}=\dfrac{\mathrm dg(x)}{\mathrm dx}\)。
故 \(f'(x)e^{f(x)}=\dfrac{\mathrm dg(x)}{\mathrm dx}\)。
则 \(f'(x)g(x)=g'(x)\)。
那么展开就得到 \(g'_i=\sum\limits_{j+k=i}f'_jg_k\)。
将导数复原即得 \((i+1)g_{i+1}=\sum\limits_{j+k=i}(j+1)f_{j+1}g_k\)。
那么显然就可以递推辣。
- \(\ln\)。
存在 \(\ln\) 的级数 \(\mathsf F\) 须保证 \(f(\varnothing)=1\)。
这样之后就定义 \(\ln \mathsf F=\sum\limits_{k=1}\dfrac{(-1)^{k+1}(\mathsf F-x^{\varnothing})^k}k\)。
同理也应该 FMT 后对每一位分开 \(\ln\)。
考虑如何 \(n^2\) \(\ln\)。
因为 \(\ln\) 是 \(\exp\) 的逆运算，所以我们找到之前的式子 \((i+1)g_{i+1}=\sum\limits_{j+k=i}(j+1)f_{j+1}g_k\)。
则我们要做的就是由 \(g\) 逆推 \(f\)。这是简单的。
具体而言，\((i+1)g_{i+1}=f_{i+1}g_0+\sum\limits_{j=0}^{i-1}(j+1)f_{j+1}g_{i-j}\)。
因为 \(g_0=1\)，所以有 \((i+1)f_{i+1}=(i+1)g_{i+1}-\sum\limits_{j=0}^{i-1}(j+1)f_{j+1}g_{i-j}\)。
改变一下定义就是 \(f_i=g_i-\dfrac1i\sum\limits_{j=1}^{i-1}jf_jg_{i-j}\)。
- 导数。
定义某个集合幂级数 \(x^{\mathbb S}\) 的导数为 \(\sum\limits_{i\in\mathbb S}x^{\mathbb S\setminus i}\)。
同理其具有线性性（即，导数之和等于和之导数）。
因此若 \(\mathsf F=\sum\limits_{\mathbb S}f(\mathbb S)x^{\mathbb S}\)，则 \(\mathsf F'=\sum\limits_{\mathbb S}f(\mathbb S)\sum\limits_{i\in \mathbb S}x^{\mathbb S\setminus i}\)。
- 积分。
计算反导数即可。
- 开根。
FMT 后每一位分开开根。
这个 \(n^2\) 递推是简单的。
I.[NOI Online #3 提高组] 优秀子序列
考虑计数每种集合出现的次数。
然后就会发现答案为 \(\prod\limits_{i=1}^n(x^{\varnothing}+x^{\{a_i\}})\)。
考虑老套路是 \(\ln\) 后 \(\exp\)。
于是 \(\exp\Big(\sum\limits_{i=1}^n\ln(x^{\varnothing}+x^{\{a_i\}})\Big)\)。
按照定义，\(\ln(x^{\varnothing}+x^{\{a_i\}})=\sum\limits_{k=1}\dfrac{(-1)^{k+1}(x^{\{a_i\}})^k}k\)。
而众所周知，\(x^{\{a_i\}}\) 做一次以上子集卷积就会变成 \(0\)。
于是 \(\ln(x^{\varnothing}+x^{\{a_i\}})=x^{\{a_i\}}\)。
于是计算 \(\exp\Big(\sum\limits_{i=1}^nx^{\{a_i\}}\Big)\) 即可。
需要注意的是 \(a_i\) 有可能为 \(0\)。那么这个选不选均可，最后乘上一个 \(2\) 的次幂即可。
代码：
#include<bits/stdc++.h>
using namespace std;
const int N=18;
const int lim=1<<N;
const int mod=1e9+7;
void FMT(int*a,int tp){for(int i=0;i<N;i++)for(int j=0;j<lim;j++)if(j&(1<<i))(a[j]+=(mod+tp*a[j^(1<<i)])%mod)%=mod;}
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
int n,a[19][lim],b[19][lim],res,zer;
#define bit(x) __builtin_popcount(x)
int pri[lim+1],phi[lim+1];
void sieve(){
phi[1]=1;
for(int i=2;i<=lim;i++){
if(!pri[i])pri[++pri[0]]=i,phi[i]=i-1;
for(int j=1;j<=pri[0]&&i*pri[j]<=lim;j++){
pri[i*pri[j]]=true;
if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]];
else{phi[i*pri[j]]=phi[i]*pri[j];break;}
}
}
}
int main(){
scanf("%d",&n),sieve();
for(int i=0,x;i<n;i++){
scanf("%d",&x);
if(!x)zer++;else a[bit(x)][x]++;
}
for(int i=0;i<=N;i++)FMT(a[i],1);
for(int i=0;i<lim;i++){
b[0][i]=1;
for(int j=1;j<=N;j++){
for(int k=0;k<j;k++)(b[j][i]+=1ll*(k+1)*a[k+1][i]%mod*b[j-1-k][i]%mod)%=mod;
b[j][i]=1ll*b[j][i]*ksm(j)%mod;
}
}
for(int i=0;i<=N;i++)FMT(b[i],-1);
for(int i=0;i<lim;i++)(res+=1ll*b[bit(i)][i]*phi[i+1]%mod)%=mod;
while(zer--)(res<<=1)%=mod;
printf("%d\n",res);
return 0;
}
II.[LOJ#154]集合划分计数
有标号集合划分计数等价于 \(\exp\)，这可以参照上一题。
还有一种证明方法就是使用 \(\exp\) 的 \(n^2\) 计算式 \((i+1)g_{i+1}=\sum\limits_{j+k=i}(j+1)f_{j+1}g_k\)。
因为 \(f_i\) 中必然包含 \(i\) 个元素，而这些元素本应该是有顺序、但此处却是无序的（明显集合内部的元素间无序），所以只需要将 \(f,g\) 乘上一个 \(i!\) 就能得到有顺序的OGF。不妨设乘完 \(i!\) 的函数为 \(F,G\)。
故 \(\dfrac{G_{i+1}}{i!}=\sum\limits_{j+k=i}\dfrac{F_{j+1}}{j!}\dfrac{G_k}{k!}\)，即 \(G_{i+1}=\sum\limits_{j+k=i}\dbinom ijF_{j+1}G_k\)。
这个式子非常有哲理，因为其组合意义就是枚举 \(i+1\) 号点与哪些点分在一块。
不超过 \(k\) 个元素的有标号集合划分计数等价于 \(g(x)=\sum\limits_{i=0}^k\dfrac{f^i(x)}{i!}\)——这是显然的。
考虑计算。首先将两边同时求导得到 \(g'=f'(x)\times\sum\limits_{i=0}^{k-1}\dfrac{f^i(x)}{i!}\)。
然后明显右式等于 \(f'(g-f^k/k!)\)。可以简单计算。
复杂度 \(n^22^n\)——可以 \(\exp+\ln\) 计算 \(f^k\)。
人傻常数大，需要狠优化。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int n,m,K,f[22][1<<21],g[22],lim,fac[22],inv[22],INV[22],h[22],H[22],res;
#define bit(x) __builtin_popcount(x)
void ADD(int&x,int y){if((x+=y)>=mod)x-=mod;}
void FMT(int*a){for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))ADD(a[j],a[j^(1<<i)]);}
int ksm(int x,int y=mod-2){int z=1;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)z=1ll*z*x%mod;return z;}
void kexp(int i){
memset(h,0,sizeof(h));
int t=0;while(t<=n&&!f[t][i])t++;
if(t*K>=n)return;
int val=f[t][i],iv=ksm(val),kv=ksm(val,K);
for(int j=t;j<=n;j++)f[j][i]=1ll*f[j][i]*iv%mod;
memset(H,0,sizeof(H));
for(int j=1;j+t<=n;j++){
H[j]=1ll*f[j+t][i]*j%mod;
for(int k=1;k<j;k++)ADD(H[j],mod-1ll*H[k]*f[t+j-k][i]%mod);
}
for(int j=t;j<=n;j++)f[j][i]=1ll*f[j][i]*val%mod;
for(int j=1;j+t<=n;j++)H[j]=1ll*H[j]*K%mod;
h[0]=1;
for(int j=1;j+t<=n;j++){
for(int k=1;k<=j;k++)ADD(h[j],1ll*H[k]*h[j-k]%mod);
h[j]=1ll*h[j]*INV[j]%mod;
}
for(int j=n;j>=t*K;j--)h[j]=h[j-t*K];
for(int j=0;j<t*K;j++)h[j]=0;
for(int j=t*K;j<=n;j++)h[j]=1ll*h[j]*kv%mod;
}
int read(){
int x=0;
char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x;
}
int main(){
// freopen("cspartition_5_1.in","r",stdin);
lim=1<<(n=read()),m=read(),K=read();
fac[0]=1;for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;
inv[n]=ksm(fac[n]);for(int i=n;i;i--)inv[i-1]=1ll*inv[i]*i%mod;
for(int i=1;i<=n;i++)INV[i]=1ll*inv[i]*fac[i-1]%mod;
for(int i=0,x;i<m;i++)x=read(),f[bit(x)][x]++;
for(int i=0;i<=n;i++)FMT(f[i]);
for(int i=0;i<lim;i++){
kexp(i);
for(int j=0;j<=n;j++)h[j]=1ll*h[j]*inv[K]%mod,f[j][i]=1ll*f[j][i]*j%mod;
memset(g,0,sizeof(g)),g[0]=1;
for(int j=1;j<=n;j++){
for(int k=1;k<=j;k++)ADD(g[j],1ll*f[k][i]*(g[j-k]+mod-h[j-k])%mod);
g[j]=1ll*g[j]*INV[j]%mod;
}
if((n-bit(i))&1)ADD(res,mod-g[n]);else ADD(res,g[n]);
}
printf("%d\n",res);
return 0;
}
III.[UOJ#94]【集训队互测2015】胡策的统计
考虑令集合幂级数 \(\mathsf F\) 表示点集为 \(\mathbb S\) 时的联通子图个数，\(\mathsf G\) 则为全部子图个数。
显然 \(\mathsf G\) 是好求的（\(2^{\text{内部边数}}\)）。
依照上文结论，有 \(\exp\mathsf F=\mathsf G\)。那么 \(\mathsf F=\ln\mathsf G\)。
考虑答案。显然 \(k\) 个连通块的子图数即为 \(\dfrac{\mathsf F^k}{k!}\)。则答案为 \(\sum\dfrac{\mathsf F^k}{k!}\times k!=\dfrac{1}{1-\mathsf F}\)。
这是易解决的。
我卡这倒霉玩意卡了很久，最后循环顺序一改，3s+瞬间变1s+。？？？
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
typedef long long ll;
const int N=1050000;
int n,m,pov[410],lim,f[N][22],INV[22],s[N],bit[N];
ll res;
int SUM(register int x,register int y){return x+y>=mod?x+y-mod:x+y;}
void ADD(register int&x,register int y){if((x+=y)>=mod)x-=mod;}
int read(){
register int x=0;
register char c=getchar();
while(c>'9'||c<'0')c=getchar();
while(c>='0'&&c<='9')x=(x<<3)+(x<<1)+(c^48),c=getchar();
return x;
}
int main(){
lim=1<<(n=read()),m=read();
pov[0]=1;for(register int i=1;i<=m;i++)pov[i]=(pov[i-1]<<1)%mod,s[(1<<(read()-1))+(1<<(read()-1))]++;
for(register int i=0;i<n;i++)for(register int j=0;j<lim;j++)if(j&(1<<i))s[j]+=s[j^(1<<i)];
for(register int i=0;i<lim;i++)bit[i]=bit[i>>1]+(i&1);
for(register int i=0;i<lim;i++)f[i][bit[i]]=pov[s[i]];
for(register int i=0;i<n;i++)for(register int j=0;j<lim;j++)if(j&(1<<i))for(register int k=1;k<=n;k++)ADD(f[j][k],f[j^(1<<i)][k]);
INV[1]=1;for(register int i=2;i<=n;i++)INV[i]=1ll*INV[mod%i]*(mod-mod/i)%mod;
for(register int i=0;i<lim;i++){
static ll g[22],h[22];
// for(int j=0;j<=n;j++)printf("%d ",f[j][i]);puts("");
for(register int j=0;j<=n;j++)h[j]=f[i][j],g[j]=0;
for(register int j=1;j<=n;j++){for(register int k=1;k<j;k++)g[j]=(g[j]+g[k]*h[j-k])%mod;g[j]=(h[j]*j+mod-g[j])%mod;}
for(register int j=1;j<=n;j++)g[j]=g[j]*INV[j]%mod;
// for(int j=0;j<=n;j++)printf("%d ",g[j][i]);puts("");
for(register int j=1;j<=n;j++)g[j]=mod-g[j];
g[0]=h[0]=1;
for(register int j=1;j<=n;j++){h[j]=0;for(register int k=0;k<j;k++)h[j]=(h[j]+h[k]*g[j-k])%mod;h[j]=mod-h[j];}
if((n-bit[i])&1)res-=h[n];else res+=h[n];
}
printf("%d\n",(res%mod+mod)%mod);
return 0;
}
IV.[LOJ#3165][CEOI2019]游乐园
Observation 1.对于一种合法方案，将所有边取反后与另一种合法状态一一对应。且对应的两者对答案的贡献之和定为边数。
那么我们只需要计数所有的合法方案数，然后乘以边数再除以 \(2\) 即可。
现在我们考虑设计DP来计数。
一个显然的想法是设 \(f(\mathbb S)\) 表示 \(\mathbb S\) 集合内部合法方案数。
现在考虑转移。
发现从中剪除所有入度为 \(0\)（出度亦可，反正边的方向全翻转后也没人管）的点会使得计数不重不漏。
发现，只要是在集合中的一组独立集（即内部无边的集）都可能成为上述剪除的点集。
但是为了避免被重复计算必须容斥一下。
于是转移就有了：\(f(\mathbb S)=\sum\limits_{(\mathbb{T\subseteq S})\land(\mathbb T\text{是独立集})\land(\mathbb T\neq\varnothing)}(-1)^{|\mathbb T|-1}f(\mathbb S\setminus\mathbb T)\)。
从另一种方面来说就是，如果设 \(g(\mathbb T)=[\mathbb T\text{是独立集}][\mathbb T\neq\varnothing](-1)^{|\mathbb T|-1}\)，那么 \(g\) 与 \(f\) 似乎可以通过卷积之类搞出来。
那不就好办了嘛。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
const int inv2=499122177;
int n,m,lim,f[1<<18][19],g[1<<18][19],res;
bool s[1<<18];
#define bit(x) __builtin_popcount(x)
int main(){
scanf("%d%d",&n,&m),lim=1<<n;
for(int i=1,x,y;i<=m;i++)scanf("%d%d",&x,&y),s[(1<<(x-1))|(1<<(y-1))]=true;
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))s[j]|=s[j^(1<<i)];
for(int i=1;i<lim;i++){
if(s[i])continue;
if(bit(i)&1)f[i][bit(i)]=1;
else f[i][bit(i)]=mod-1;
}
for(int j=0;j<n;j++)for(int k=0;k<lim;k++)if(k&(1<<j))for(int i=1;i<=n;i++)(f[k][i]+=f[k^(1<<j)][i])%=mod;
for(int i=0;i<lim;i++){
g[i][0]=1;
for(int j=0;j<=n;j++)for(int k=1;j+k<=n;k++)(g[i][j+k]+=1ll*g[i][j]*f[i][k]%mod)%=mod;
if((n-bit(i))&1)(res+=mod-g[i][n])%=mod;else(res+=g[i][n])%=mod;
}
printf("%d\n",1ll*res*m%mod*inv2%mod);
return 0;
}
V.[LOJ#6673]EntropyIncreaser 与山林
题目要求出欧拉连通子图数。
那么如果我们求出欧拉子图数，对其求 \(\ln\) 即可得到联通子图。
于是问题转换为对每个子点集求出其欧拉子图数。
欧拉子图的数数方法：
考虑求出该点集的导出子图（即所有两端都在该点集之内的边构成的子图）的任一生成森林，并任意选择或不选剩下不在森林内的边。
现在考虑决定森林中的边选择与否。可以发现，使用类树形DP的方案，我们可以对每一组其它边的选择方案唯一确定一组森林边选择方案，使得所有点的度数均为偶。
故欧拉子图的数量即为 \(2^{\text{边数-森林边数}}\)。
代码：
#include<bits/stdc++.h>
using namespace std;
const int mod=998244353;
int n,m,lim,X[410],Y[410],f[1<<20][21],dsu[20],pov[410],INV[21],g[21],res;
int find(int x){return dsu[x]==x?x:dsu[x]=find(dsu[x]);}
bool merge(int x,int y){if((x=find(x))==(y=find(y)))return false;dsu[x]=y;return true;}
#define bit(x) __builtin_popcount(x)
int main(){
scanf("%d%d",&n,&m),lim=1<<n;
pov[0]=1;for(int i=1;i<=m;i++)scanf("%d%d",&X[i],&Y[i]),X[i]--,Y[i]--,pov[i]=(pov[i-1]<<1)%mod;
for(int i=0;i<lim;i++){
for(int j=0;j<n;j++)if(i&(1<<j))dsu[j]=j;
int sum=0;
for(int j=1;j<=m;j++)if((i&(1<<X[j]))&&(i&(1<<Y[j])))sum+=!merge(X[j],Y[j]);
// printf("%d:%d\n",i,sum);
f[i][bit(i)]=pov[sum];
}
INV[1]=1;for(int i=2;i<=n;i++)INV[i]=1ll*INV[mod%i]*(mod-mod/i)%mod;
for(int i=0;i<n;i++)for(int j=0;j<lim;j++)if(j&(1<<i))for(int k=0;k<=n;k++)(f[j][k]+=f[j^(1<<i)][k])%=mod;
for(int i=0;i<lim;i++){
g[0]=0;
for(int j=1;j<=n;j++){
g[j]=0;
for(int k=1;k<j;k++)(g[j]+=1ll*g[k]*f[i][j-k]%mod)%=mod;
g[j]=(1ll*j*f[i][j]+mod-g[j])%mod;
}
for(int j=1;j<=n;j++)g[j]=1ll*g[j]*INV[j]%mod;
if((n-bit(i))&1)(res+=mod-g[n])%=mod;else(res+=g[n])%=mod;
}
printf("%d\n",res);
return 0;
}
浙公网安备 33010602011771号