专栏文章

位运算卷积(FWT) & 集合幂级数

算法·理论参与者 47已保存评论 55

文章操作

快速查看文章及其快照的属性,并进行相关操作。

当前评论
55 条
当前快照
1 份
快照标识符
@mhz5u0c5
此快照首次捕获于
2025/11/15 01:57
3 个月前
此快照最后确认于
2025/11/29 05:25
3 个月前
查看原文
2020.10:\rm2020.10: 进行了小修小补。
2020.11:\rm2020.11: 重制了更严谨的 “k进制不进位加法卷积” 部分。
2020.11:\rm2020.11: 加更集合幂级数。

0x01.FWT\rm FWT 概论

  • 位运算卷积
众所周知,多项式乘法是加法卷积,因为第 ii 项和第 jj 项的乘积贡献到第 i+ji+j 项。
类似地定义位运算卷积 : 第 ii 项和第 jj 项的乘积贡献到第 iji⊕j 项。其中 是某种位运算
S[k]=ij=kA[i]B[j]S[k]=\sum\limits_{i⊕j=k}A[i]B[j] ,记作卷积式 AB=SA*B=S
  • 构造的尝试
众所周知,FFT\rm FFT 把多项式转换成点值之后,从卷积变为了直接点积。
我们自然也期望把位运算卷积转化成点积。
FWT(A)FWT(A) 是幂级数 AA 经过 FWT\rm FWT 变换之后得到的幂级数。
我们需要令其满足 : AB=CFWT(A)FWT(B)=FWT(C)A*B=C \Longleftrightarrow FWT(A)·FWT(B)=FWT(C) (点积)。
FFT\rm FFT 是一个线性变换,我们也希望 FWT\rm FWT 变换是线性的。
我们还不知道怎么变换,于是设 c(i,j)c(i,j) 为变换系数,即 A[j]A[j]FWT(A)[i]FWT(A)[i] 的贡献系数。
FWT(A)[i]=j=0n1c(i,j)AjFWT(A)[i]=\sum\limits_{j=0}^{n-1}c(i,j)A_j
根据 FWT(A)FWT(B)=FWT(C)FWT(A)·FWT(B)=FWT(C) ,得到 :
FWT(A)[i]FWT(B)[i]=FWT(C)[i]FWT(A)[i]FWT(B)[i]=FWT(C)[i]
j=0n1c(i,j)A[j]k=0n1c(i,k)B[k]=p=0n1c(i,p)C[p]\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]
j=0n1k=0n1c(i,j)c(i,k)A[j]B[k]=p=0n1c(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]
根据 AB=CA*B=C ,又能得到 :
C[p]=t1t2=pA[t1]B[t2]C[p]=\sum\limits_{t_1⊕t_2=p}A[t_1]B[t_2]
p=0n1c(i,p)C[p]=p=0n1c(i,p)t1t2=pA[t1]B[t2]\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]
j=0n1k=0n1c(i,j)c(i,k)A[j]B[k]=p=0n1c(i,p)t1t2=pA[t1]B[t2]\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]
j=0n1k=0n1c(i,j)c(i,k)A[j]B[k]=p=0n1t1t2=pA[t1]B[t2]c(i,t1t2)=t1=0n1t2=0n1A[t1]B[t2]c(i,t1t2)\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,jk)c(i,j)c(i,k)=c(i,j⊕k) ,只需要满足这个就好了。
另外,由于位运算每一位的独立性, c(i,j)c(i,j) 也有一个重要性质: 可以分位考虑。
假设我们已经(根据真值表)求出满足要求的 c([0,1],[0,1])c([0,1],[0,1]),我们能这样构造出所有的 cc :
设二进制数 aa 的每一位分别为 : a0,a1,a2...a_0,a_1,a_2...
则有 c(i,j)=c(i0,j0)c(i1,j1)c(i2,j2)...c(i,j)=c(i_0,j_0)c(i_1,j_1)c(i_2,j_2)... ,就是把每一位的变换系数乘起来
那么 : 对于每个 tt ,都有 c(it,jt)c(it,kt)=c(it,jtkt)c(i,j)c(i,k)=c(i,jk)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)
这就符合我们的条件。
好了,现在假设我们已经有了符合要求的 cc ,如何快速求解 FWT\rm FWT 变换呢?
FWT(A)[i]=j=0n1c(i,j)A[j]FWT(A)[i]=\sum\limits_{j=0}^{n-1}c(i,j)A[j]
根据这个式子直接求和,复杂度至少是 O(n2)O(n^2) 的,并没有起到优化作用。
我们考虑按位拆半
=j=0(n/2)1c(i,j)A[j]+j=(n/2)n1c(i,j)A[j]=\sum\limits_{j=0}^{(n/2)-1}c(i,j)A[j]+\sum\limits_{j=(n/2)}^{n-1}c(i,j)A[j]
ii'ii 去掉二进制首位剩下的数。
在首位分开考虑 :
=j=0(n/2)1c(i0,j0)c(i,j)A[j]+j=(n/2)n1c(i0,j0)c(i,j)A[j]=\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(i0,0)j=0(n/2)1c(i,j)A[j]+c(i0,1)j=(n/2)n1c(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)c(i',j') 就是去掉最高位的子变换,这里规模减半了。
A0A_0 为幂级数下标首位为 00 的部分,类似地有 A1A_1
i0=0i_0=0 ,则有 :
FWT(A)[i]=c(0,0)FWT(A0)[i]+c(0,1)FWT(A1)[i](i[0,n/2))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)
i0=1i_0=1 ,则有 :
FWT(A)[i+(n/2)]=c(1,0)FWT(A0)[i]+c(1,1)FWT(A1)[i](i[0,n/2))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)O(n) 的代价,根据上列式子合并两个规模为 n/2n/2 的子变换。
所以,若 n=2mn=2^m ,需要合并 mm 次,复杂度为 O(m2m)O(m2^m)
( 可能有点抽象,但是您如果写过FFT,看到代码就会懂了 )
此外,逆变换 IFWT\rm IFWT 就是对 cc 矩阵求个逆,具体见下文。
( 一个重要的地方是,这个构造出来的 cc 矩阵一定要有逆,否则就变换不回去TAT )

0x02.基础位运算卷积

针对不同的位运算,根据 c(i,j)c(i,k)=c(i,jk)c(i,j)c(i,k)=c(i,j⊕k) 构造出 c([0,1],[0,1])c([0,1],[0,1]) 即可。
我们把这个矩阵称为位矩阵
构造的过程可能有些繁杂,可以直接记结论,或者去后面看扩展版的。

1.1 Or1.1\ \rm Or 卷积

设位矩阵为 c=[c(0,0)c(0,1)c(1,0)c(1,1)]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,jk)c(i,j)c(i,k)=c(i,j|k)
  • c(0,0)c(0,0)=c(0,00)c(0,0)c(0,0)=c(0,0|0)
    c(0,0)=1\Rightarrow c(0,0)=100
同理不难推知 c(_,_){0,1}c(\_,\_)∈\{0,1\}
  • c(0,1)c(0,0)=c(0,10)c(0,1)c(0,0)=c(0,1|0)
    c(0,1)=0\Rightarrow c(0,1)=0c(0,0)=c(0,1)=1c(0,0)=c(0,1)=1
  • c(1,1)c(1,0)=c(1,10)c(1,1)c(1,0)=c(1,1|0)
    c(1,1)=0\Rightarrow c(1,1)=0c(1,0)=c(1,1)=1c(1,0)=c(1,1)=1
首先,如果有一排0或者一列0那么这个矩阵就没有逆,那么可以构造出:
两种情况 : [1110]\begin{bmatrix}1&1\\1&0\end{bmatrix}[1011]\begin{bmatrix}1&0\\1&1\end{bmatrix}
Tips : OrOr 卷积的上面第二个矩阵 FWTFWT 相当于子集求和。
原因:第二个矩阵相当于 c(i,j)=[i&j=j]c(i,j)=[i\&j=j]
Ai=i&j=jAiA'_i=\sum\limits_{i\&j=j}A_i等价于Ai=jiAiA'_i=\sum\limits_{j∈i}A_i
(也可以使用高维前缀和来推导)
(下面采用第二个矩阵)
FWT(A)[i]=FWT(A0)[i]FWT(A)[i]=FWT(A_0)[i]
FWT(A)[i+(n/2)]=FWT(A0)[i]+FWT(A1)[i]FWT(A)[i+(n/2)]=FWT(A_0)[i]+FWT(A_1)[i]
对于逆变换,把矩阵求个逆可得 [1011]\begin{bmatrix}1&0\\-1&1\end{bmatrix}
IFWT(A)[i]=IFWT(A0)[i]IFWT(A)[i]=IFWT(A_0)[i]
IFWT(A)[i+(n/2)]=IFWT(A1)[i]IFWT(A0)[i]IFWT(A)[i+(n/2)]=IFWT(A_1)[i]-IFWT(A_0)[i]

1.2 And1.2\ \rm And 卷积

起点 : c(i,j)c(i,k)=c(i,j&k)c(i,j)c(i,k)=c(i,j\&k)
同上,容易得到 c(_,_){0,1}c(\_,\_)∈\{0,1\}
  • c(0,1)c(0,0)=c(0,1&0)c(0,1)c(0,0)=c(0,1\&0)
    c(0,0)=0\Rightarrow c(0,0)=0c(0,0)=c(0,1)=1c(0,0)=c(0,1)=1
  • c(1,1)c(1,0)=c(1,1&0)c(1,1)c(1,0)=c(1,1\&0)
    c(1,0)=0\Rightarrow c(1,0)=0c(1,0)=c(1,1)=1c(1,0)=c(1,1)=1
还是老样子,如果有一排 00 或者一列 00 那么这个矩阵就没有逆,那么可以构造出:
两种情况 : [0111]\begin{bmatrix}0&1\\1&1\end{bmatrix}[1101]\begin{bmatrix}1&1\\0&1\end{bmatrix},下面采用第二种。
FWT(A)[i]=FWT(A0)[i]+FWT(A1)[i]FWT(A)[i]=FWT(A_0)[i]+FWT(A_1)[i]
FWT(A)[i+(n/2)]=FWT(A1)[i]FWT(A)[i+(n/2)]=FWT(A_1)[i]
把矩阵求个逆可得[1101]\begin{bmatrix}1&-1\\0&1\end{bmatrix}:
IFWT(A)[i]=IFWT(A0)[i]IFWT(A1)[i]IFWT(A)[i]=IFWT(A_0)[i]-IFWT(A_1)[i]
IFWT(A)[i+(n/2)]=IFWT(A1)[i]IFWT(A)[i+(n/2)]=IFWT(A_1)[i]

1.3 Xor1.3\ \rm Xor 卷积

起点 : c(i,j)c(i,k)=c(i,j xor k)c(i,j)c(i,k)=c(i,j\ \text{xor}\ k)
  • 对于任意的 x,yx,y ,均有 c(0,0)c(x,y)=c(x,y xor 0)=c(x,y)c(0,0)c(x,y)=c(x,y\ \text{xor}\ 0)=c(x,y)
    c(0,0)=1\Rightarrow c(0,0)=1.
  • c(1,1)c(1,1)=c(1,0)c(1,1)c(1,1)=c(1,0)
    此时若 c(1,1)=c(1,0)=0c(1,1)=c(1,0)=0,则一行为 00 ,矩阵无逆。
    所以 c(1,1),c(1,0)c(1,1),c(1,0) 必然都非 00
  • c(1,0)c(1,1)=c(1,1)c(1,0)c(1,1)=c(1,1)
    刚才说c(1,1)c(1,1)00,所以此处 c(1,0)c(1,0) 一定是1.
  • c(0,1)c(0,1)=c(0,0)c(0,1)c(0,1)=c(0,0)
    c(0,1){1,1}\Rightarrow c(0,1)∈\{-1,1\}
两种情况 : [1111]\begin{bmatrix}1&1\\-1&1\end{bmatrix}[1111]\begin{bmatrix}1&1\\1&-1\end{bmatrix} ,下面采用第二种。
附 :不难观察出 c(i,j)=(1)i&jc(i,j)=(-1)^{|i\&j|}
FWT(A)i=FWT(A0)i+FWT(A1)iFWT(A)_i=FWT(A_0)_i+FWT(A_1)_i
FWT(A)i+(n/2)=FWT(A0)iFWT(A1)iFWT(A)_{i+(n/2)}=FWT(A_0)_i-FWT(A_1)_{i}
求逆可得[0.50.50.50.5]\begin{bmatrix}0.5&0.5\\0.5&-0.5\end{bmatrix}
IFWT(A)i=IFWT(A0)i+IFWT(A1)i2IFWT(A)_i=\dfrac{IFWT(A_0)_i+IFWT(A_1)_i}{2}
IFWT(A)i+(n/2)=IFWT(A0)iIFWT(A1)i2IFWT(A)_{i+(n/2)}=\dfrac{IFWT(A_0)_i-IFWT(A_1)_i}{2}

1.41.4 模板题 & Code:

  • FWT :
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.FWT\rm FWT 变换的一些性质

前面讲了, FWT\rm FWT 变换的本质是线性变换。
这意味着 FWT(A+B)=FWT(A)+FWT(B)FWT(A+B)=FWT(A)+FWT(B) ,且 FWT(cA)=cFWT(A)FWT(cA)=cFWT(A)
此外,如果需要卷积的向量只有少数非 00 项( 2,32,3 个)可能会有分类讨论的神奇的解法。
例题(我只能做到这些了):
题做多了再来填坑吧。

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

0x04.位值域的扩展

其实位运算的本质是对一个 nn0101 向量的运算。
或运算就是每一维取 max\max。且运算就是每一维取 minmin
异或运算则是每一维对应相加再 mod 2\bmod\ 2
位运算有个特点 : 向量的每一位都是独立的。
我们把{0,1}\{0,1\}扩展到 [0,k)Z[0,k)∩Z 也就是扩展到 kk 进制,看看会得到什么?

每一维取 max\bf max

对应原来的按位 or\rm or
可得 c(x,y)c(x,y)=c(x,max(y,y))=c(x,y)c(x,y)c(x,y)=c(x,max(y,y))=c(x,y)
所以整个矩阵中只有 0011
又可得 c(x,y1)c(x,y2)=c(x,max(y1,y2))c(x,y1)c(x,y2)=c(x,max(y1,y2))
假如有 c(x,y1)=1, c(x,y2)=0c(x,y1)=1,\ c(x,y2)=0,可得 c(x,y1)c(x,y2)=c(x,max(y1,y2))c(x,y1)c(x,y2)=c(x,max(y1,y2))
又因为 10=01*0=0 所以 max(y1,y2)=y2max(y1,y2)=y2 ,可得 y1<y2y1<y2.
也就是说,每一行以内, 11 总是在 00 的前面。
接下来,除了矩阵有逆之外,再没有别的要求了。
如果要有逆的话,每一行都必须不同,那么 11 的个数分别就是 1...n1...n ,随意排列都可以。
例子 : k=4k=4 的情形
[1000110011101111]\begin{bmatrix}1&0&0&0\\1&1&0&0\\1&1&1&0\\1&1&1&1\end{bmatrix}或者[1111100011101100]\begin{bmatrix}1&1&1&1\\1&0&0&0\\1&1&1&0\\1&1&0&0\end{bmatrix}等等4!4!种。
为了方便,一般取用第一种,求逆可得[1000110001100011]\begin{bmatrix}1&0&0&0\\-1&1&0&0\\0&-1&1&0\\0&0&-1&1\end{bmatrix}
暴力做的话 O(kn+1n)O(k^{n+1}n)
这个矩阵本质上就是前缀和和差分……所以可以优化到 O(knn)O(k^nn)
(宏观上就是高维前缀和)

每一维取 min\bf min

对应原来的按位 and\rm and
类似的,整个矩阵中只有 0011 ,每一行以内, 11 总是在 00 的后面。
例子 : k=4k=4 的情形
[1111011100110001]\begin{bmatrix}1&1&1&1\\0&1&1&1\\0&0&1&1\\0&0&0&1\end{bmatrix}或者[1111001101110001]\begin{bmatrix}1&1&1&1\\0&0&1&1\\0&1&1&1\\0&0&0&1\end{bmatrix}等等4!4!种。
为了方便,一般取用第一种,求逆可得[1100011000110001]\begin{bmatrix}1&-1&0&0\\0&1&-1&0\\0&0&1&-1\\0&0&0&1\end{bmatrix}
这个矩阵本质上就是后缀和和差分……同样能优化到 O(knn)O(k^nn)
(宏观上就是高维后缀和)

每一维加起来 mod k\bf mod\ k

对应原来的按位 xor\rm xor
这玩意就比较复杂了,需要动用单位根……
我们知道 c(x,y1)c(x,y2)=c(x,(y1+y2)mod k)c(x,y1)c(x,y2)=c(x,(y1+y2)mod\ k)
c(x,y)=ωkyc(x,y)=\omega_k^y 就能满足要求: c(x,y1)c(x,y2)=ωky1ωky2=ωk(y1+y2)mod kc(x,y1)c(x,y2)=\omega_k^{y1}\omega_k^{y2}=\omega_k^{(y1+y2)mod\ k}
但是,每一行都一样的话,矩阵就没有逆。
你会发现直接把 FFT 的范德蒙德矩阵拿来用就好了:
[111...11wk1wk2...wkk11wk2wk4...wk2(k1)...............1wkk1wk2(k1)...wk(k1)(k1)]\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}
逆矩阵我们也知道,那就是:
1k[111...11wk1wk2...wk(k1)1wk2wk4...wk2(k1)...............1wk(k1)wk2(k1)...wk(k1)(k1)]\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(kn+1n)O(k^{n+1}n) 的。
可以使用 FFT\rm FFT 优化到 O(knnlogk)O(k^nn\log k)
但是,单位根在模意义很可能不存在。考虑扩充我们“数”的表示。
首先想到的是,人为地定义代数对象 xx 满足 xk=1x^k=1 (且 xx 的阶恰为 kk),用其代替单位根 wkw_k,不难验证其满足前文构造矩阵的条件。
此举相当于用 (modxk1)\pmod{x^{k}-1} 下的循环卷积多项式替换了“数”。
我们的位矩阵长这样 :
C1=[111...11x1x2...xk11x2x4...x2(k1)...............1xk1x2(k1)...x(k1)(k1)]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}
C2=1k[111...11x1x2...x(k1)1x2x4...x2(k1)...............1x(k1)x2(k1)...x(k1)(k1)]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}
接下来的推导可能需要抽象代数和高等代数相关知识。详情可见 : 近世代数乱编
问题在于,此时 (modxk1)\pmod{x^{k}-1} 下的循环卷积多项式并不是域,其存在零因子。
这样就会导致两个位矩阵并非严格互逆。
具体而言,我们原本希望 C1C2=IC_1*C_2=I ,这需要满足 i=0k1C1[a][i]C2[i][b]=[a=b]\sum\limits_{i=0}^{k-1}C_1[a][i]C_2[i][b]=[a=b]
a=ba=b 时,有 i=0n1C1[a][i]C2[i][a]=1ki=0k1xaixai=1\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 ,这不需要消去律也成立,这部分贡献是没有差错的。
aba≠b 时,有 i=0n1C1[a][i]C2[i][b]=1ki=0k1xaixbi=1ki=0k1x(ab)i\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}
xx 的阶恰为 kk ,则 x(ab)x^{(a-b)} 必不为 11。在经典的范德蒙德矩阵证明中,由此可得上式 =0=0
但是,此时由于无法做除法,等比数列求和公式不再生效,无法证明上式 =0=0
但是,仍然恒有 0=(xk1)=(x1)(1+x+x2+...+xk)0=(x^k-1)=(x-1)(1+x+x^2+...+x^k) 存在,也就是说,上面求和的结果虽然可能非 00 ,但必然是零因子。
也就是说,我们算得的结果是由 (真实答案)+(零因子的线性组合) 构成的。
接下来要找到一种合适的扩域方法,这样就能避免零因子问题。
不妨取分圆多项式 Φk(x)\Phi_k(x) ,下面不加证明地给出两个定理 :
(相关知识可见《近世代数乱编》)
  • ① 分圆多项式在 QQ 上不可约。
    这保证了 (modΦk(x))\pmod {\Phi_k(x)} 意义下的“数”是个域,不存在零因子。
    这里注意,还要验证该多项式在 FpF_p 下是否能够分解,一般情况下是不能的,此时可以正常使用。
    若计算时没有利用 FpF_p 的性质(如求逆元),则可以看作大整数计算,最后将答案取模,此时不需要检查 FpF_p 下是否能够分解。
  • ② 在 (modΦk(x))\pmod {\Phi_k(x)} 意义下, xx 的阶恰好为 kk
    这又保证了我们构造的前提成立。
这样,只需在 (modΦk(x))\pmod {\Phi_k(x)} 下计算,即可得到满足题意的答案。
但是,模多项式的计算较为繁琐,常数较大且不易优化。
而更好的消息是,由于 Φk(x)(xk1)\Phi_k(x)|(x^k-1) ,则有 F(x)modΦk(x)=F(x)mod(xk1)modΦk(x)F(x)\bmod \Phi_k(x)=F(x)\bmod (x^k-1)\bmod \Phi_k(x)
这表明我们可以先用“循环卷积多项式”代替严格的扩域,到最后再对 Φk(x)\Phi_k(x) 取模。
此时任意两个数的乘法就变成 O(k2)O(k^2) 多项式循环卷积,复杂度升高了 O(k2)O(k^2)
观察矩阵可得,每次乘的都是一个单项式,复杂度就只需升高 O(k)O(k)
暴力做的话,总的复杂度是 O(kn+2n)O(k^{n+2}n)
可以用 FFT\rm FFT 优化到 O(kn+1nlogk)O(k^{n+1}n\log k) ,不过大多数时候不实用。

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

0x04.子集卷积 & 集合幂级数

作用 : 求 S[s]=isF[i]G[s xor i]S[s]=\sum\limits_{i\subseteq s}F[i]G[s\ {\rm xor}\ i]
另一种说法 : S[s]=ij=si&j=0F[i]G[j]S[s]=\sum\limits_{\tiny\begin{matrix}i|j=s\\i\&j=0\end{matrix}}F[i]G[j]
组合意义就是每次挑选不交的两个集合拼成新的集合。
好像在某些 DP\rm DP 题目里面见过? 这确实可以用来优化某些毒瘤子集 DP\rm DP
使用枚举子集的技巧即可做到 O(3n)O(3^n) 计算。考虑使用卷积知识加速计算。
直接套用前面学习的位运算卷积只能求 S[s]=ij=sF[i]G[j]S[s]=\sum\limits_{i|j=s}F[i]G[j] ,无法再满足 i&j=0i\&j=0
考虑到 ij=s, i&j=0i|j=s,\ i\&j=0 等价于 ij=s,i+j=si|j=s,|i|+|j|=|s|
( 其中 i|i| 表示 ii 在二进制下的 11 个数 )
我们可以将原来的 F[k]F[k] 替换成占位多项式 F[k]xkF[k]x^{|k|}
这样,计算 or\rm or 卷积之后,利用 xx 上的加法卷积, [xp]S[k][x^p]S[k] 就是满足 ij=s,i+j=pi|j=s,|i|+|j|=p(i,j)(i,j) 的贡献。
我们取出 [xk]S[k][x^{|k|}]S[k] 即为我们想要的答案。
考虑计算的复杂度,多项式加法和数乘的复杂度是 O(n)O(n) 的,故 FWT\rm FWT 变换的复杂度变为 O(2nn2)O(2^nn^2)
中间还需要点乘,多项式乘法用暴力实现,复杂度也是 O(2nn2)O(2^nn^2)
在具体实现中,可以把 xx 维度放在前面以减小常数。(代码中并没有这样做)
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;
}
这就是一道很好的模板题 (不过值域只有 2172^{17} 所以很多人 3173^{17} 艹过去了)。
cnt[i]=j[s[j]=i]cnt[i]=\sum\limits_{j}[s[j]=i]
求出 cntcnt 和本身的子集卷积 F1F1cntcnt 和本身的异或卷积 F2F2
G1[i]=fib(F1[i])G1[i]=fib(F1[i])G2G2 类似, cnt[i]=fib(cnt[i])cnt'[i]=fib(cnt[i])
最后把 G1,G2,cntG1,G2,cnt' 卷在一起就好。
  • 半在线子集卷积
    给出集合幂级数 W,CW,C ,求幂级数 SS 使得 :
    S[s]=C[s]tsS[t]W[st]S[s]=C[s]\sum\limits_{\small t\subsetneq s}S[t]W[s-t]
对于 S[s]S[s] ,需要得知所有为 ss 的真子集的位置的 SS 才能计算。
我们可以按照 s|s| 的顺序计算 S[s]S[s] ,这样就能保证需要的值已经被计算好了。
此时为了方便需要将 xx 维度放在前面。具体实现见例题。
  • 集合幂级数运算进阶
    需熟知如何在 O(n2)O(n^2) 内做一系列多项式操作,可见 NTT与多项式全家桶
    下面以 F(x)=sF[s]xsF(x)=\sum\limits_{s}F[s]x^s 来描述一个集合幂级数。
    注意,其中 xx 的上标是一个集合。
    • 求逆
      变换后将占位多项式求逆,然后逆变换回来即可。
      集合幂级数卷积的单位元是 xx^{\varnothing} ,可以将其视作常数项 11
    • exp,ln\bf exp,ln
      设有集合幂级数 FF ,模仿多项式级数来定义其 exp,ln\exp,\ln
      定义 expF=k=0Fkk!\exp F=\sum\limits_{k=0}\dfrac{F^k}{k!} ,要求 [x]F=0[x^{\varnothing}]F=0
      定义 lnF\ln Fexp\exp 的逆运算,要求 [x]F=1[x^{\varnothing}]F=1
      也有定义式 lnF=k=1(1)k+1(Fx)kk\ln F=\sum\limits_{k=1}\dfrac{(-1)^{k+1}(F-x^{\varnothing})^k}{k}
      也只需要变换后对占位多项式进行 exp,ln\exp,\ln 即可计算。
      集合幂级数的运算也有和生成函数类似的组合意义。
不难发现,对于序列中的 A[i]A[i] 写出集合幂级数 (x+xA[i])(x^{\varnothing}+x^{A[i]}) ,我们的目标是出所有集合幂级数的卷积。
考虑模仿快速 0101 背包的套路,对其先取 ln\lnexp\exp
P(x)=i=1n(x+xA[i])P(x)=\prod\limits_{i=1}^n(x^{\varnothing}+x^{A[i]})
=expi=1nln(x+xA[i])=\exp\sum\limits_{i=1}^n\ln(x^{\varnothing}+x^{A[i]})
=expi=1nk=1(1)k+1(xA[i])kk=\exp\sum\limits_{i=1}^n\sum\limits_{k=1}\dfrac{(-1)^{k+1}(x^{A[i]})^k}{k}
(xA[i])k(x^{A[i]})^k 显然只有 k=1k=1 时为 xA[i]x^{A[i]} ,其他情况为 00。(注意这是子集卷积)
=expi=1nxA[i]=\exp\sum\limits_{i=1}^nx^{A[i]}
似乎变成 exp\exp 板子了呢……
不过要注意对 A[i]=0A[i]=0 的情况特殊处理。
不难发现“划分”就对应生成无序集合,这即为 exp\exp 的组合意义。
本题要求划分的集合数 k\leq k ,则是一个部分 exp\exp
i=0kF(x)kk!\sum\limits_{i=0}^k\dfrac{F(x)^k}{k!}
对占位多项式做部分 exp\exp 即可。
ODE\rm ODE : G=i=0kFkk!G=F(G(Fk/k!))G=\sum\limits_{i=0}^k\dfrac{F^k}{k!}\Rightarrow G'=F'*\big(G-(F^k/k!)\big)
常数莫名其妙的大,跑不过 O(2nn3)O(2^nn^3) 实锤……
G[s]G[s] 为点集 ss 的联通生成子图个数,F[s]F[s] 为点集 ss 的生成子图个数。
FF 是易求的,有 F[s]=2s内部边数F[s]=2^{\text{s内部边数}}
根据 exp\exp 的组合意义显然有 F=expGF=\exp G ,则 G=lnFG=\ln F
P[s]P[s] 为点集 ss 的生成子图的连通值之和。
枚举连通块数目则有 P=k=0Gkk!k!=11GP=\sum\limits_{k=0}\dfrac{G^k}{k!}k!=\dfrac{1}{1-G}
求逆即可。
能够发现无环图所有边反向之后还是无环图,互反的一对图翻转边数和恰为 mm。所以答案为方案数乘以 m/2m/2
F[s]F[s] 表示点集 ss 内适当反转后无环的方案数。
考虑拆 DAG\rm DAG 为子问题,每次去除入度为 00 的所有点。
钦定入度为 00 的点集 tst\subseteq stt 内部必须没有边,即 tt 为独立集。
tt 和其余点 sts-t 之间的边方向确定。所以方案数为 [t是独立集]F[st][t\text{是独立集}]F[s-t]
由于我们只能钦定一些点入度为 00 ,所以还需容斥。
可得 F[s]=ts(1)t1[t是独立集]F[st]F[s]=\sum\limits_{t\subseteq s}(-1)^{|t|-1}[t\text{是独立集}]F[s-t]
G=(1)t1[t是独立集]G=(-1)^{|t|-1}[t\text{是独立集}]
F=FG+1F=F*G+1。可得 F=11GF=\dfrac{1}{1-G},求逆即可。
可能要先去 多项式计数杂谈 看看怎么数欧拉图。
首先要求出偶度图的集合幂级数,然后求 ln\ln 即可得到欧拉图。
现在问题变为对每个点集求生成偶度图个数。
找出任意一个生成森林,将其余的边随意选取,此时能得到目前每个点的奇偶性。
由于森林无环的性质,每种电度奇偶性都能通过恰当地选取森林中的边得到,且方案数恰好为 11
因此,偶度图个数为 2边数生成森林边数2^{\text{边数}-\text{生成森林边数}}

评论

55 条评论,欢迎与作者交流。

正在加载评论...