位运算卷积(FWT) & 集合幂级数
command_block
2019-05-04 23:24:57
$\rm2020.10:$ 进行了小修小补。
$\rm2020.11:$ 重制了更严谨的 “k进制不进位加法卷积” 部分。
$\rm2020.11:$ 加更集合幂级数。
# 0x01.$\rm FWT$ 概论
- **位运算卷积**
众所周知，多项式乘法是加法卷积，因为第 $i$ 项和第 $j$ 项的乘积贡献到第 $i+j$ 项。
类似地定义位运算卷积 : 第 $i$ 项和第 $j$ 项的乘积贡献到第 $i⊕j$ 项。其中 $⊕$ 是某种**位运算**。
即 $S[k]=\sum\limits_{i⊕j=k}A[i]B[j]$ ，记作卷积式 $A*B=S$ 。
- **构造的尝试**
众所周知，$\rm FFT$ 把多项式转换成点值之后，从卷积变为了直接点积。
我们自然也期望把位运算卷积转化成点积。
设 $FWT(A)$ 是幂级数 $A$ 经过 $\rm FWT$ 变换之后得到的幂级数。
我们需要令其满足 : $A*B=C \Longleftrightarrow FWT(A)·FWT(B)=FWT(C)$ (点积)。
$\rm FFT$ 是一个**线性变换**，我们也希望 $\rm FWT$ 变换是线性的。
我们还不知道怎么变换，于是设 $c(i,j)$ 为变换系数，即 $A[j]$ 对 $FWT(A)[i]$ 的贡献系数。
则 $FWT(A)[i]=\sum\limits_{j=0}^{n-1}c(i,j)A_j$
根据 $FWT(A)·FWT(B)=FWT(C)$ ,得到 :
$$FWT(A)[i]FWT(B)[i]=FWT(C)[i]$$
$$\sum\limits_{j=0}^{n-1}c(i,j)A[j]\sum\limits_{k=0}^{n-1}c(i,k)B[k]=\sum\limits_{p=0}^{n-1}c(i,p)C[p]$$
$$\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum\limits_{p=0}^{n-1}c(i,p)C[p]$$
根据 $A*B=C$ ,又能得到 :
$$C[p]=\sum\limits_{t_1⊕t_2=p}A[t_1]B[t_2]$$
$$\sum\limits_{p=0}^{n-1}c(i,p)C[p]=\sum\limits_{p=0}^{n-1}c(i,p)\sum\limits_{t_1⊕t_2=p}A[t_1]B[t_2]$$
$$\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum\limits_{p=0}^{n-1}c(i,p)\sum\limits_{t_1⊕t_2=p}A[t_1]B[t_2]$$
$$\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}c(i,j)c(i,k)A[j]B[k]=\sum\limits_{p=0}^{n-1}\sum\limits_{t_1⊕t_2=p}A[t_1]B[t_2]c(i,t_1⊕t_2)=\sum\limits_{t_1=0}^{n-1}\sum\limits_{t_2=0}^{n-1}A[t_1]B[t_2]c(i,t_1⊕t_2)$$
对比左右两边，得到 $c(i,j)c(i,k)=c(i,j⊕k)$ ，只需要满足这个就好了。
另外,由于位运算每一位的独立性, $c(i,j)$ 也有一个重要性质: 可以分位考虑。
假设我们已经(根据真值表)求出满足要求的 $c([0,1],[0,1])$，我们能这样构造出所有的 $c$ :
设二进制数 $a$ 的每一位分别为 : $a_0,a_1,a_2...$
则有 $c(i,j)=c(i_0,j_0)c(i_1,j_1)c(i_2,j_2)...$ ，就是**把每一位的变换系数乘起来**。
那么 : 对于每个 $t$ ，都有 $c(i_t,j_t)c(i_t,k_t)=c(i_t,j_t⊕k_t)\Longleftrightarrow c(i,j)c(i,k)=c(i,j⊕k)$
这就符合我们的条件。
好了,现在假设我们已经有了符合要求的 $c$ ,如何快速求解 $\rm FWT$ 变换呢?
$FWT(A)[i]=\sum\limits_{j=0}^{n-1}c(i,j)A[j]$
根据这个式子直接求和,复杂度至少是 $O(n^2)$ 的,并没有起到优化作用。
我们考虑**按位拆半**：
$=\sum\limits_{j=0}^{(n/2)-1}c(i,j)A[j]+\sum\limits_{j=(n/2)}^{n-1}c(i,j)A[j]$
设 $i'$ 为 $i$ 去掉二进制首位剩下的数。
在首位分开考虑 :
$=\sum\limits_{j=0}^{(n/2)-1}c(i_0,j_0)c(i',j')A[j]+\sum\limits_{j=(n/2)}^{n-1}c(i_0,j_0)c(i',j')A[j]$
$=c(i_0,0)\sum\limits_{j=0}^{(n/2)-1}c(i',j')A[j]+c(i_0,1)\sum\limits_{j=(n/2)}^{n-1}c(i',j')A[j]$
考虑到 $c(i',j')$ 就是去掉最高位的子变换,这里规模减半了。
设 $A_0$ 为幂级数下标首位为 $0$ 的部分,类似地有 $A_1$。
若 $i_0=0$ ,则有 :
$FWT(A)[i]=c(0,0)FWT(A_0)[i]+c(0,1)FWT(A_1)[i]\quad \big(i\in[0,n/2)\big)$
若 $i_0=1$ ,则有 :
$FWT(A)[i+(n/2)]=c(1,0)FWT(A_0)[i]+c(1,1)FWT(A_1)[i]\quad \big(i\in[0,n/2)\big)$
我们就能以 $O(n)$ 的代价,根据上列式子合并两个规模为 $n/2$ 的子变换。
所以,若 $n=2^m$ ,需要合并 $m$ 次,复杂度为 $O(m2^m)$。
( 可能有点抽象,但是您如果写过FFT,看到代码就会懂了 )
此外,逆变换 $\rm IFWT$ 就是对 $c$ 矩阵求个逆,具体见下文。
( 一个重要的地方是,这个构造出来的 $c$ 矩阵**一定要有逆**,否则就变换不回去TAT )
# 0x02.基础位运算卷积
针对不同的位运算,根据 $c(i,j)c(i,k)=c(i,j⊕k)$ 构造出 $c([0,1],[0,1])$ 即可。
我们把这个矩阵称为**位矩阵**。
构造的过程可能有些繁杂，可以直接记结论，或者去后面看扩展版的。
## $1.1\ \rm Or$ 卷积
设位矩阵为 $c=\begin{bmatrix}c(0,0)&c(0,1)\\c(1,0)&c(1,1)\end{bmatrix}$
**起点** : $c(i,j)c(i,k)=c(i,j|k)$
- $c(0,0)c(0,0)=c(0,0|0)$
$\Rightarrow c(0,0)=1$ 或 $0$。
同理不难推知 $c(\_,\_)∈\{0,1\}$
- $c(0,1)c(0,0)=c(0,1|0)$
$\Rightarrow c(0,1)=0$ 或 $c(0,0)=c(0,1)=1$
- $c(1,1)c(1,0)=c(1,1|0)$
$\Rightarrow c(1,1)=0$ 或 $c(1,0)=c(1,1)=1$
首先,如果有一排0或者一列0那么这个矩阵就没有逆,那么可以构造出:
两种情况 : $\begin{bmatrix}1&1\\1&0\end{bmatrix}$ 或 $\begin{bmatrix}1&0\\1&1\end{bmatrix}$。
> **Tips** : $Or$ 卷积的上面第二个矩阵 $FWT$ 相当于子集求和。
> 原因:第二个矩阵相当于 $c(i,j)=[i\&j=j]$
> $A'_i=\sum\limits_{i\&j=j}A_i$等价于$A'_i=\sum\limits_{j∈i}A_i$。
> (也可以使用高维前缀和来推导)
(下面采用第二个矩阵)
$FWT(A)[i]=FWT(A_0)[i]$
$FWT(A)[i+(n/2)]=FWT(A_0)[i]+FWT(A_1)[i]$
对于逆变换,把矩阵求个逆可得 $\begin{bmatrix}1&0\\-1&1\end{bmatrix}$ 。
$IFWT(A)[i]=IFWT(A_0)[i]$
$IFWT(A)[i+(n/2)]=IFWT(A_1)[i]-IFWT(A_0)[i]$
## $1.2\ \rm And$ 卷积
**起点** : $c(i,j)c(i,k)=c(i,j\&k)$。
同上,容易得到 $c(\_,\_)∈\{0,1\}$。
- $c(0,1)c(0,0)=c(0,1\&0)$
$\Rightarrow c(0,0)=0$ 或 $c(0,0)=c(0,1)=1$
- $c(1,1)c(1,0)=c(1,1\&0)$
$\Rightarrow c(1,0)=0$ 或 $c(1,0)=c(1,1)=1$
还是老样子,如果有一排 $0$ 或者一列 $0$ 那么这个矩阵就没有逆,那么可以构造出:
两种情况 : $\begin{bmatrix}0&1\\1&1\end{bmatrix}$ 或 $\begin{bmatrix}1&1\\0&1\end{bmatrix}$,下面采用第二种。
$FWT(A)[i]=FWT(A_0)[i]+FWT(A_1)[i]$
$FWT(A)[i+(n/2)]=FWT(A_1)[i]$
把矩阵求个逆可得$\begin{bmatrix}1&-1\\0&1\end{bmatrix}$:
$IFWT(A)[i]=IFWT(A_0)[i]-IFWT(A_1)[i]$
$IFWT(A)[i+(n/2)]=IFWT(A_1)[i]$
## $1.3\ \rm Xor$ 卷积
**起点** : $c(i,j)c(i,k)=c(i,j\ \text{xor}\ k)$
- 对于任意的 $x,y$ ,均有 $c(0,0)c(x,y)=c(x,y\ \text{xor}\ 0)=c(x,y)$
$\Rightarrow c(0,0)=1$.
- $c(1,1)c(1,1)=c(1,0)$
此时若 $c(1,1)=c(1,0)=0$,则一行为 $0$ ,矩阵无逆。
所以 $c(1,1),c(1,0)$ 必然都非 $0$。
- $c(1,0)c(1,1)=c(1,1)$
刚才说$c(1,1)$ 非 $0$,所以此处 $c(1,0)$ 一定是1.
- $c(0,1)c(0,1)=c(0,0)$
$\Rightarrow c(0,1)∈\{-1,1\}$
两种情况 : $\begin{bmatrix}1&1\\-1&1\end{bmatrix}$ 或 $\begin{bmatrix}1&1\\1&-1\end{bmatrix}$ ,下面采用第二种。
附 ：不难观察出 $c(i,j)=(-1)^{|i\&j|}$
$FWT(A)_i=FWT(A_0)_i+FWT(A_1)_i$
$FWT(A)_{i+(n/2)}=FWT(A_0)_i-FWT(A_1)_{i}$
求逆可得$\begin{bmatrix}0.5&0.5\\0.5&-0.5\end{bmatrix}$
$IFWT(A)_i=\dfrac{IFWT(A_0)_i+IFWT(A_1)_i}{2}$
$IFWT(A)_{i+(n/2)}=\dfrac{IFWT(A_0)_i-IFWT(A_1)_i}{2}$
## $1.4$ 模板题 & Code:
[P4717 【模板】快速沃尔什变换 (FWT)](https://www.luogu.com.cn/problem/P4717)
- **FWT** :
[评测记录](https://www.luogu.com.cn/record/34249999)
```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#define Maxn 135000
#define ll long long
using namespace std;
const int mod=998244353,inv2=499122177;
inline int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
const ll
Cor[2][2] ={{1,0},{1,1}},
Cand[2][2] ={{1,1},{0,1}},
Cxor[2][2] ={{1,1},{1,mod-1}},
ICor[2][2] ={{1,0},{mod-1,1}},
ICand[2][2]={{1,mod-1},{0,1}},
ICxor[2][2]={{inv2,inv2},{inv2,mod-inv2}};
void FWT(ll *F,const ll c[2][2],int n)
{
for (int len=1;len<n;len<<=1)
for (int p=0;p<n;p+=len+len)
for (int i=p;i<p+len;i++){
ll sav0=F[i];
F[i]=(c[0][0]*F[i]+c[0][1]*F[i+len])%mod;
F[i+len]=(c[1][0]*sav0+c[1][1]*F[i+len])%mod;
}
}
void bitmul(ll *F,ll *G,const ll C[2][2],const ll IC[2][2],int n)
{
FWT(F,C,n);FWT(G,C,n);
for (int i=0;i<n;i++)F[i]=F[i]*G[i]%mod;
FWT(F,IC,n);
}
void print(ll *s,int n){
for (int i=0;i<n;i++)
printf("%lld ",s[i]);puts("");
}
#define cpy(f,g,n) memcpy(f,g,sizeof(ll)*(n))
ll f[Maxn],g[Maxn],a[Maxn],b[Maxn];
int n;
int main()
{
n=(1<<read());
for (int i=0;i<n;i++)f[i]=read();
for (int i=0;i<n;i++)g[i]=read();
cpy(a,f,n);cpy(b,g,n);
bitmul(a,b,Cor,ICor,n);
print(a,n);
cpy(a,f,n);cpy(b,g,n);
bitmul(a,b,Cand,ICand,n);
print(a,n);
cpy(a,f,n);cpy(b,g,n);
bitmul(a,b,Cxor,ICxor,n);
print(a,n);
return 0;
}
```
- **分治乘法** :
留个代码以备将来研究……
与前后文无关。
```cpp
#include<algorithm>
#include<cstdio>
#define Maxn 140000
#define mod 998244353
#define ll long long
using namespace std;
inline int read()
{
register int X=0;
register char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
int n,pn,inv2;
ll f[Maxn],g[Maxn];
ll a[Maxn],b[Maxn],c[Maxn];
void mulor(ll *a,ll *b,ll *c,int len)
{
if (!(len>>=1)){
c[0]=(a[0]*b[0])%mod;
return ;
}for (int i=0;i<len;i++){
a[i+len]=(a[i+len]+a[i])%mod;
b[i+len]=(b[i+len]+b[i])%mod;
}mulor(a,b,c,len);
mulor(a+len,b+len,c+len,len);
for (int i=0;i<len;i++)
c[i+len]=(c[i+len]-c[i]+mod)%mod;
}
void muland(ll *a,ll *b,ll *c,int len)
{
if (!(len>>=1)){
c[0]=(a[0]*b[0])%mod;
return ;
}for (int i=0;i<len;i++){
a[i]=(a[i+len]+a[i])%mod;
b[i]=(b[i+len]+b[i])%mod;
}muland(a,b,c,len);
muland(a+len,b+len,c+len,len);
for (int i=0;i<len;i++)
c[i]=(c[i]-c[i+len]+mod)%mod;
}
void mulxor(ll *a,ll *b,ll *c,int len)
{
if (!(len>>=1)){
c[0]=(a[0]*b[0])%mod;
return ;
}for (int i=0;i<len;i++){
ll sava=a[i],savb=b[i];
a[i]=(a[i+len]-a[i]+mod)%mod;
b[i]=(b[i+len]-b[i]+mod)%mod;
a[i+len]=(a[i+len]+sava)%mod;
b[i+len]=(b[i+len]+savb)%mod;
}mulxor(a,b,c,len);
mulxor(a+len,b+len,c+len,len);
for (int i=0;i<len;i++){
ll savc=c[i];
c[i]=(savc+c[i+len])*inv2%mod;
c[i+len]=(c[i+len]-savc+mod)*inv2%mod;
}
}
void print(ll *s)
{
for (int i=0;i<n;i++)
printf("%lld ",s[i]);puts("");
}
void copy(ll *f,ll *g)
{for (int i=0;i<n;i++)f[i]=g[i];}
int main()
{
inv2=(mod+1)/2;
scanf("%d",&n);n=(1<<n);
for (int i=0;i<n;i++)f[i]=read();
for (int i=0;i<n;i++)g[i]=read();
copy(a,f);copy(b,g);
mulor(a,b,c,n);
print(c);
copy(a,f);copy(b,g);
muland(a,b,c,n);
print(c);
copy(a,f);copy(b,g);
mulxor(a,b,c,n);
print(c);
return 0;
}
```
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
# 0x03.$\rm FWT$ 变换的一些性质
前面讲了, $\rm FWT$ 变换的本质是线性变换。
这意味着 $FWT(A+B)=FWT(A)+FWT(B)$ ,且 $FWT(cA)=cFWT(A)$
此外,如果需要卷积的向量只有少数非 $0$ 项( $2,3$ 个)可能会有分类讨论的神奇的解法。
例题(~~我只能做到这些了~~):
[CF449D Jzzhu and Numbers](https://www.luogu.org/problem/CF449D) + [题解Link](https://www.luogu.org/blog/command-block/post-shuo-xue-ji-lu-cf449d-jzzhu-and-numbers)
[CF1119H Triple](https://www.luogu.org/problem/CF1119H) + [题解Link](https://www.luogu.org/problemnew/solution/CF1119H)
题做多了再来填坑吧。
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
# 0x04.位值域的扩展
其实**位运算的本质**是对一个 $n$ 维 $01$ 向量的运算。
或运算就是每一维取 $\max$。且运算就是每一维取 $min$。
异或运算则是每一维对应相加再 $\bmod\ 2$。
位运算有个特点 : 向量的每一位都是独立的。
我们把$\{0,1\}$扩展到 $[0,k)∩Z$ 也就是扩展到 $k$ 进制，看看会得到什么？
## 每一维取 $\bf max$
对应原来的按位 $\rm or$。
可得 $c(x,y)c(x,y)=c(x,max(y,y))=c(x,y)$
所以整个矩阵中只有 $0$ 或 $1$。
又可得 $c(x,y1)c(x,y2)=c(x,max(y1,y2))$
假如有 $c(x,y1)=1,\ c(x,y2)=0$,可得 $c(x,y1)c(x,y2)=c(x,max(y1,y2))$
又因为 $1*0=0$ 所以 $max(y1,y2)=y2$ ,可得 $y1<y2$.
也就是说，每一行以内， $1$ 总是在 $0$ 的前面。
接下来，除了矩阵有逆之外,再没有别的要求了。
如果要有逆的话，每一行都必须不同，那么 $1$ 的个数分别就是 $1...n$ ,随意排列都可以。
例子 ： $k=4$ 的情形
$\begin{bmatrix}1&0&0&0\\1&1&0&0\\1&1&1&0\\1&1&1&1\end{bmatrix}$或者$\begin{bmatrix}1&1&1&1\\1&0&0&0\\1&1&1&0\\1&1&0&0\end{bmatrix}$等等$4!$种。
为了方便,一般取用第一种,求逆可得$\begin{bmatrix}1&0&0&0\\-1&1&0&0\\0&-1&1&0\\0&0&-1&1\end{bmatrix}$
暴力做的话 $O(k^{n+1}n)$。
这个矩阵本质上就是前缀和和差分……所以可以优化到 $O(k^nn)$。
(宏观上就是高维前缀和)
## 每一维取 $\bf min$
对应原来的按位 $\rm and$。
类似的，整个矩阵中只有 $0$ 或 $1$ ,每一行以内, $1$ 总是在 $0$ 的后面。
例子 ： $k=4$ 的情形
$\begin{bmatrix}1&1&1&1\\0&1&1&1\\0&0&1&1\\0&0&0&1\end{bmatrix}$或者$\begin{bmatrix}1&1&1&1\\0&0&1&1\\0&1&1&1\\0&0&0&1\end{bmatrix}$等等$4!$种。
为了方便，一般取用第一种,求逆可得$\begin{bmatrix}1&-1&0&0\\0&1&-1&0\\0&0&1&-1\\0&0&0&1\end{bmatrix}$
这个矩阵本质上就是后缀和和差分……同样能优化到 $O(k^nn)$
(宏观上就是高维后缀和)
## 每一维加起来 $\bf mod\ k$
对应原来的按位 $\rm xor$。
这玩意就比较复杂了,需要动用单位根……
我们知道 $c(x,y1)c(x,y2)=c(x,(y1+y2)mod\ k)$
$c(x,y)=\omega_k^y$ 就能满足要求： $c(x,y1)c(x,y2)=\omega_k^{y1}\omega_k^{y2}=\omega_k^{(y1+y2)mod\ k}$
但是，每一行都一样的话，矩阵就没有逆。
你会发现直接把 `FFT` 的范德蒙德矩阵拿来用就好了:
$\begin{bmatrix}1& 1 & 1& ... & 1\\ 1& w_k^1& w_k^2& ... & w_k^{k - 1}\\ 1& w_k^2 & w_k^4& ... & w_k^{2(k - 1)}\\ ...& ...& ...& ...& ...\\ 1& w_k^{k - 1}& w_k^{2(k - 1)} & ... & w_k^{(k - 1)(k - 1)}\end{bmatrix}$
逆矩阵我们也知道，那就是:
$\dfrac{1}{k}\begin{bmatrix}1& 1 & 1& ... & 1\\ 1& w_k^{-1}& w_k^{-2}& ... & w_k^{-(k - 1)}\\ 1& w_k^{-2} & w_k^{-4}& ... & w_k^{-2(k - 1)}\\ ...& ...& ...& ...& ...\\ 1& w_k^{-(k - 1)}& w_k^{-2(k - 1)} & ... & w_k^{-(k - 1)(k - 1)}\end{bmatrix}$
暴力计算显然是 $O(k^{n+1}n)$ 的。
可以使用 $\rm FFT$ 优化到 $O(k^nn\log k)$。
但是，单位根在模意义很可能不存在。考虑扩充我们“数”的表示。
首先想到的是，人为地定义代数对象 $x$ 满足 $x^k=1$ （且 $x$ 的阶恰为 $k$），用其代替单位根 $w_k$，不难验证其满足前文构造矩阵的条件。
此举相当于用 $\pmod{x^{k}-1}$ 下的循环卷积多项式替换了“数”。
我们的位矩阵长这样 ：
$C_1=\begin{bmatrix}1& 1 & 1& ... & 1\\ 1& x^1& x^2& ... & x^{k - 1}\\ 1& x^2 & x^4& ... & x^{2(k - 1)}\\ ...& ...& ...& ...& ...\\ 1& x^{k - 1}& x^{2(k - 1)} & ... & x^{(k - 1)(k - 1)}\end{bmatrix}$
$C_2=\dfrac{1}{k}\begin{bmatrix}1& 1 & 1& ... & 1\\ 1& x^{-1}& x^{-2}& ... & x^{-(k - 1)}\\ 1& x^{-2} & x^{-4}& ... & x^{-2(k - 1)}\\ ...& ...& ...& ...& ...\\ 1& x^{-(k - 1)}& x^{-2(k - 1)} & ... & x^{-(k - 1)(k - 1)}\end{bmatrix}$
接下来的推导可能需要抽象代数和高等代数相关知识。详情可见 ： [近世代数乱编](https://www.luogu.com.cn/blog/command-block/jin-shi-dai-shuo-luan-bian)
问题在于，此时 $\pmod{x^{k}-1}$ 下的循环卷积多项式并不是域，其存在零因子。
这样就会导致两个位矩阵并非严格互逆。
具体而言，我们原本希望 $C_1*C_2=I$ ,这需要满足 $\sum\limits_{i=0}^{k-1}C_1[a][i]C_2[i][b]=[a=b]$。
当 $a=b$ 时，有 $\sum\limits_{i=0}^{n-1}C_1[a][i]C_2[i][a]=\dfrac{1}{k}\sum\limits_{i=0}^{k-1}x^{ai}x^{-ai}=1$ ，这不需要消去律也成立，这部分贡献是没有差错的。
当 $a≠b$ 时，有 $\sum\limits_{i=0}^{n-1}C_1[a][i]C_2[i][b]=\dfrac{1}{k}\sum\limits_{i=0}^{k-1}x^{ai}x^{-bi}=\dfrac{1}{k}\sum\limits_{i=0}^{k-1}x^{(a-b)i}$
若 $x$ 的阶恰为 $k$ ，则 $x^{(a-b)}$ 必不为 $1$。在经典的范德蒙德矩阵证明中，由此可得上式 $=0$
但是，此时由于无法做除法，等比数列求和公式不再生效，无法证明上式 $=0$。
但是，仍然恒有 $0=(x^k-1)=(x-1)(1+x+x^2+...+x^k)$ 存在，也就是说，上面求和的结果虽然可能非 $0$ ，但必然是零因子。
也就是说，我们算得的结果是由 (真实答案)+(零因子的线性组合) 构成的。
接下来要找到一种合适的扩域方法，这样就能避免零因子问题。
不妨取分圆多项式 $\Phi_k(x)$ ，下面不加证明地给出两个定理 ：
（相关知识可见《近世代数乱编》）
- ① 分圆多项式在 $Q$ 上不可约。
这保证了 $\pmod {\Phi_k(x)}$ 意义下的“数”是个域，不存在零因子。
这里注意，还要验证该多项式在 $F_p$ 下是否能够分解，一般情况下是不能的，此时可以正常使用。
若计算时没有利用 $F_p$ 的性质（如求逆元），则可以看作大整数计算，最后将答案取模，此时不需要检查 $F_p$ 下是否能够分解。
- ② 在 $\pmod {\Phi_k(x)}$ 意义下， $x$ 的阶恰好为 $k$。
这又保证了我们构造的前提成立。
这样，只需在 $\pmod {\Phi_k(x)}$ 下计算，即可得到满足题意的答案。
但是，模多项式的计算较为繁琐，常数较大且不易优化。
而更好的消息是，由于 $\Phi_k(x)|(x^k-1)$ ，则有 $F(x)\bmod \Phi_k(x)=F(x)\bmod (x^k-1)\bmod \Phi_k(x)$
这表明我们可以先用“循环卷积多项式”代替严格的扩域，到最后再对 $\Phi_k(x)$ 取模。
此时任意两个数的乘法就变成 $O(k^2)$ 多项式循环卷积，复杂度升高了 $O(k^2)$。
观察矩阵可得，每次乘的都是一个单项式，复杂度就只需升高 $O(k)$。
暴力做的话,总的复杂度是 $O(k^{n+2}n)$ 。
可以用 $\rm FFT$ 优化到 $O(k^{n+1}n\log k)$ ，不过大多数时候不实用。
- **例题①** : [CF1103E Radix sum](https://www.luogu.org/problem/CF1103E)
$k$ 进制不进位加法卷积快速幂。
- **例题②** : [P5577 [CmdOI2019]算力训练](https://www.luogu.org/problem/P5577)
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
------------
# 0x04.子集卷积 & 集合幂级数
- [P6097 【模板】子集卷积](https://www.luogu.com.cn/problem/P6097)
作用 : 求 $S[s]=\sum\limits_{i\subseteq s}F[i]G[s\ {\rm xor}\ i]$
另一种说法 : $S[s]=\sum\limits_{\tiny\begin{matrix}i|j=s\\i\&j=0\end{matrix}}F[i]G[j]$
组合意义就是每次挑选不交的两个集合拼成新的集合。
好像在某些 $\rm DP$ 题目里面见过? 这确实可以用来优化某些毒瘤子集 $\rm DP$。
使用枚举子集的技巧即可做到 $O(3^n)$ 计算。考虑使用卷积知识加速计算。
直接套用前面学习的位运算卷积只能求 $S[s]=\sum\limits_{i|j=s}F[i]G[j]$ ,无法再满足 $i\&j=0$ 。
考虑到 $i|j=s,\ i\&j=0$ 等价于 $i|j=s,|i|+|j|=|s|$。
( 其中 $|i|$ 表示 $i$ 在二进制下的 $1$ 个数 )
我们可以将原来的 $F[k]$ 替换成占位多项式 $F[k]x^{|k|}$。
这样，计算 $\rm or$ 卷积之后，利用 $x$ 上的加法卷积， $[x^p]S[k]$ 就是满足 $i|j=s,|i|+|j|=p$ 的 $(i,j)$ 的贡献。
我们取出 $[x^{|k|}]S[k]$ 即为我们想要的答案。
考虑计算的复杂度，多项式加法和数乘的复杂度是 $O(n)$ 的，故 $\rm FWT$ 变换的复杂度变为 $O(2^nn^2)$。
中间还需要点乘，多项式乘法用暴力实现，复杂度也是 $O(2^nn^2)$。
在具体实现中，可以把 $x$ 维度放在前面以减小常数。（代码中并没有这样做）
```cpp
#include<algorithm>
#include<cstdio>
#define MaxN 1050000
using namespace std;
const int mod=1000000009;
int read(){
int X=0;char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
int m;
struct Poly
{
int a[21];
void operator += (const Poly &B){
for (int i=0;i<=m;i++)
a[i]=(a[i]+B.a[i])%mod;
}
void operator -= (const Poly &B){
for (int i=0;i<=m;i++)
a[i]=(a[i]-B.a[i])%mod;
}
Poly operator * (const Poly &B) const{
Poly R;
for (int i=0;i<=m;i++)R.a[i]=0;
for (int i=0;i<=m;i++)
for (int j=0;i+j<=m;j++)
R.a[i+j]=(R.a[i+j]+1ll*a[i]*B.a[j])%mod;
return R;
}
};
void DWT(Poly *f,int n)
{
for (int l=1;l<n;l<<=1)
for (int p=0;p<n;p+=l+l)
for (int k=0;k<l;k++)
f[p|l|k]+=f[p|k];
}
void IDWT(Poly *f,int n)
{
for (int l=1;l<n;l<<=1)
for (int p=0;p<n;p+=l+l)
for (int k=0;k<l;k++)
f[p|l|k]-=f[p|k];
}
Poly F[MaxN],G[MaxN],T[MaxN];
int n,c[MaxN];
int main()
{
m=read();n=(1<<m);
for (int i=1;i<n;i++)c[i]=c[i>>1]+(i&1);
for (int i=0;i<n;i++)F[i].a[c[i]]=read();
for (int i=0;i<n;i++)G[i].a[c[i]]=read();
DWT(F,n);DWT(G,n);
for (int i=0;i<n;i++)
F[i]=F[i]*G[i];
IDWT(F,n);
for (int i=0;i<n;i++)
printf("%d ",(F[i].a[c[i]]+mod)%mod);
return 0;
}
```
**例题** ：[CF914G Sum the Fibonacci](https://www.luogu.org/problem/CF914G)
这就是一道很好的模板题 (不过值域只有 $2^{17}$ 所以很多人 $3^{17}$ 艹过去了)。
令 $cnt[i]=\sum\limits_{j}[s[j]=i]$
求出 $cnt$ 和本身的子集卷积 $F1$ ， $cnt$ 和本身的异或卷积 $F2$。
令 $G1[i]=fib(F1[i])$ ， $G2$ 类似， $cnt'[i]=fib(cnt[i])$
最后把 $G1,G2,cnt'$ 卷在一起就好。
- **半在线子集卷积**
给出集合幂级数 $W,C$ ，求幂级数 $S$ 使得 ：
$$S[s]=C[s]\sum\limits_{\small t\subsetneq s}S[t]W[s-t]$$
对于 $S[s]$ ，需要得知所有为 $s$ 的真子集的位置的 $S$ 才能计算。
我们可以按照 $|s|$ 的顺序计算 $S[s]$ ，这样就能保证需要的值已经被计算好了。
此时为了方便需要将 $x$ 维度放在前面。具体实现见例题。
**例题** : [P4221 [WC2018]州区划分](https://www.luogu.org/problem/P4221) | [Solution](https://www.luogu.com.cn/blog/command-block/post-shuo-xue-ji-lu-p4221-wc2018-zhou-ou-hua-fen)
- **集合幂级数运算进阶**
需熟知如何在 $O(n^2)$ 内做一系列多项式操作，可见 [NTT与多项式全家桶](https://www.luogu.com.cn/blog/command-block/ntt-yu-duo-xiang-shi-quan-jia-tong)。
下面以 $F(x)=\sum\limits_{s}F[s]x^s$ 来描述一个集合幂级数。
注意，其中 $x$ 的上标是一个集合。
- **求逆**
变换后将占位多项式求逆，然后逆变换回来即可。
集合幂级数卷积的单位元是 $x^{\varnothing}$ ，可以将其视作常数项 $1$。
- $\bf exp,ln$
设有集合幂级数 $F$ ，模仿多项式级数来定义其 $\exp,\ln$。
定义 $\exp F=\sum\limits_{k=0}\dfrac{F^k}{k!}$ ,要求 $[x^{\varnothing}]F=0$ 。
定义 $\ln F$ 为 $\exp$ 的逆运算，要求 $[x^{\varnothing}]F=1$ 。
也有定义式 $\ln F=\sum\limits_{k=1}\dfrac{(-1)^{k+1}(F-x^{\varnothing})^k}{k}$
也只需要变换后对占位多项式进行 $\exp,\ln$ 即可计算。
集合幂级数的运算也有和生成函数类似的组合意义。
- **例题①** ： [P6570 [NOI Online #3 提高组]优秀子序列](https://www.luogu.com.cn/problem/P6570)
不难发现，对于序列中的 $A[i]$ 写出集合幂级数 $(x^{\varnothing}+x^{A[i]})$ ，我们的目标是出所有集合幂级数的卷积。
考虑模仿快速 $01$ 背包的套路，对其先取 $\ln$ 后 $\exp$。
$P(x)=\prod\limits_{i=1}^n(x^{\varnothing}+x^{A[i]})$
$=\exp\sum\limits_{i=1}^n\ln(x^{\varnothing}+x^{A[i]})$
$=\exp\sum\limits_{i=1}^n\sum\limits_{k=1}\dfrac{(-1)^{k+1}(x^{A[i]})^k}{k}$
而 $(x^{A[i]})^k$ 显然只有 $k=1$ 时为 $x^{A[i]}$ ，其他情况为 $0$。（注意这是子集卷积）
$=\exp\sum\limits_{i=1}^nx^{A[i]}$
似乎变成 $\exp$ 板子了呢……
不过要注意对 $A[i]=0$ 的情况特殊处理。
- **例题②** ： [Loj#154. 集合划分计数](https://loj.ac/p/154)
不难发现“划分”就对应生成无序集合，这即为 $\exp$ 的组合意义。
本题要求划分的集合数 $\leq k$ ，则是一个部分 $\exp$ ：
$$\sum\limits_{i=0}^k\dfrac{F(x)^k}{k!}$$
对占位多项式做部分 $\exp$ 即可。
有 $\rm ODE$ : $G=\sum\limits_{i=0}^k\dfrac{F^k}{k!}\Rightarrow G'=F'*\big(G-(F^k/k!)\big)$
常数莫名其妙的大，跑不过 $O(2^nn^3)$ 实锤……
[评测记录](https://loj.ac/s/1001124)
- **例题③** ： [Uoj#94. 【集训队互测2015】胡策的统计](https://uoj.ac/problem/94)
设 $G[s]$ 为点集 $s$ 的联通生成子图个数，$F[s]$ 为点集 $s$ 的生成子图个数。
$F$ 是易求的，有 $F[s]=2^{\text{s内部边数}}$。
根据 $\exp$ 的组合意义显然有 $F=\exp G$ ，则 $G=\ln F$。
设 $P[s]$ 为点集 $s$ 的生成子图的连通值之和。
枚举连通块数目则有 $P=\sum\limits_{k=0}\dfrac{G^k}{k!}k!=\dfrac{1}{1-G}$。
求逆即可。
[评测记录](https://uoj.ac/submission/441128)
- **例题④** ： [Loj#3165. 「CEOI2019」游乐园](https://loj.ac/p/3165)
能够发现无环图所有边反向之后还是无环图，互反的一对图翻转边数和恰为 $m$。所以答案为方案数乘以 $m/2$。
设 $F[s]$ 表示点集 $s$ 内适当反转后无环的方案数。
考虑拆 $\rm DAG$ 为子问题，每次去除入度为 $0$ 的所有点。
钦定入度为 $0$ 的点集 $t\subseteq s$ ，$t$ 内部必须没有边，即 $t$ 为独立集。
$t$ 和其余点 $s-t$ 之间的边方向确定。所以方案数为 $[t\text{是独立集}]F[s-t]$
由于我们只能钦定一些点入度为 $0$ ，所以还需容斥。
可得 $F[s]=\sum\limits_{t\subseteq s}(-1)^{|t|-1}[t\text{是独立集}]F[s-t]$。
设 $G=(-1)^{|t|-1}[t\text{是独立集}]$
即 $F=F*G+1$。可得 $F=\dfrac{1}{1-G}$，求逆即可。
[评测记录](https://loj.ac/s/1001418)
- **例题⑤** ： [#6673. EntropyIncreaser 与山林](https://loj.ac/p/6673)
可能要先去 [多项式计数杂谈](https://www.luogu.com.cn/blog/command-block/sheng-cheng-han-shuo-za-tan) 看看怎么数欧拉图。
首先要求出偶度图的集合幂级数，然后求 $\ln$ 即可得到欧拉图。
现在问题变为对每个点集求生成偶度图个数。
找出任意一个生成森林，将其余的边随意选取，此时能得到目前每个点的奇偶性。
由于森林无环的性质，每种电度奇偶性都能通过恰当地选取森林中的边得到，且方案数恰好为 $1$。
因此，偶度图个数为 $2^{\text{边数}-\text{生成森林边数}}$。
[评测记录](https://loj.ac/s/1032156)