Skip to content

快速沃尔什变换(FWT)

2023-12-31

一个很邪门的算法。

二进制卷积 ​

快速傅里叶变换可以高效地计算两个多项式的乘积,进而加速了卷积运算。

INFO

序列 A 和序列 B 的卷积记作 C=00A∗B,其中

Ck=∑i+j=kAiBj

构造多项式

A(x)=A0+A1x+A2x2+⋯B(x)=B0+B1x+B2x2+⋯

我们知道 A(x)⋅B(x) 中 xk 前的系数就是 Ck,而 A(x)⋅B(x) 的系数序列可由 FFT 得到。因此 FFT 也被认为是快速卷积算法。

现在考虑一类奇怪的卷积:

∑iorj=kAiBj∑iandj=kAiBj∑ixorj=kAiBj

其中 or,00and,00xor 是二进制按位或、与、异或。下图附上了一位二进制运算的真值表。

面对这类奇葩的卷积时,FFT 就束手无策了。

ABAorBAandBAxorB
00000
01101
10101
11110

还有高手? ​

快速沃尔什变换(Fast Walsh-Hadamard Transform,FWT)是加速 二进制卷积 运算的算法,其大方向与 FFT 非常相似。

回忆一下 FFT 如何加速多项式乘法:

EXAMPLE
  1. 对 A(x) 和 B(x) 使用 DFT,得到两个点集:{(xi,A(xi))∣0≤i<N}{(xi,B(xi))∣0≤i<N}
  2. 计算出 C(x)=00A(x)⋅B(x) 的点值表示法:{(xi,A(xi)⋅B(xi))∣0≤i<N}
  3. 使用 IDFT 将其转化为系数表示法。

我们在 DFT 和 IDFT 中使用 FFT 的思想来加速这个过程。

FWT 也是同样的套路:

EXAMPLE
  1. 对序列 A,00B 使用 FWT,得到两个中间态(中间态也是序列):FWT(A),FWT(B)
  2. 使用点对点乘法,计算出 C 的中间态:FWT(C)={FWT(A)i×FWT(B)i | i=0,1,2,⋯}
  3. 使用 UFWT(也称 IFWT)将其还原为序列 C。

只不过 FWT 的中间态略显抽象,并且 or,00and,00xor 卷积所对应的中间态还各不相同。

OR 卷积 ​

在 or 卷积中,序列 X 对应的中间态为与之等长的序列 FWT(X),其中

FWT(X)k=∑iork=kXi

🤔沃尔什究竟是怎么想出来的?

设有两个序列 A,00B,待求的序列是 C,其中 Ck=00∑iorj=kAiBj。

可以证明,对 A 和 B 的中间态进行点对点乘法可以得出 C 的中间态。

证

容易发现

aorc=c,borc=c⇒(aorb)orc=c

可以从集合意义理解:若 a∈c 且 b∈c,则 a 和 c 的并集也是 c 的子集。

接下来的证明步骤会用到该结论。

 FWT(A)k×FWT(B)k= (∑iork=kAi)(∑jork=kBj)= ∑iork=k,jork=kAiBj= ∑(iorj)ork=kAiBj=====x=iorj=∑xork=k ∑iorj=xAiBj= ∑xork=kCx= FWT(C)k

证毕。

现在推导一下中间态的求法。

类似于 FFT,此处也规定所有序列的长度 n 都是 2 的幂次。

我们将序列 A=00{A0,A1,⋯,An−1} 均分为两段

Aprev={A0,A1,⋯,An2−1}Aback={An2,An2+1,⋯,An−1}

然后仔细盯真 A 的中间态:

FWT(A)k=∑iork=kAi
  • 当 0≤k≤n2−001 时,由 or 运算的性质可知,iork=00k 中 i 的大小不可能超过 k。因此 i 也落在 0⋯n2−001 的范围内,也就是说 Ai 只能来自 Aprev,即 ∑iork=kAi=00∑iork=k(Aprev)i,即FWT(A)k=FWT(Aprev)k
  • 当 n2≤k≤n−001 时,Ai 既可以来自 Aprev,又可以来自 Aback。由于 FWT(A)k 是一个求和式,我们可以顺理成章地把它拆成两个部分的和FWT(A)k=FWT(Aprev)k−n2+FWT(Aback)k−n2注:此处对 k 减去 n2 是为了调整下标,使其适应相应的分段序列。或者说,Aprev 和 Aback 的长度都只有 n2,作此调整可以防止数组访问越界。

即

FWT(A)k={FWT(Aprev)k,0≤k≤n2−1FWT(Aprev)k−n2+FWT(Aback)k−n2,n2≤k≤n−1

更进一步地,以序列为基本单位看问题:

(1)FWT(A)={A,n=1merge(FWT(Aprev),FWT(Aback)⊕FWT(Aprev)),n>1

其中 merge 是拼接两个序列的操作,⊕ 表示两个序列的点对点加法。


上面我们已经剖析了 FWT 的算法,现在还需要讨论其反演 UFWT。

其实只要把点对点加法 ⊕ 改为点对点减法 ⊖ 就可以得到 UFWT 的递推式。

(2)UFWT(A)={A,n=1merge(UFWT(Aprev),UFWT(Aback)⊖UFWT(Aprev)),n>1

容易证明 (2) 式是 (1) 式的逆变换。

证

此论断等价于 UFWT(FWT(A))=00A,这里采用数学归纳法证明。

设序列 A 的长度为 n。当 n=001 时命题显然成立。

由于 n 始终是 2 的幂次,即只需证明当 n>1 时

命题对n2 成立⇒命题对n 成立

现假设命题对 n2 成立。由于 Aprev 和 Aback 都是长为 n2 的序列,因此有

UFWT(FWT(Aprev))=Aprev,UFWT(FWT(Aback))=Aback

由式 (1) (2) 得(将 merge 拆开):

{FWT(A)prev=FWT(Aprev)FWT(A)back=FWT(Aback)⊕FWT(Aprev){UFWT(A)prev=UFWT(Aprev)UFWT(A)back=UFWT(Aback)⊖UFWT(Aprev)

因此

UFWT(FWT(A))prev=UFWT(FWT(A)prev)=UFWT(FWT(Aprev))=AprevUFWT(FWT(A))back=UFWT(FWT(A)back)⊖UFWT(FWT(A)prev)=UFWT(FWT(Aback)⊕FWT(Aprev))⊖UFWT(FWT(Aprev))=UFWT(FWT(Aback))⊕UFWT(FWT(Aprev))⊖UFWT(FWT(Aprev))=Aback⊕Aprev⊖Aprev=Aback

因此

UFWT(FWT(A))=merge(Aprev,Aback)=A

即命题对 n 成立。证毕。

模板 ​

cpp
void FWT_OR(vector<int>& a, int l, int r, bool inv) {
    if (l == r)
        return;
    int mid = (l + r) >> 1;
    FWT_OR(a, l, mid, inv);
    FWT_OR(a, mid + 1, r, inv);
    for (int i = 0; i <= mid - l; i ++) {
        if (!inv)
            a[mid + 1 + i] += a[l + i];
        else
            a[mid + 1 + i] -= a[l + i];
    }
}

AND 卷积 ​

and 卷积和 or 卷积在本质上是类似的。

FWT(X)k=∑iandk=kXi

类似可证 FWT(A)k×FWT(B)k=00FWT(C)k。

由于 and 运算的性质,递推式的点对点加减法操作都集中在前面。

FWT(A)={A,n=1merge(FWT(Aprev)⊕FWT(Aback),FWT(Aback)),n>1UFWT(A)={A,n=1merge(UFWT(Aprev)⊖UFWT(Aback),UFWT(Aback)),n>1

模板 ​

cpp
void FWT_AND(vector<int>& a, int l, int r, bool inv) {
    if (l == r)
        return;
    int mid = (l + r) >> 1;
    FWT_AND(a, l, mid, inv);
    FWT_AND(a, mid + 1, r, inv);
    for (int i = 0; i <= mid - l; i ++) {
        if (!inv)
            a[l + i] += a[mid + 1 + i];
        else
            a[l + i] -= a[mid + 1 + i];
    }
}

XOR 卷积 ​

xor 卷积的构造相比前两个更加麻烦。

设序列 X(长为 n)的中间态为 FWT(X),其中

FWT(X)k=∑i=0n−1f(i,k)Xi

其中 f(x,00y) 是待确定的函数,需要通过以下过程反推出来。


要使 FWT(A)k×FWT(B)k=00FWT(C)k,即

∑i=0n−1f(i,k)Ai×∑j=0n−1f(j,k)Bj=∑x=0n−1f(x,k)Cx

将 Cx=00∑ixorj=xAiBj 代入

∑i=0n−1f(i,k)Ai×∑j=0n−1f(j,k)Bj=∑x=0n−1f(x,k)∑ixorj=xAiBj∑i=0n−1∑j=0n−1f(i,k)f(j,k)AiBj=∑x=0n−1∑ixorj=xf(x,k)AiBj

在 0⋯n−001 范围内,对于每个可能的 i,00j 组合,总有一个 x 满足 ixorj=00x,所以 ∑x=0n−1∑ixorj=x 实际上遍历了所有的 i 和 j,因此

∑i=0n−1∑j=0n−1f(i,k)f(j,k)AiBj=∑i=0n−1∑j=0n−1f(ixorj,k)AiBj

即得

(i)f(i,k)f(j,k)=f(ixorj,k)

注意到

(ii)axorb=c⇒(−1)bitcount(a)(−1)bitcount(b)=(−1)bitcount(c)

其中 bitcount(x) 表示 x 的二进制中 1 的个数。

TIP

(ii) 式是在说:axorb=00c,bitcount(a)+00bitcount(b) 的奇偶性和 bitcount(c) 相同。

这个结论是显然的。在二进制按位异或的过程中:

  • 若 a,00b 当前位均为 1,则 c 的该位为 0,相当于少了两个 1,bitcount 的奇偶性没变;
  • 若一个是 1,另一个是 0,则 c 的该位为 1,相当于没差;
  • 若都是 0 也没差。

既然奇偶性相同,那么 (−1)bitcount(a)+bitcount(b)=00(−1)bitcount(c),变形可得 (ii) 式。

又注意到

(iandk)xor(jandk)=(ixorj)andk

代入 (ii) 式得

(−1)bitcount(iandk)(−1)bitcount(jandk)=(−1)bitcount((ixorj)andk)

于是便可以确定函数 f 的形式:

f(x,y)=(−1)bitcount(xandy)

该形式完美符合 (i) 式的结论。

从而序列 X 的中间态的全貌应该是这样的:

FWT(X)k=∑i=0n−1(−1)bitcount(iandk)Xi

现对 FWT(A)k 进行分治。

  • 当 0≤k≤n2−001 时,直接把求和式拆分成两半:FWT(A)k=∑i=0n−1(−1)bitcount(iandk)Ai=∑i=0n2−1(−1)bitcount(iandk)Ai+∑i=n2n−1(−1)bitcount(iandk)Ai=∑i=0n2−1(−1)bitcount(iandk)(Aprev)i+∑i=0n2−1(−1)bitcount(iandk)(Aback)i=FWT(Aprev)k+FWT(Aback)k
  • 当 n2≤k≤n−001 时,需要对下标 k 进行偏移,偏移量为 n2,理由同 OR 卷积。 需要注意的是,当 i 取 n2⋯n−001 时,bitcount(iand(k−n2))=00−00bitcount(iandk),下标偏移导致了此处符号的改变。FWT(A)k=∑i=0n−1(−1)bitcount(iandk)Ai=∑i=0n2−1(−1)bitcount(iand(k−n2))Ai−∑i=n2n−1(−1)bitcount(iand(k−n2))Ai=∑i=0n2−1(−1)bitcount(iand(k−n2))(Aprev)i−∑i=0n2−1(−1)bitcount(iand(k−n2))(Aback)i=FWT(Aprev)k−n2−FWT(Aback)k−n2

即

FWT(A)k={FWT(Aprev)k+FWT(Aback)k,0≤k≤n2−1FWT(Aprev)k−n2−FWT(Aback)k−n2,n2≤k≤n−1

进一步地

FWT(A)={A,n=1merge(FWT(Aprev)⊕FWT(Aback),FWT(Aprev)⊖FWT(Aback)),n>1

对应地

UFWT(A)={A,n=1merge(UFWT(Aprev)⊕UFWT(Aback)2,UFWT(Aprev)⊖UFWT(Aback)2),n>1

其中对序列的除法运算也是点对点除法。

模板 ​

cpp
void FWT_XOR(vector<int>& a, int l, int r, bool inv) {
    if (l == r)
        return;
    int mid = (l + r) >> 1;
    FWT_XOR(a, l, mid, inv);
    FWT_XOR(a, mid + 1, r, inv);
    for (int i = 0; i <= mid - l; i ++) {
        int u = a[l + i], v = a[mid + 1 + i];
        a[l + i] = u + v;
        a[mid + 1 + i] = u - v;
    }
    if (inv) {
        for (int i = l; i <= r; i ++) {
            a[i] /= 2;
        }
    }
}