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$ 为序列的长度。在密码学场景下，就是为了加速多项式乘法。

首先，我们考虑多项式:

$$A(x) = a_n x^n + a_{n-1} x^{n-1} + a_{n-2} x^{n-2} + \cdots + a_1 x + a_0 B(x) = b_n x^n + b_{n-1} x^{n-1} + b_{n-2} x^{n-2} + \cdots + b_1 x + b_0$$

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

$$C(x) = \sum_{i=0}^{2n} c_i x^i = \sum_{i=0}^{2n} \left( \sum_{j=0}^{i} a_j b_{i-j} \right) x^i$$

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

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_k=A\left(w_n^k\right)=\sum_{i=0}^{n-1} a_i w_n^{k i}$$

其中向量 $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) B(x)=\left(1+2 x+x^2\right)\left(1-2 x+x^2\right)=1-2 x^2+x^4$$

多项式点值表示法

设 $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)$。

$$T = ((x_0,A(x_0)),(x_1,A(x_1)),\ldots,(x_d,A(x_d)))$$

在点值表示法下，多项式的乘法可以转化为点值的乘法，即 $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$，我们可以将其分为奇数项和偶数项：

$$F(x)=a_0+a_1 x+a_2 x^2+\ldots+a_{2n-1} x^{2n-1}=a_0+a_2 x^2+\ldots+a_{n-2} x^{n-2}+x(a_1+a_3 x+\ldots+a_{n-1} x^{n-2})$$

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

$$FL(x)=a_0+a_2 x+\ldots+a_{n-2} x^{n/2-1}\\ FR(x)=a_1+a_3 x+\ldots+a_{n-1} x^{n/2-1}\\ F(x)=FL(x^2)+x FR(x^2)$$

如果$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)$。

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

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

$$\begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 & \cdots & 1\\ 1 & w_n & \cdots & w_n^{n-1}\\ \vdots & \vdots & \ddots & \vdots\\ w_n^{n-1} & w_n^{2(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix}$$

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

$$\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} 1 & 1 & \cdots & 1\\ 1 & w_n^{-1} & \cdots & w_n^{-(n-1)}\\ \vdots & \vdots & \ddots & \vdots\\ w_n^{-(n-1)} & w_n^{-2(n-1)} & \cdots & w_n^{-(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{n-1} \end{bmatrix}$$

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

def IFFT_in(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

def IFFT(P):
    n = len(P)
    y = IFFT_in(P)
    for i in range(n):
        y[i] /= n
    return y

蝴蝶变换

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

def rev(i, n):
    result = 0
    for j in range(n.bit_length() - 1):
        result = (result << 1) | (i & 1)
        i >>= 1
    return result

def FFT(P, IFFT=False):
    n = len(P)
    log_n = n.bit_length() - 1
    
    # Reorder P based on reversed bit indices
    for i in range(n):
        j = rev(i, n)
        if i < j:
            P[i], P[j] = P[j], P[i]
    
    # FFT or IFFT computation
    l = 2
    while l <= n:
        angle = (2 * cmath.pi / l) * (-1 if IFFT else 1)
        wn = cmath.exp(angle * 1j)
        
        for i in range(0, n, l):
            w = 1
            for k in range(i, i + l // 2):
                u = P[k]
                v = w * P[k + l // 2]
                P[k] = u + v
                P[k + l // 2] = u - v
                w *= wn
        l *= 2
    
    # Normalize if it's IFFT
    if IFFT:
        for i in range(n):
            P[i] /= n

    return P

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