FFT & NTT

$e^{i \theta}=\cos (\theta)+i * \sin (\theta)$

FFT-Fast Fourier Transformation

快速傅里叶变换(Fast Fourier Transform, FFT),是快速计算序列的离散傅立叶变换(DFT)及其逆变换的方法,将计算复杂度从 $O(n^2)$ 降低到 $O(n\log n)$ ,其中 $n$ 为序列的长度。在密码学场景下,就是为了加速多项式乘法。

首先,我们考虑多项式:

我们要计算 $C(x) = A(x) \cdot B(x)$ ,即两个多项式的乘积,那么我们可以想到的最简单的方法就是$O(N^2)$:将每一项依次分别相乘再相加。

这样的时间复杂度肯定是不可接受的,所以我们需要寻找更快的方法。

Priliminary

单位根

首先,我们需要引入一个复数 $w$ ,满足 $w^n = 1$ ,即 $w$ 是 $n$ 次单位根。考虑欧拉公式$e^{i \theta}=\cos (\theta)+i * \sin (\theta)$,我们可以得到 $w_n^k = e^{2\pi i k / n}, k\in [n]$ ,可以发现,复数的 $n$ 次单位根平分单位圆。
单位根还具有一些特殊的性质:

  1. $w_{dn}^{dk} = \cos (2\pi d k / dn)+i * \sin (2\pi d k / dn) = w_n^k$
  2. $w_n^{k+n/2} = \cos (2\pi k / n + \pi)+i * \sin (2\pi k / n + \pi) = - \cos (2\pi k / n ) - i \sin (2\pi k / n + \pi) = -w_n^{k}$
  3. $w_n^{k+n} = \cos (2\pi k / n + 2\pi)+i * \sin (2\pi k / n + 2\pi) = w_n^{k}$
  4. $w_n^0 = 1$

DFT

已知 $A(x)$ 的系数为 $\left(a_0, a_1, \ldots, a_{n-1}\right)$, 对于 $k=0,1, \ldots, n-1$, 定义:

其中向量 $y=(y_0, y_1, \ldots, y_{n-1})$ 是系数向量 $a=\left(a_0, a_1, \ldots, a_{n-1}\right)$ 的离散傅立叶变换。

多项式系数表示法

设 $A(x)$ 表示一个d次多项式, 则 $A(x)=a_0+a_1x+\ldots,+a_d x^d$
利用这种方法计算多项式卷积复杂度为 $O\left(d^2\right)$, 其实就是直接对应相乘(暴力)。
例如: $A(x)=1+2 x+x^2, B(x)=1-2 x+x^2$

多项式点值表示法

设 $A(x)$ 表示一个d次多项式, 我们知道该多项式在 $d+1$ 个点的值就可以唯一确定这个多项式。
因此,我们任取 $d+1$ 个点 $x_0,x_1,\ldots,x_d$,则 $A(x)$ 在这些点的值为 $A(x_0),A(x_1),\ldots,A(x_d)$,我们可以用这些点的值来表示多项式 $A(x)$。

在点值表示法下,多项式的乘法可以转化为点值的乘法,即 $A(x)B(x)$ 在 $d+1$ 个点的值为 $A(x_i)B(x_i)$ 的乘积,计算复杂度是 $O(d)$。但是要注意一个问题$A(x)B(x)$是一个 $2d$ 次多项式,因此我们需要找到 $2d+1$ 个点来表示这个多项式。

因此,我们想到了一个方法,将多项式转化为点值表示,进行$O(d)$的乘法,然后再将点值表示的结果转化为多项式。但是这个转化的过程是如何进行的呢?如果我们采用普通求值的方法,复杂度仍然是 $O(d^2)$ 的。接下来,我们就开始逐步进行优化,引入FFT算法。

FFT Algorithm

我们现在的首要问题是如何将多项式快速转化为点值表示,并且如何将点值表示的结果句快速转化为多项式,这就是单位根的作用了。我们观察我们手中的多项式 $F(x)$ ,假设$n=2^k$,我们可以将其分为奇数项和偶数项:

将其转化为两个$n/2$次的多项式$FL(x)$和$FR(x)$

如果$k< n/2$,把$w_n^k$代入$F(x)$,我们可以得到:$F(w_n^k)=FL(w_{n}^{2k})+w_n^k FR(w_{n}^{2k})$,即$F(w_n^k)=FL(w_{n/2}^{k})+w_n^k FR(w_{n/2}^{k})$;如果$k> n/2$,则有$F(w_n^k)=FL(w_{n/2}^{k-n/2})- w_n^{k-n/2} FR(w_{n/2}^{k-n/2})$。
如果我们知道了$FL(x)$和$FR(x)$在$w_{n/2}^{k}$处的值,我们就可以通过上述公式计算出$F(x)$在$w_n^k$处的值。
接着我们就可以利用分治思想,通过一个递归的算法来计算出$F(x)$ 在 $w_n^0, w_n^1, \ldots, w_n^{n-1}$ 处的点值表示。复杂度为$O(n\log n)$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def FFT(P):
n = len(P)
if n == 1:
return P
w = 1
wn = cmath.exp(2 * cmath.pi * 1j / n)
P0 = [P[i] for i in range(0, n, 2)]
P1 = [P[i] for i in range(1, n, 2)]
y0 = FFT(P0)
y1 = FFT(P1)
y = [0] * n
for k in range(n // 2):
y[k] = y0[k] + w * y1[k]
y[k + n // 2] = y0[k] - w * y1[k]
w *= wn
return y

IFFT Algorithm

对于点值计算,其实就是计算一个矩阵乘法:

中间的范德蒙矩阵就是一个DFT矩阵,有了正向的DFT,我们就可以通过逆矩阵来计算逆DFT,即IFFT,即

这样我们就可以通过IFFT来将点值表示的多项式转化为系数表示的多项式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def IFFT(P):
n = len(P)
if n == 1:
return P
w = 1
wn = cmath.exp(-2 * cmath.pi * 1j / n)
P0 = [P[i] for i in range(0, n, 2)]
P1 = [P[i] for i in range(1, n, 2)]
y0 = IFFT(P0)
y1 = IFFT(P1)
y = [0] * n
for k in range(n // 2):
y[k] = y0[k] + w * y1[k]
y[k + n // 2] = y0[k] - w * y1[k]
w *= wn
return y

蝴蝶变换

蝴蝶变换是一个很简单的东西,其实就是我们分治的顺序转换为二进制之后,发现每一个数都是其二进制反过来。这样我们就可以通过一个循环来计算出所有的点值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def rev(i):
return int(bin(i)[:1:-1], 2)

def FFT(P,IFFT):
n = len(P)
for i in range(n):
j = rev(i)
if i < j:
P[i], P[j] = P[j], P[i]
l = 2
while l <= n:
wn = cmath.exp(2 * cmath.pi * 1j / l)
if IFFT:
wn = cmath.exp(-2 * cmath.pi * 1j / l)
for i in range(0, n, l):
w = 1
k = i
while k < i + l // 2:
u = P[k]
v = w * P[k + l // 2]
P[k] = u + v
P[k + l // 2] = u - v
w *= wn
k += 1
l *= 2

NTT-Number Theoretic Transform

NTT是FFT的一种变种,其实质是在模数为 $p$ 的情况下,将FFT的复数域转化为整数域。就是将单位根 $w_n$ 替换为 $g$,其中 $g$ 是模数 $p$ 的原根。这样我们就可以在整数域上进行FFT的计算,此处不再赘述。

Reference

[1] https://www.cnblogs.com/Ning-H/p/12072626.html
[2] https://www.cnblogs.com/Ning-H/p/13693580.html
[3] https://www.cnblogs.com/pam-sh/p/15976275.html