跳转至

快速沃尔什变换

简介

沃尔什转换(Walsh Transform)1是在频谱分析上作为离散傅立叶变换的替代方案的一种方法,在信号处理中应用很广泛.FFT 是 double 类型的,但是 Walsh 把信号在不同震荡频率方波下拆解,因此所有的系数都是绝对值大小相同的整数,这使得不需要作浮点数的乘法运算,提高了运算速度.

所以,FWT 和 FFT 的核心思想应该是相同的,都是对数组的变换.我们记对数组 A 进行快速沃尔什变换后得到的结果为 FWT[A].

那么 FWT 核心思想就是:

我们需要一个新序列 C,由序列 A 和序列 B 经过某运算规则得到,即 C=A⋅B;

我们先正向得到序列 FWT[A],FWT[B],再根据 FWT[C]=FWT[A]⋅FWT[B] 在 O(n) 的时间复杂度内求出 FWT[C],其中 ⋅ 是序列对应位置相乘;

然后逆向运算得到原序列 C.时间复杂度为 O(nlog⁡n).

在算法竞赛中,FWT 是用于解决对下标进行位运算卷积问题的方法.

公式:Ci=∑i=j⊕kAjBk

(其中 ⊕ 是二元位运算中的某一种)

下面我们举 ∪(按位或)、∩(按位与)和 ⊕(按位异或)为例.

FWT 的运算

或运算

如果有 k=i∪j,那么 i 的二进制位为 1 的位置和 j 的二进制位为 1 的位置肯定是 k 的二进制位为 1 的位置的子集.

现在要得到 FWT[C]=FWT[A]⋅FWT[B],我们就要构造这个 FWT 的规则.

我们按照定义,显然可以构造 FWT[A]i=Ai′=∑i=i∪jAj,来表示 j 满足二进制中 1 为 i 的子集.

那么有:

FWT[A]i⋅FWT[B]i=(∑i∪j=iAj)(∑i∪k=iBk)=∑i∪j=i∑i∪k=iAjBk=∑i∪(j∪k)=iAjBk=FWT[C]i

那么我们接下来看 FWT[A] 怎么求.

首先肯定不能枚举了,复杂度为 O(n2).既然不能整体枚举,我们就考虑分治.

我们把整个区间二分,其实二分区间之后,下标写成二进制形式是有规律可循的.

我们令 A0 表示 A 的前一半,A1 表示区间的后一半,那么 A0 就是 A 下标最大值的最高位为 0,他的子集就是他本身的子集(因为最高位为 0 了),但是 A1 的最高位是 1,他满足条件的子集不仅仅是他本身,还包最高位为 0 的子集,即

FWT[A]=merge(FWT[A0],FWT[A0]+FWT[A1])

其中 merge 表示像字符串拼接一样把两个数组拼起来,+ 就是普通加法,表示对应二进制位相加.

这样我们就通过二分能在 O(log⁡n) 的时间复杂度内完成拼接,每次拼接的时候要完成一次运算,也就是说在 O(nlog⁡n) 的时间复杂度得到了 FWT[A].

接下来就是反演了,其实反演是很简单的,既然知道了 A0 的本身的子集是他自己(A0=FWT[A0]),A1 的子集是 FWT[A0]+FWT[A1],那就很简单的得出反演的递推式了:

UFWT[A′]=merge(UFWT[A0′],UFWT[A1′]−UFWT[A0′])

下面我们给出代码实现.容易发现顺变换和逆变换可以合并为一个函数,顺变换时 type=1,逆变换时 type=−1.

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
void Or(ll *a, ll type) {  // 迭代实现,常数更小
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j + k] += a[i + j] * type) %= P;
      }
    }
  }
}

与运算

与运算类比或运算可以得到类似结论

FWT[A]=merge(FWT[A0]+FWT[A1],FWT[A1]) UFWT[A′]=merge(UFWT[A0′]−UFWT[A1′],UFWT[A1′])

下面我们给出代码实现.顺变换时 type=1,逆变换时 type=−1.

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
void And(ll *a, ll type) {
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j] += a[i + j + k] * type) %= P;
      }
    }
  }
}

异或运算

异或的卷积是基于如下原理:

若我们令 x∘y 表示 x∩y 中 1 数量的奇偶性,即 x∘y=popcnt(x∩y)mod2,那么容易有 (x∘y)⊕(x∘z)=x∘(y⊕z).

对于 FWT[A] 的运算其实也很好得到.

设 FWT[A]i=∑i∘j=0Aj−∑i∘j=1Aj.我们来证一下 FWT[C]=FWT[A]⋅FWT[B] 的正确性:

FWT[A]iFWT[B]i=(∑i∘j=0Aj−∑i∘j=1Aj)(∑i∘k=0Bk−∑i∘k=1Bk)=(∑i∘j=0Aj∑i∘k=0Bk+∑i∘j=1Aj∑i∘k=1Bk)−(∑i∘j=0Aj∑i∘k=1Bk+∑i∘j=1Aj∑i∘k=0Bk)=∑(j⊕k)∘i=0AjBk−∑(j⊕k)∘i=1AjBk=FWT[C]i

来看看怎么快速计算 A,B 的值,依旧是分治:

对于 i 在当前位为 0 的子数列 FWT[A0],进行 ∘ 运算时发现它和 0 计算或和 1 计算结果都不会变(因为 0∩0=0,0∩1=0),所以 FWT[A]=∑i∘j=0Aj−∑i∘j=1Aj 中的 ∑i∘j=1Aj=0.

对于 i 在当前位为 1 的子数列 A1,进行 ∘ 运算时发现它和 0 计算结果是 0,和 1 计算结果是 1(因为 1∩0=0,1∩1=1).

综上,有:

FWT[A]=merge((FWT[A0]+FWT[A1])−0,FWT[A0]−FWT[A1])

也就是:

FWT[A]=merge(FWT[A0]+FWT[A1],FWT[A0]−FWT[A1])

逆变换易得:

UFWT[A′]=merge(UFWT[A0′]+UFWT[A1′]2,UFWT[A0′]−UFWT[A1′]2)

给出代码,顺变换时 type=1,逆变换时 type=12.

实现
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
void Xor(ll *a, ll type) {
  for (ll x = 2; x <= n; x <<= 1) {
    ll k = x >> 1;
    for (ll i = 0; i < n; i += x) {
      for (ll j = 0; j < k; j++) {
        (a[i + j] += a[i + j + k]) %= P;
        (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
        (a[i + j] *= type) %= P;
        (a[i + j + k] *= type) %= P;
      }
    }
  }
}

同或运算

类比异或运算给出公式:

FWT[A]i=∑C1Aj−∑C2Aj(C1 表示 popcnt(x∪y)mod2 为 0,C2 表示 popcnt(x∪y)mod2 为 1)

FWT[A]=merge(FWT[A1]−FWT[A0],FWT[A1]+FWT[A0]) UFWT[A′]=merge(UFWT[A1′]−UFWT[A0′]2,UFWT[A1′]+UFWT[A0′]2)

另一个角度的 FWT

我们设 c(i,j) 是 Aj 对 FWT[A]i 的贡献系数.我们可以重新描述 FWT 变换的过程:

FWT[A]i=∑j=0n−1c(i,j)Aj

因为有:

FWT[A]i⋅FWT[B]i=FWT[C]i

所以我们可以通过简单的证明得到:c(i,j)c(i,k)=c(i,j⊙k).其中 ⊙ 是任意一种位运算.

同时,c 函数还有一个重要的性质,它可以按位处理.

举个例子,我们变换的时候:

FWT[A]i=∑j=0n−1c(i,j)Aj

这么做是比较劣的,我们将其拆分:

FWT[A]i=∑j=0n/2−1c(i,j)Aj+∑j=n/2n−1c(i,j)Aj

考虑前面的式子和后面的式子 i,j 的区别,发现只有最高位不同.

所以我们将 i,j 去除最高位的值为 i′,j′,并记 i0 为 i 的最高位.有:

FWT[A]i=c(i0,0)∑j=0n/2−1c(i′,j′)Aj+c(i0,1)∑j=n/2n−1c(i′,j′)Aj

如果 i0=0,则有:

FWT[A]i=c(0,0)∑j=0n/2−1c(i′,j′)Aj+c(0,1)∑j=n/2n−1c(i′,j′)Aj

i0=1 则有:

FWT[A]i=c(1,0)∑j=0n/2−1c(i′,j′)Aj+c(1,1)∑j=n/2n−1c(i′,j′)Aj

也就是说,我们只需要:

[c(0,0)c(0,1)c(1,0)c(1,1)]

四个数就可以完成变换了.我们称这个矩阵为位矩阵.

如果我们要进行逆变换,则需要上面的位矩阵的逆矩阵.

若逆矩阵为 c−1,可以通过类似操作得到原数:

Ai=∑j=0nc−1(i,j)FWT[A]j

逆矩阵不一定存在,比如如果有一排 0 或者一列 0 那么这个矩阵就没有逆,我们在构造时需要格外小心.

按位或

我们可以构造:

[1011]

这样满足 c(i,j)c(i,k)=c(i,j∪k).我们发现,这和我们前面推出的 FWT[A]=merge(FWT[A0],FWT[A0]+FWT[A1]) 一模一样!同理,下面也是一个满足这个条件的矩阵,但我们一般使用上面这个:

[1110]

虽然下面这个矩阵也满足 c(i,j)c(i,k)=c(i,j∪k),但这个矩阵存在一排 0,不存在逆,所以不合法:

[0011]

如果我们要进行逆变换,则需要对矩阵求逆,以 最上面 这个矩阵为例,得:

[10−11]

然后按照顺变换的方法,把逆变换矩阵代入即可.

按位与

我们可以构造:

[1101]

这样满足 c(i,j)c(i,k)=c(i,j∩k).

逆矩阵:

[1−101]

按位异或

我们可以构造:

[111−1]

这样满足 c(i,j)c(i,k)=c(i,j⊕k).

逆矩阵:

[0.50.50.5−0.5]

FWT 是线性变换

FWT 是线性变换.也就是说,它满足:

FWT[A+B]=FWT[A]+FWT[B]

以及:

FWT[c⋅A]=c⋅FWT[A]

K 维 FWT

其实位运算的本质是对一个 n 维 {0,1} 向量的运算.或运算就是每一维取 max.且运算就是每一维取 min.异或运算则是每一维对应相加再 mod2.

位运算有个特点:向量的每一位都是独立的.

我们把 {0,1} 扩展到 [0,K)∩Z 也就是扩展到 K 进制,看看会得到什么?

max 运算

我们将 ∪ 运算拓展到 K 进制,定义 i∪j 表示按位取 max,有:

c(i,j)c(i,k)=c(i,j∪k)

若 j=k,那么上式又是:

c(i,j)c(i,j)=c(i,j)

也就是说,每一行的 1 必定只能在 0 的前面,如果在后面则不合法了.手玩一下可以发现一组合法构造:

[1000110011101111]

求逆可得:

[1000−11000−11000−11]

min 运算

我们将 ∩ 运算拓展到 K 进制,定义 i∩j 表示按位取 min,有:

c(i,j)c(i,k)=c(i,j∩k)

若 j=k,那么上式又是:

c(i,j)c(i,j)=c(i,j)

也就是说,每一行的 1 必定只能在 0 的后面,如果在前面则不合法了.手玩一下可以发现一组合法构造:

[1111011100110001]

求逆可得:

[1−10001−10001−10001]

前两者用得比较少,用得比较多的是:

不进位加法

我们将 ⊕ 运算拓展到 K 进制,定义 i⊕j 表示按位相加再 modK,有:

c(i,j)c(i,k)=c(i,j⊕k)

我们构造 c(i,j)=ωKj,就可以满足要求了:

ωKjωKk=ωKj⊕k

但是每一行都一样矩阵也没有逆,所以我们可以构造 c(i,j)=ωK(i−1)j 即可.

有下面这个矩阵:

[111⋯11ωK1ωK2⋯ωKk−11ωK2ωK4⋯ωK2(k−1)1ωK3ωK6⋯ωK3(k−1)⋮⋮⋮⋱⋮1ωKk−1ωK2(k−1)⋯ωK(k−1)(k−1)]

此即为 范德蒙德矩阵,求逆可得:

1K[111⋯11ωK−1ωK−2⋯ωK−(k−1)1ωK−2ωK−4⋯ωK−2(k−1)1ωK−3ωK−6⋯ωK−3(k−1)⋮⋮⋮⋱⋮1ωK−(k−1)ωK−2(k−1)⋯ωK−(k−1)(k−1)]

如果我们题目给出的模数是存在单位根的,我们就可以简单实现.

但是 单位根在模意义下可能不存在,所以我们考虑扩域,就是人为地定义一个 x,满足 xK=1,然后直接把 x 代入计算,这样每个数都是一个关于 x 的 k−1 次多项式.我们只需要在 modxK−1 下计算即可.那么矩阵可以这么表示:

[111⋯11x1x2⋯xk−11x2x4⋯x2(k−1)1x3x6⋯x3(k−1)⋮⋮⋮⋱⋮1xk−1x2(k−1)⋯x(k−1)(k−1)]

但是这么做可能会存在零因子,也就是 一个数有多种表示方法,我们无法确定一个数的真实值.

我们考虑不 modxK−1 了,我们 mod 分圆多项式 ΦK(x),他满足 x 的阶为 k,且在 Q 上不可约.所以我们定义上面的计算是在 modΦK(x) 下进行即可.

还有一个问题是,modΦK(x) 常数大(因为 Φ 本身就是一个多项式).但是因为 ΦK(x)∣xk−1,我们只需要在计算时 modxk−1,最后再 modΦK(x) 即可.

例题

「CF 1103E」Radix sum

给定一个长度为 n 的序列 a1,a2,...,an,对于每一个 p∈[0,n−1],求满足下列条件的整数序列 i1,i2,...,in 的方案数,对 258 取模:

  • ∀j∈[1,n],ij∈[1,n];
  • ∑j=1naij=p,这里的加法定义为十进制不进位加法.

n≤105,ai≤105

题解

我们可以想到 dp:设计状态 fi,s 表示考虑到第 i 个数,当前加法状态为 s.因为 FWT 变换是线性的,可以先变换为 FWT 点值表示法,然后变成自己的 n 次幂,最后再变换回来.

上面是平凡的,但是题目给出了模数 258.发现没有单位根,所以考虑扩域.

这里的分圆多项式 Φ10(x)=x4−x3+x2−x+1.

然而我们发现 UFWT 时,需要除去进制 10,然而我们发现 10 在 258 下没有逆元.实际上我们发现 5 在 258 下是有逆元的:57646075230342349,我们只需要再除去一个 2 就可以了.设已经除以了 5 的答案为 x,真正的答案为 y,也就是 25y≡x(mod264),显然,我们有 y≡x25(mod264−5),也就是 y≡x25(mod259),所以直接将最后的答案除以 25 即可.虽然出题人不知道为什么要模 258,但再取下模即可.

【CF103329F】【XXII Opencup, Grand Prix of XiAn】The Struggle

给出一个椭圆 E,其中所有整点的坐标均在 [1,4⋅106] 之间.求 ∑(x,y)∈E(x⊕y)33x−2y−1mod109+7 的值.

题解

这是一道比较不裸的题,出题人提供了详细的英文题解,具体请见 此链接.

参考资料