Fast Fourier Transforms

概述

FFT是数论中的一个非常有趣的算法，也是许多算法的关键组成部分，包括大整数乘法、多项式乘法和擦除码的生成和恢复

除了容错数据存储和恢复中有在用之外，擦除码还可以确保可扩展区块链和STARK中的数据可用性

背景

原本的傅里叶变换是一个数学操作，将数据从频域转换到时域

这意味着如果你有一些数据，运行傅里叶变换可以得到一系列不同频率和振幅的正弦波，将这些正弦波相加可以近似的得到原始数据

傅里叶变化可以用于生成一个和大象有关的方程

也就是说，反向使用傅里叶变换只需要将得到的正弦波相加，就可以在想要采样的任意多个点上计算结果

本文讨论的傅里叶变换是一种类似的算法，但区别在于他不是实数或复数域上连续的傅里叶变换，而是有限域上离散的傅里叶变换，且主要讨论下列两种不同的操作，而非原本傅里叶变换中的频域和时域之间的转换

• 多点多项式评估：评估$n$个不同的点构成的度小于$n$的多项式及其逆运算
• 多项式插值：给定$n$个不同的点，恢复一个度小于$n$的多项式

举个例子，比如我们想在模5的域上进行计算，假设我们有一个多项式$f(x)=x^2+3$

出于方便，我们可以仅保存多项式的系数，按$x$的幂次由低到高排列，我们得到系数向量$[3,0,1]$，此时我们在三个不同的点$[0,1,2]$计算多项式的值，可以得到$[3,4,2]$（注意模5的域上计算）

利用向量$[0,1,2]$和$[3,4,2]$，我们可以恢复出$f(x)$的系数向量为$[3,0,1]$

该方法的复杂度为$O(n^2)$，为点的个数的平方阶，具体的伪代码如下

def eval_poly_at(self, poly, x, modulus):
    y = 0
    power_of_x = 1
    for coefficient in poly:
        y += power_of_x * coefficient
        power_of_x *= x
    return y % modulus

每次循环需要$n$次操作，一共有$n$个点，因此$n$次循环，复杂度为$O(n^2)$

不过拉式插值法更为复杂一些，我们可以这样，对于一个域$D$和一个点$x$，我们可以构造一个多项式，使得该多项式在点$x$处的值为1，在其他点处的值均为0

比如$D=\{1,2,3,4\}$和$x=1$，我们可以构造这样一个多项式

$$ P(x)=\frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)} $$

不难看出，$P(1)=1$，$P(2)=P(3)=P(4)=0$

此外我们还可以通过将这些多项式相乘和相加来恢复给定域上任何期望输出集的多项式

比如上面这个多项式记为$P_1$，在$x=2,x=3,x=4$处等于1的多项式分别记为$P_2,P_3,P_4$，这样我们可以构造一个新的多项式

$$ P=3P_1+P_2+4P_3+P_4 $$

该多项式在$D=\{1,2,3,4\}$上的输出向量为$[3,1,4,1]$

构造这样的多项式期望复杂度为$O(n^2)$，因为构造一个$P_i$需要$O(n^2)$，计算线性组合也需要$O(n^2)$，总的期望为$O(n^2)$

快速傅里叶变换的目的是让上面这个过程更快一些

快速傅里叶变换

虽然快，但肯定是有代价的，代价是不能选择任意的域和集合，不过通过拉式插值法，可以选择任意需要的$x$和$y$的坐标以及任意想要的域，并得到一个与这些值对应的多项式

对于FFT而言，必须使用有限域，且选择的集合必须是这个有限域一个乘法子群（比如某个生成元的幂生成的群）

举个例子，模$p=337$的素数域，选择子群$D=\{1,85,148,111,336,252,189,226\}$，该群的生成元为85

此外还有一个要求，子群的大小必须为$2^n$，上面这个大小为8就很好，如果是模$p=59$的素域，其只有一个乘法子群的大小为3，显然不为$2^n$

Note：其实大小为$2^n$是为了构建一个递归算法（模快速幂）来加速计算

比如说我们想要用FFT来计算某个集合中$x$的2的次幂，其中$x^{2^k}=1$，如果用前面的那个模337的子群$D=\{1,85,148,111,336,252,189,226\}$，我们可以得到$85^{2^3}=85^8=1$

此时假设我们有一个多项式$P(x)=6x^7+2x^6+9x^5+5x^4+x^3+4x^2+x+3$，用系数向量表示为$p=[3,1,4,1,5,9,2,6]$

接下来我们希望在域上所有的点计算该多项式的值，首先将这个多项式分为两部分，$x$的幂为奇数和偶数两个部分，分别记为$P_{odd}=6x^3+9x^2+x+1$和$P_{even}=2x^3+5x^2+4x+3$，对应的系数向量分别为$p_{odd}=[1,1,9,6],p_{even}=[3,4,5,2]$

此时我们利用$P_{odd}$和$P_{even}$可以恢复出原始的多项式$P$

$$ \begin{aligned}P(x)&=x\cdot P_{odd}(x^2)+P_{even}(x^2)\\P(-x)&=-x\cdot P_{odd}(x^2)+P_{even}(x^2)\end{aligned} $$

不难看出，$P_{odd}$和$P_{even}$的大小均为多项式$P$的一半，且需要计算的点由原来的$x$变为了$x^2$，这意味着点的数量也减少了一半

然后观察子群$D=\{1,85,148,111,336,252,189,226\}$，不难看出里面有四对在模337下互为相反数，1和336，85和252，148和189，111和226，也即有四对$(x,-x)$，且这四对相反数的平方相同

因此我们利用FFT计算所有的$P_{even}(x)$，只需要计算一个更小的集合$\{1,148,189,336\}$即可，$P_{odds}$同理，从而我们将问题的规模减半了

给一个简单的Python代码，如下

def fft(vals, modulus, domain):
    if len(vals) == 1:
        return vals
    L = fft(vals[::2], modulus, domain[::2])
    R = fft(vals[1::2], modulus, domain[::2])
    o = [0 for i in vals]
    for i, (x, y) in enumerate(zip(L, R)):
        y_times_root = y*domain[i]
        o[i] = (x+y_times_root) % modulus
        o[i+len(L)] = (x-y_times_root) % modulus
    return o

然后根据上述例子，输入测试样例

>>> fft([3,1,4,1,5,9,2,6], 337, [1, 85, 148, 111, 336, 252, 189, 226])
[31, 70, 109, 74, 334, 181, 232, 4]

这样一来，我们可以利用FFT快速的计算数的乘法，举个例子，比如我们想计算$1253\times1895$，我们可以将这两个数视为两个多项式，系数向量分别为$[3,5,2,1],[5,9,8,1]$，然后利用前面的FFT算法计算两个数的乘积，最后将得到的多项式的系数向量再转换为数字就行

二元域的FFT

事实上不只有素数域，比如二元域也是可以的

二元域中只有元素0和1，比如$f(x)=x^3+x+1$，所有二元域上的加减乘法都需要模2

举一些简单的例子，比如模二元域上的不可约多项式$f(x)=x^4+x+1$，计算$(x^2+1)\cdot (x^3+1)$，计算结果为$x^5+x^3+x^2+1=(x^4+x+1)\cdot x+x^3+x+1 \equiv x^3+x+1\ mod \ f(x)$

同样的，这个例子也可以用系数向量来表示，同样的系数向量中多项式的系数以$x$的阶数升序表示，因此上述乘法可以表示为$[1,0,1,0]\cdot [1,0,0,1]=[1,1,0,1]$

采用二元域的原因有两个

• 编码二元数据很有用，n bit的数据就是n bit长，素数域通常无法完成，因为素数除了2以外都不是2的次幂

如果是模65537，可以每2 Bytes编码为一个模65537的有限域内的元素，但是如果采用FFT的话，结果可能会输出65536，这个数无法用2 Bytes表示

• 二元域还有一个好处就是加法等于减法，且不包含任何偶数项（模2下偶数项均为0），奇数项的系数一定是1（奇数模2一定是1）

如果说需要利用二元域进行编码，啧二元域上的FFT一定得看看，但是二元域上的FFT有个问题，二元域没有阶为$2^n$的非平凡乘法群，因为乘法群的阶均为$2^n-1$

举个例子，如果说在模$f(x)=x^4+x+1$条件下计算$x+1$的幂次，则会在15次幂得到1，而非16次，这是因为有限域内的元素个数为16，且其中有一个为0，对于一个非零的数进行幂运算是无法得到零的（前提是模的是不可约多项式），因此$x+1$的幂次会在15次开始循环

还记得我们之前提到了对子集的结构有特殊的要求，我们要求子集必须是大小为2的整数次幂的乘法群，因为这样每次平方都可以使得大小减半，但是二元域会如何呢

给定一个包含$2^k$个元素的集合$D$，其中包括零元，此时我们可以构造一个只有原集合一半大小的集合$D'$，只需要选择一个特定的$k$，对原本$D$中的$x$，乘一个$x+k$即可

因为原来的$D$是一个子空间，且$k$也在$D$中，因此$x+k$也在$D$中，从而我们得到了一个二对一的映射

举个例子，假设$k=1$，我们可以得到这个列表

如果我们希望用FFT，和之前一样，将一个大小为$n$的多项式和大小为$n$的集合划分为两个大小为$\frac{n}{2}$的多项式和集合，不过和之前在素数域不同的是，我们采用不同的等式

这里我们将多项式$P$划分为$P_{odd},P_{even}$，满足

$$ P(x)=P_{odd}(x\cdot (k-x))+x\cdot P_{odd}(x\cdot (k-x)) $$

对于我们找到的$P_{odd},P_{even}$，其同样满足

$$ P(x+k)=P_{odd}(x\cdot (k-x))+(x+k)\cdot P_{odd}(x\cdot (k-x)) $$

因此我们可以递归的采用FFT来缩小集合$D$

将$P$转换为$P_{odd},P_{even}$会导致其本身是非平凡的，这样一来最笨的办法的期望为$O(n^2)$，在二元域上我们可以利用一个事实$(x^2-kx)^2=x^4-k^2x^2$，或者更一般的$(x^2-kx)^{2^i}=x^{2^{i+1}}-k^{2^i}x^{2^i}$，这样一来期望复杂度就降低至了$O(nlog n)$

如果需要逆向FFT来做插值法，需要以相反的顺序运行算法中的步骤，具体可以看$[2]$或者针对算法优化的论文$[3]$

至此，我们得到了一种高效、可扩展的方法来计算和插值多项式，如果我们想使用快速傅立叶变换来恢复缺失部分的编码数据，那么也存在相应的算法，但是这些算法会比FFT效率低一些

References

$[1]$ Fast Fourier Transforms

$[2]$ Github - Binary_fft

$[3]$ Fast Reed-Solomon Interactive Oracle Proofs of Proximity