Skip to content

快速沃尔什变换(FWT)

2023-12-31

一个很邪门的算法。

二进制卷积

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

INFO

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

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=kAiBjiandj=kAiBjixorj=kAiBj

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

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

ABAorBAandBAxorB
00000
01101
10101
11110

还有高手?

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

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

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

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

FWT 也是同样的套路:

EXAMPLE
  1. 对序列 A,B 使用 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,and,xor 卷积所对应的中间态还各不相同。

OR 卷积

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

FWT(X)k=iork=kXi

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

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

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

容易发现

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

可以从集合意义理解:若 acbc,则 ac 的并集也是 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={A0,A1,,An1} 均分为两段

Aprev={A0,A1,,An21}Aback={An2,An2+1,,An1}

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

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

FWT(A)k={FWT(Aprev)k,0kn21FWT(Aprev)kn2+FWT(Aback)kn2,n2kn1

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

(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))=A,这里采用数学归纳法证明。

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

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

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

现假设命题对 n2 成立。由于 AprevAback 都是长为 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))=AbackAprevAprev=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=FWT(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=0n1f(i,k)Xi

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


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

i=0n1f(i,k)Ai×j=0n1f(j,k)Bj=x=0n1f(x,k)Cx

Cx=ixorj=xAiBj 代入

i=0n1f(i,k)Ai×j=0n1f(j,k)Bj=x=0n1f(x,k)ixorj=xAiBji=0n1j=0n1f(i,k)f(j,k)AiBj=x=0n1ixorj=xf(x,k)AiBj

0n1 范围内,对于每个可能的 i,j 组合,总有一个 x 满足 ixorj=x,所以 x=0n1ixorj=x 实际上遍历了所有的 ij,因此

i=0n1j=0n1f(i,k)f(j,k)AiBj=i=0n1j=0n1f(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=cbitcount(a)+bitcount(b) 的奇偶性和 bitcount(c) 相同。

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

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

既然奇偶性相同,那么 (1)bitcount(a)+bitcount(b)=(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=0n1(1)bitcount(iandk)Xi

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

  • 0kn21 时,直接把求和式拆分成两半:FWT(A)k=i=0n1(1)bitcount(iandk)Ai=i=0n21(1)bitcount(iandk)Ai+i=n2n1(1)bitcount(iandk)Ai=i=0n21(1)bitcount(iandk)(Aprev)i+i=0n21(1)bitcount(iandk)(Aback)i=FWT(Aprev)k+FWT(Aback)k
  • n2kn1 时,需要对下标 k 进行偏移,偏移量为 n2,理由同 OR 卷积。 需要注意的是,当 in2n1 时,bitcount(iand(kn2))=bitcount(iandk),下标偏移导致了此处符号的改变。FWT(A)k=i=0n1(1)bitcount(iandk)Ai=i=0n21(1)bitcount(iand(kn2))Aii=n2n1(1)bitcount(iand(kn2))Ai=i=0n21(1)bitcount(iand(kn2))(Aprev)ii=0n21(1)bitcount(iand(kn2))(Aback)i=FWT(Aprev)kn2FWT(Aback)kn2

FWT(A)k={FWT(Aprev)k+FWT(Aback)k,0kn21FWT(Aprev)kn2FWT(Aback)kn2,n2kn1

进一步地

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;
        }
    }
}