位运算卷积(FWT) & 集合幂级数
command_block · · 算法·理论
根据
根据
对比左右两边,得到
另外,由于位运算每一位的独立性,
假设我们已经(根据真值表)求出满足要求的
设二进制数
则有
那么 : 对于每个
这就符合我们的条件。
好了,现在假设我们已经有了符合要求的
根据这个式子直接求和,复杂度至少是
我们考虑按位拆半:
设
在首位分开考虑 :
考虑到
设
若
若
我们就能以
所以,若
( 可能有点抽象,但是您如果写过FFT,看到代码就会懂了 )
此外,逆变换
( 一个重要的地方是,这个构造出来的
0x02.基础位运算卷积
针对不同的位运算,根据
我们把这个矩阵称为位矩阵。
构造的过程可能有些繁杂,可以直接记结论,或者去后面看扩展版的。
1.1\ \rm Or 卷积
设位矩阵为
起点 :
-
c(0,0)c(0,0)=c(0,0|0)
同理不难推知
-
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那么这个矩阵就没有逆,那么可以构造出:
两种情况 :
Tips :
Or 卷积的上面第二个矩阵FWT 相当于子集求和。原因:第二个矩阵相当于
c(i,j)=[i\&j=j] (也可以使用高维前缀和来推导)
(下面采用第二个矩阵)
对于逆变换,把矩阵求个逆可得
1.2\ \rm And 卷积
起点 :
同上,容易得到
-
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
还是老样子,如果有一排
两种情况 :
把矩阵求个逆可得
1.3\ \rm Xor 卷积
起点 :
-
对于任意的
x,y ,均有c(0,0)c(x,y)=c(x,y\ \text{xor}\ 0)=c(x,y) -
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\}
两种情况 :
附 :不难观察出
求逆可得
1.4 模板题 & Code:
P4717 【模板】快速沃尔什变换 (FWT)
- FWT :
评测记录
#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;
}
- 分治乘法 :
留个代码以备将来研究……
与前后文无关。
#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 变换的一些性质
前面讲了,
这意味着
此外,如果需要卷积的向量只有少数非
例题(我只能做到这些了):
CF449D Jzzhu and Numbers + 题解Link
CF1119H Triple + 题解Link
题做多了再来填坑吧。
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
0x04.位值域的扩展
其实位运算的本质是对一个
或运算就是每一维取
异或运算则是每一维对应相加再
位运算有个特点 : 向量的每一位都是独立的。
我们把
每一维取 \bf max
对应原来的按位
可得
所以整个矩阵中只有
又可得
假如有
又因为
也就是说,每一行以内,
接下来,除了矩阵有逆之外,再没有别的要求了。
如果要有逆的话,每一行都必须不同,那么
例子 :
暴力做的话
这个矩阵本质上就是前缀和和差分……所以可以优化到
(宏观上就是高维前缀和)
每一维取 \bf min
对应原来的按位
类似的,整个矩阵中只有
例子 :
这个矩阵本质上就是后缀和和差分……同样能优化到
(宏观上就是高维后缀和)
每一维加起来 \bf mod\ k
对应原来的按位
这玩意就比较复杂了,需要动用单位根……
我们知道
但是,每一行都一样的话,矩阵就没有逆。
你会发现直接把 FFT
的范德蒙德矩阵拿来用就好了:
逆矩阵我们也知道,那就是:
暴力计算显然是
可以使用
但是,单位根在模意义很可能不存在。考虑扩充我们“数”的表示。
首先想到的是,人为地定义代数对象
此举相当于用
我们的位矩阵长这样 :
接下来的推导可能需要抽象代数和高等代数相关知识。详情可见 : 近世代数乱编
问题在于,此时
这样就会导致两个位矩阵并非严格互逆。
具体而言,我们原本希望
当
当
若
但是,此时由于无法做除法,等比数列求和公式不再生效,无法证明上式
但是,仍然恒有
也就是说,我们算得的结果是由 (真实答案)+(零因子的线性组合) 构成的。
接下来要找到一种合适的扩域方法,这样就能避免零因子问题。
不妨取分圆多项式
(相关知识可见《近世代数乱编》)
-
① 分圆多项式在
Q 上不可约。这保证了
\pmod {\Phi_k(x)} 意义下的“数”是个域,不存在零因子。这里注意,还要验证该多项式在
F_p 下是否能够分解,一般情况下是不能的,此时可以正常使用。若计算时没有利用
F_p 的性质(如求逆元),则可以看作大整数计算,最后将答案取模,此时不需要检查F_p 下是否能够分解。 -
② 在
\pmod {\Phi_k(x)} 意义下,x 的阶恰好为k 。这又保证了我们构造的前提成立。
这样,只需在
但是,模多项式的计算较为繁琐,常数较大且不易优化。
而更好的消息是,由于
这表明我们可以先用“循环卷积多项式”代替严格的扩域,到最后再对
此时任意两个数的乘法就变成
观察矩阵可得,每次乘的都是一个单项式,复杂度就只需升高
暴力做的话,总的复杂度是
可以用
-
例题① : CF1103E Radix sum
-
例题② : P5577 [CmdOI2019]算力训练
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
0x04.子集卷积 & 集合幂级数
- P6097 【模板】子集卷积
作用 : 求
另一种说法 :
组合意义就是每次挑选不交的两个集合拼成新的集合。
好像在某些
使用枚举子集的技巧即可做到
直接套用前面学习的位运算卷积只能求
考虑到
( 其中
我们可以将原来的
这样,计算
我们取出
考虑计算的复杂度,多项式加法和数乘的复杂度是
中间还需要点乘,多项式乘法用暴力实现,复杂度也是
在具体实现中,可以把
#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
这就是一道很好的模板题 (不过值域只有
令
求出
令
最后把
-
半在线子集卷积
给出集合幂级数
W,C ,求幂级数S 使得 :S[s]=C[s]\sum\limits_{\small t\subsetneq s}S[t]W[s-t]
对于
我们可以按照
此时为了方便需要将
例题 : P4221 [WC2018]州区划分 | Solution
-
集合幂级数运算进阶
需熟知如何在
O(n^2) 内做一系列多项式操作,可见 NTT与多项式全家桶。下面以
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 提高组]优秀子序列
不难发现,对于序列中的
考虑模仿快速
而
似乎变成
不过要注意对
- 例题② : Loj#154. 集合划分计数
不难发现“划分”就对应生成无序集合,这即为
本题要求划分的集合数
对占位多项式做部分
有
常数莫名其妙的大,跑不过
评测记录
- 例题③ : Uoj#94. 【集训队互测2015】胡策的统计
设
由于我们只能钦定一些点入度为
可得
设
即
评测记录
- 例题⑤ : #6673. EntropyIncreaser 与山林
可能要先去 多项式计数杂谈 看看怎么数欧拉图。
首先要求出偶度图的集合幂级数,然后求
现在问题变为对每个点集求生成偶度图个数。
找出任意一个生成森林,将其余的边随意选取,此时能得到目前每个点的奇偶性。
由于森林无环的性质,每种电度奇偶性都能通过恰当地选取森林中的边得到,且方案数恰好为
因此,偶度图个数为
评测记录