此划水文全为结论、板子,证明还得看大爷证明。
我们可以用FFT、NTT计算多项式乘法:
C
(
k
)
=
∑
i
+
j
=
k
A
(
i
)
∗
B
(
j
)
C(k)=\sum_{i+j=k}A(i)*B(j)
C(k)=i+j=k∑A(i)∗B(j)
但是计算不了这种玩意:
C
(
k
)
=
∑
i
∣
j
=
k
A
(
i
)
∗
B
(
j
)
C(k)=\sum_{i|j=k}A(i)*B(j)
C(k)=i∣j=k∑A(i)∗B(j)
或者这种东东:
C
(
k
)
=
∑
i
&
j
=
k
A
(
i
)
∗
B
(
j
)
C(k)=\sum_{i\&j=k}A(i)*B(j)
C(k)=i&j=k∑A(i)∗B(j)
或者这种:
C
(
k
)
=
∑
i
x
o
r
j
=
k
A
(
i
)
∗
B
(
j
)
C(k)=\sum_{i\ xor\ j=k}A(i)*B(j)
C(k)=i xor j=k∑A(i)∗B(j)
所以一个算法应着时代的需求诞生了:
多项式位运算
前言
(WTF) FWT是一个神奇的算法,它的名字叫做“快速沃尔什变换”。
虽然不知道百度百科写的是smg,但是网上博客已经足够多以至于让人学懂了。
话说我以前还不知不解地打过这个板子。现在学感觉当时浪费了这么好的一个算法。
或运算
过程
多项式或运算是用来快速做这条式子的:
C
(
i
)
=
∑
i
=
k
∣
j
A
(
j
)
∗
B
(
k
)
C(i)=\sum_{i=k|j}A(j)*B(k)
C(i)=i=k∣j∑A(j)∗B(k)
其具体思想其实和FFT差不多,首先先把多项式A转化为点值表达,然后再快速运算。
然后再转回插值。这样完成一个算法流程。然鹅我们如何快速DFT这个丑陋的式子呢?
拿出流程图:
一样三部曲:
如何找出类似于FFT的点值乘法
我们这里并不考虑去什么点值插值,我们只需要一些简单的数论变换即可得到我们所要求的东东。
同样的思路,考虑把A转化成另外一个东东,假设叫做
F
W
T
(
A
)
FWT(A)
FWT(A),其满足可以快速求出
F
W
T
(
C
)
=
F
W
T
(
A
)
∗
F
W
T
(
B
)
FWT(C)=FWT(A)*FWT(B)
FWT(C)=FWT(A)∗FWT(B)。算出来之后还可以快速转化回
A
A
A,也就是
I
F
W
T
IFWT
IFWT逆变换。
首先我们观察柿子:
i
∣
j
=
k
i|j=k
i∣j=k,我们发现,由于是或运算,看到或运算就应该拆位(这个套路真的好用)。
那么拆完位就可以得到一个结论:如果将k拆位完后,得到1位置集合。那么i的1位置集合与j的1位置集合即为k的1位置集合的子集。
得到这个结论后,我们即可构造:
F
W
T
(
A
)
=
∑
i
=
i
∣
j
A
(
j
)
FWT(A)=\sum_{i=i|j}A(j)
FWT(A)=i=i∣j∑A(j)
这条柿子的意义即为:j的1的位置集合为i的1的位置集合的子集。
(其实为什么这样构造,理由是要去进行一波反演得到的。这里就8推了,记结论)
那么如果把
F
W
T
(
A
)
FWT(A)
FWT(A)和
F
W
T
(
B
)
FWT(B)
FWT(B)乘起来会变成什么呢?
F
W
T
(
A
)
∗
F
W
T
(
B
)
=
∑
i
=
i
∣
j
A
(
j
)
∗
∑
i
=
i
∣
k
B
(
k
)
=
∑
i
=
i
∣
j
∑
i
=
i
∣
k
A
(
j
)
∗
B
(
k
)
=
∑
i
=
i
∣
(
j
∣
k
)
A
(
j
)
∗
B
(
k
)
=
F
W
T
(
C
)
FWT(A)*FWT(B)=\sum_{i=i|j}A(j)*\sum_{i=i|k}B(k)\\=\sum_{i=i|j}\sum_{i=i|k}A(j)*B(k)\\=\sum_{i=i|(j|k)}A(j)*B(k)\\=FWT(C)
FWT(A)∗FWT(B)=i=i∣j∑A(j)∗i=i∣k∑B(k)=i=i∣j∑i=i∣k∑A(j)∗B(k)=i=i∣(j∣k)∑A(j)∗B(k)=FWT(C)
那么就有方向了。如果现在已经得到一个新的东东
F
W
T
(
A
)
FWT(A)
FWT(A),那么答案就可以在
O
(
n
)
O(n)
O(n)时间内求出。
问题是如何快速求出
F
W
T
(
A
)
FWT(A)
FWT(A)。
如何找出FWT正变换
同样考虑像FFT一样分治求之。
那么递归可能会慢。怎么办?能不能像FFT一样把下标取出掘其规律?
那就试试。(逝世)
还记得这张神图把。
考虑设a0表示当前分治到的当前位为0的序列,a1表示当前分治到的当前位为1的序列。
再回顾
F
W
T
(
A
)
FWT(A)
FWT(A)的定义:
F
W
T
(
A
)
=
∑
i
=
i
∣
j
A
(
j
)
FWT(A)=\sum_{i=i|j}A(j)
FWT(A)=i=i∣j∑A(j)
那么可以知道在对应位置的a0永远是a1的子集。
所以可以写成这个式子:
F
W
T
(
A
)
=
(
F
W
T
(
A
0
)
,
F
W
T
(
A
1
+
A
0
)
)
FWT(A)=(FWT(A0),FWT(A1+A0))
FWT(A)=(FWT(A0),FWT(A1+A0))其中
A
1
+
A
0
A1+A0
A1+A0表示对应位置相加。
于是我们的FWT正变换就完成了。
如何找出FWT逆变换(IFWT)
这个直接理解还挺简单的。我们发现,由于FWT是求子集和,现在要求的是当前值。
所以把子集减去即可。
于是写成这个式子:
F
W
T
(
A
)
=
(
F
W
T
(
A
0
)
,
F
W
T
(
A
1
−
A
0
)
)
FWT(A)=(FWT(A0),FWT(A1-A0))
FWT(A)=(FWT(A0),FWT(A1−A0))
后话
话说其实FWT可以有多种方法来看。
网上有的说这其实是FMT(快速莫比乌斯变换)(其实这玩意我真的没搞懂它和反演有什么区别),我一开始看Vfleaking的反演觉得是个子集反演用分治来做。
然鹅比较简单的理解就是从位运算的意义上理解,可以推柿子来比较严谨地去证明,当然也可以我这样直接理解地去推结论。
一个问题有多种角度去看,有时是一件好事,有时却是一件坏事,因为可能会莫衷一是,一开始学会陷入理解的困境。
不管怎样,还是挺好玩的。就好像3blue1brown让我学会的一个道理,去发现数学的美,而不是去死板地认为这个就是一个工具。
诶我tm怎么莫名其妙就放了怎么多没用的屁话
结论
F
W
T
(
A
)
=
(
F
W
T
(
A
0
)
,
F
W
T
(
A
1
+
A
0
)
)
FWT(A)=(FWT(A0),FWT(A1+A0))
FWT(A)=(FWT(A0),FWT(A1+A0))
I
F
W
T
(
A
)
=
(
I
F
W
T
(
A
0
)
,
I
F
W
T
(
A
1
−
A
0
)
)
IFWT(A)=(IFWT(A0),IFWT(A1-A0))
IFWT(A)=(IFWT(A0),IFWT(A1−A0))
代码
void orfwt(int a[],int inv)
{
long long op,oq;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
op=inv*a[k];
oq=a[k+mid];
op=(op+oq+mo)%mo;
a[k+mid]=op;
}
}
}
}
和运算
过程
这个过程和或运算的过程其实长得基本一样,就不再推一遍了。
结论
F
W
T
(
A
)
=
(
F
W
T
(
A
0
+
A
1
)
,
F
W
T
(
A
1
)
)
FWT(A)=(FWT(A0+A1),FWT(A1))
FWT(A)=(FWT(A0+A1),FWT(A1))
I
F
W
T
(
A
)
=
(
I
F
W
T
(
A
0
−
A
1
)
,
I
F
W
T
(
A
1
)
)
IFWT(A)=(IFWT(A0-A1),IFWT(A1))
IFWT(A)=(IFWT(A0−A1),IFWT(A1))
代码
void andfwt(int a[],int inv)
{
long long op,oq;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
op=a[k];
oq=inv*a[k+mid];
op=(op+oq+mo)%mo;
a[k]=op;
}
}
}
}
异或运算
过程
这个其实和上面的思路差不多,但是构造的东东还是很不同的。
如何找出类似于FFT的点值乘法
考虑设一个函数
f
(
x
)
f(x)
f(x)表示当前
x
x
x在二进制表示下1的数量的奇偶性。
则有:
F
W
T
(
A
)
=
∑
f
(
i
&
j
)
=
0
A
(
j
)
−
∑
f
(
i
&
j
)
=
1
A
(
j
)
FWT(A)=\sum_{f(i\&j)=0}A(j)-\sum_{f(i\&j)=1}A(j)
FWT(A)=f(i&j)=0∑A(j)−f(i&j)=1∑A(j)
这玩意当然满足
F
W
T
(
C
)
=
F
W
T
(
A
)
∗
F
W
T
(
B
)
FWT(C)=FWT(A)*FWT(B)
FWT(C)=FWT(A)∗FWT(B)
理由就是把这条柿子带进去,然后再瞎换一下即可得到。
当然,在推柿子的时候可能要用到一个结论:
f
(
i
&
j
)
x
o
r
f
(
i
&
k
)
=
f
(
i
&
(
j
x
o
r
k
)
)
f(i\&j)\ xor\ f(i\&k)=f(i\&(j\ xor\ k))
f(i&j) xor f(i&k)=f(i&(j xor k))
这里8推了。
如何找出FWT正变换
考虑如何求
F
W
T
(
A
)
FWT(A)
FWT(A)
观察柿子:
F
W
T
(
A
)
=
∑
f
(
i
&
j
)
=
0
A
(
j
)
−
∑
f
(
i
&
j
)
=
1
A
(
j
)
FWT(A)=\sum_{f(i\&j)=0}A(j)-\sum_{f(i\&j)=1}A(j)
FWT(A)=f(i&j)=0∑A(j)−f(i&j)=1∑A(j)
由于是位运算,那么继续套路,拆位。
考虑当前第i位,这一位就有4种情况:
- 0 xor 0=0,奇偶性不会改变。
- 0 xor 1=1,奇偶性不会改变。
- 1 xor 0=1,奇偶性不会改变。
- 1 xor 1=0,这时奇偶性改变。
这意味着什么呢?
假如把
F
W
T
(
A
)
FWT(A)
FWT(A)取个绝对值(当然这样计算最后的答案是不会改变的),那么就可以把奇偶性改变都看做是减去贡献,反之则为加上贡献。
那么正变换就得到了:
F
W
T
(
A
)
=
(
F
W
T
(
A
0
+
A
1
)
,
F
W
T
(
A
0
−
A
1
)
)
FWT(A)=(FWT(A0+A1),FWT(A0-A1))
FWT(A)=(FWT(A0+A1),FWT(A0−A1))
如何找出FWT逆变换(IFWT)
逆变换一样简单,把上面的计算贡献反过来即可。
I
F
W
T
(
A
)
=
(
I
F
W
T
(
A
0
+
A
1
)
2
,
I
F
W
T
(
A
0
−
A
1
)
2
)
IFWT(A)=(\frac{IFWT(A0+A1)}2,\frac{IFWT(A0-A1)}2)
IFWT(A)=(2IFWT(A0+A1),2IFWT(A0−A1))
结论
F
W
T
(
A
)
=
(
F
W
T
(
A
0
+
A
1
)
,
F
W
T
(
A
0
−
A
1
)
)
FWT(A)=(FWT(A0+A1),FWT(A0-A1))
FWT(A)=(FWT(A0+A1),FWT(A0−A1))
I
F
W
T
(
A
)
=
(
I
F
W
T
(
A
0
+
A
1
)
2
,
I
F
W
T
(
A
0
−
A
1
)
2
)
IFWT(A)=(\frac{IFWT(A0+A1)}2,\frac{IFWT(A0-A1)}2)
IFWT(A)=(2IFWT(A0+A1),2IFWT(A0−A1))
代码
void xorfwt(int a[],int inv)
{
long long op,oq,kk;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
if (inv==1)
{
op=a[k];
oq=a[k+mid];
kk=(op+oq+mo)%mo;a[k]=kk;
kk=(op-oq+mo)%mo;a[k+mid]=kk;
}
else
{
op=a[k];
oq=a[k+mid];
kk=(op+oq+mo)%mo*inv2%mo;a[k]=kk;
kk=(op-oq+mo)%mo*inv2%mo;a[k+mid]=kk;
}
}
}
}
}
总代码
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
using namespace std;
const long long mo=998244353;
const long long inv2=499122177;
int n,m,a[131072],b[131072],ja[131072],jb[131072];
void orfwt(int a[],int inv)
{
long long op,oq;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
op=inv*a[k];
oq=a[k+mid];
op=(op+oq+mo)%mo;
a[k+mid]=op;
}
}
}
}
void andfwt(int a[],int inv)
{
long long op,oq;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
op=a[k];
oq=inv*a[k+mid];
op=(op+oq+mo)%mo;
a[k]=op;
}
}
}
}
void xorfwt(int a[],int inv)
{
long long op,oq,kk;
for (int len=2;len<=m;len<<=1)
{
int mid=len/2;
for (int j=0;j<mid;j++)
{
for (int k=j;k<m;k+=len)
{
if (inv==1)
{
op=a[k];
oq=a[k+mid];
kk=(op+oq+mo)%mo;a[k]=kk;
kk=(op-oq+mo)%mo;a[k+mid]=kk;
}
else
{
op=a[k];
oq=a[k+mid];
kk=(op+oq+mo)%mo*inv2%mo;a[k]=kk;
kk=(op-oq+mo)%mo*inv2%mo;a[k+mid]=kk;
}
}
}
}
}
void solve(int a[],int b[],int ki)
{
if (ki==1) orfwt(a,1),orfwt(b,1);
else if (ki==2) andfwt(a,1),andfwt(b,1);
else xorfwt(a,1),xorfwt(b,1);
long long op,oq;
for (int i=0;i<m;i++)
{
op=a[i];
oq=b[i];
op=op*oq%mo;
a[i]=op;
}
if (ki==1) orfwt(a,-1);
else if (ki==2) andfwt(a,-1);
else xorfwt(a,-1);
for (int i=0;i<m;i++)
{
printf("%d ",a[i]);
}
printf("\n");
}
int main()
{
scanf("%d",&n);
m=1<<n;
for (int i=0;i<m;i++)
{
scanf("%d",&ja[i]);
}
for (int i=0;i<m;i++)
{
scanf("%d",&jb[i]);
}
for (int i=1;i<=3;i++)
{
memcpy(a,ja,sizeof(a));
memcpy(b,jb,sizeof(b));
solve(a,b,i);
}
}
学习资料:
http://oi-wiki.com/math/poly/fwt/
https://blog.csdn.net/zhouyuheng2003/article/details/85950280
https://blog.csdn.net/hzj1054689699/article/details/83340154
https://www.cnblogs.com/cjyyb/p/9065615.html
http://blog.leanote.com/post/rockdu/TX20
https://blog.csdn.net/zhouyuheng2003/article/details/84728063
https://zhuanlan.zhihu.com/p/41867199
https://www.cnblogs.com/wjyyy/p/FWT.html