快速傅立叶变换FFT与快速数论变换NTT

FFT与NTT

在详细讲解FFT的内容之前，我们先主要回顾一下关于多项式和原根的相关内容：
多项式
系数表达式：若一个n次多项式$A(x)=a_0+a_{1}x+\dots +a_n{x}^n$是由其系数$(a_0,a_1,\dots ,a_n)$唯一确定，则这种方式被称为多项式的系数表达形式。有时也被写为如下形式：
$$A(x)=\sum _{i=0}^n{a}_i{x}^i$$
点值表达式：若将多项式定义域内的任意一个值$x_0$代入上面的系数表达式，则可以得到多项式的值$y_0$。根据多项式插值的思想，一个多项式能够被至少n+1个不同的点(xi,yi)，i∈{0,n}(x_i,y_i)，i\in \{0,n\}(xi​,yi​)，i∈{0,n}所唯一确定,这种方式被称为多项式的点值表达式。

单位根与原根
单位根
若想将多项式的系数表达式转化为点值表达式，我们需要取n+1个不同的点来进行计算，由于复平面上的1的原根具有很好的性质，且可以满足我们的要求，故将其作为介绍目标。形式化定义为：$w^n=1$，记为$w_n$。对于任意一个正整数n，有n个n次单位根恰好均匀分布在复平面单位圆上。

有欧拉恒等式可得$e^{\pi i}+1=0$，即$e^{\pi i}=-1$。则我们可以得到$e^{2\pi i}=1$。故我们令
$$w_n^k=e^{\frac{2k\pi i}{n}}=(e^{\frac{2\pi i}{n}}{)}^k$$
具体例子如下图：

且关于单位根有下面三个重要的定理：
[1]$w_n^{k+\frac{n}{2}}=e^{\frac{2\pi i}{n}(k+\frac{n}{2})}=e^{\frac{2\pi i}{n}k}{e}^{\pi i}=-e^{\frac{2\pi i}{n}k}=-w_n^k$
[2]$w_{2n}^{2k}=w_n^k$
[3]$w_n^{-k}=1e^{\frac{2\pi i}{n}(-k)}=e^{2\pi i}{e}^{\frac{2\pi i}{n}(-k)}=e^{\frac{2\pi i}{n}(n-k)}=w_n^{n-k}$

原根
【阶】若$r,n$是互素的整数，且$r\ne 0,n>0$,使得等式$r^x\equiv 1(modn)$成立的最小正整数$x$称为$r$模$n$的阶，通常记为$Or{d}_n(r)$。
【欧拉函数】假设$n$是一个正整数，则它的欧拉函数就是小于$n$且与$n$互素正整数的个数，记为$\varphi (n)$。
【原根】若$r,n$是互素的整数，且$r\ne 0,n>0$和$Or{d}_n(r)=\varphi (n)$,则称$r$是模$n$的原根。

FFT
为了方便理解，假设我们的计算目标是C(x)=A(x)B(x),且多项式A(x)和B(x)分别是a,b。令n=a+b。
核心思想：利用拉格朗日插值的思想，实现多项式的系数表达式和点值表达式之间的转换。

主要工具：
离散傅立叶变换：(DFT)
$$X_k=\sum _{j=0}^{n-1}{x}_j{e}^{\frac{2\pi i}{n}kj}$$
离散傅立叶逆变换：(IDFT)
$$x_j=\frac{1}{n}\sum _{k=0}^{n-1}{X}_k{e}^{-\frac{2\pi i}{n}jk}$$

具体步骤：
[1]：取n个单位根：$w_n^0$，$w_n^1$，…，$w_n^{n-1}$
[2]：将n个不同的单位根代入多项式$A(x)$和$B(x)$,即计算$A(w_n^k)$和$B(w_n^k)$
[3]：计算$C(w_n^k)=A(w_n^k)B(w_n^k)$
[4]：令$D(x)=C_0+C_1{x}^1+\dots +C_{n-1}{x}^{n-1}$，并计算$D(w_n^k)$
[5]：计算多项式$C(x)$的系数$c_j=\frac{1}{n}D(w_n^{n-j})$

详细过程：由于第[1]、[3]步比较简单，我们直接跳过。
【2】计算$A(w_n^k)$、$B(w_n^k)$
最直接将n个单位根代入多项式，但是这样的时间复杂度比较高。故这里我们采用分治的思想来进行加速。举个例子：
假设有一个多项式$H(x)=h_0+h_1{x}^1+h_2{x}^2+h_3{x}^3$，将该多项式写成按奇偶项拆分成两个子多项式的和，则有：
$$H(x)=(h_0+h_2{x}^2)+(h_1{x}^1+h_3{x}^3)=(h_0+h_2{x}^2)+x(h_1+h_3{x}^2)$$
$$H(x)=H_1(x^2)+x{H}_2(x^2)$$
若将这种思想用在计算$A(w_n^k)$、$B(w_n^k)$上，我们以计算$A(w_n^k)$为例,$0\le k<\frac{n}{2}$，有：
当单位值在复平面单位圆的上部分时：
$$A(w_n^k)=A_1(w_n^{2k})+w_n^k{A}_2(w_n^{2k})=A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_2(w_{\frac{n}{2}}^k)$$
当单位值在复平面单位圆的下部分时：
$$A(w_n^{k+\frac{n}{2}})=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}{A}_2(w_n^{2k+n})=A_1(w_n^{2k})-w_n^k{A}_2(w_n^{2k})==A_1(w_{\frac{n}{2}}^k)-w_n^k{A}_2(w_{\frac{n}{2}}^k)$$
由此可见，$k$只需取一半范围的值就可以得到整个范围的计算结果。事实上，多项式$A_1$和$A_2$的计算可以利用分治的思想，这里明白即可，不再赘述。

【4】在计算$C(w_n^k)=A(w_n^k)B(w_n^k)$之后，后面的主要工作是利用IDFT将$(w_n^k,C(w_n^k))$转化为系数表达式，即求$c_j$。为了方便描述和形式化表达，我们这里引入一个n-1次多项式$D(x)=C_0+C_1{x}^1+\dots +C_{n-1}{x}^{n-1}={\sum }_{k=0}^{n-1}{C}_k{x}^k$,且$C_k=C(w_n^k)$已经在【3】步中计算完毕。
【5】由IDFT可得：
$$c_j=\frac{1}{n}\sum _{k=0}^{n-1}{C}_k{e}^{-\frac{2\pi i}{n}jk}=\frac{1}{n}\sum _{k=0}^{n-1}{C}_k(w_n^{-j}{)}^k=\frac{1}{n}\sum _{k=0}^{n-1}{C}_k(w_n^{n-j}{)}^k=\frac{1}{n}D(w_n^{n-j})$$
若将上述的过程中的$D(x)$看作一个系数已知的多项式，其求$c_j$的过程也可以使用分治的思想进行加速。

注意：关于DFT和IDFT的证明可以直接在网上搜索，比较简单，不进行叙述。

NTT
由于在密码学中，多项式的卷积并不是像FFT中的一样，它们往往是循环多项式的卷积，即在多项式的卷积后要模上一个大素数。故我们需要对FFT进行相应的改造，即NTT。

需要说明的是，在FFT中我们选取的是复平面的单位根，而在NTT中我们使用的是相应的原根。这是因为假设当$p$是一个奇素数，$g$是模$n$的原根，当$0<i,j<p$,有$g^i\equiv g^j(modp)$,当且仅当$i=j$。该性质与单位根十分类似，但由于在NTT中我们所取的不同根的个数是与多项式的次数是相关的，而上面$i,j$的最大值只能取到$p-1$,即只能取$p$个不同的值，故我们要对原根继续进行改造：

给定模数$M=c\ast 2^k+1$(这是满足密码学中多项式循环卷积的模数形式的)及其原根$g$以及正整数$n=2^l$(通常是多项式的次数)，且$l\le k$，则有：
$$M-1=c\ast 2^k$$
$$\frac{M-1}{n}=\frac{c\ast 2^k}{2^l}=c\ast 2^{k-l}$$
若令$a=g^{\frac{M-1}{n}}$，则刚好能够构造出$n$个$modM$不同值,满足$a^i\equiv (g^{\frac{M-1}{n}}{)}^i\equiv (g^{M-1}{)}^{\frac{i}{n}}\equiv 1(modM)$。
下面是原根满足的两个性质:
[1]：$a^n\equiv 1(modM)$
[2]：$a^{\frac{n}{2}}\equiv -1(modM)$

主要工具：对于次数为$N=2^l$(不足2的整数次幂的可以通过补0来扩张)且系数元素均小于$M$的多项式系数序列$x(n)$,$n$表示第n个系数,有相应的快速数论变换：
快速数论变换：(NTT)
$$X(m)\equiv \sum _{n=0}^{N-1}x(n)a^{mn}modM$$
快速数论逆变换：(INTT)
$$x(n)=\frac{1}{N}\sum _{m=0}^{N-1}X(m)a^{-mn}modM$$

蝶形操作：在FFT中，将多项式的系数表达式转化为点值表达式的过程中，我们主要采取分治的思想进行加速。而在NTT中，我们采用的是迭代算法——碟形操作。我们采用数学归纳法来进行说明，具体过程如下：
(1)当$N=2^2$时，则多项式系数为$x(0),x(1),x(2),x(3)$,同时$m=0,1,2,3$。NTT可以写为下面的形式：
$$X(m)\equiv \sum _{n=0}^{3}x(n)a^{mn}modM,m=0,1,2,3$$
而将$m,n$写成二进制表达式为：
$$m=2m_1+m_0=(m_1,m_0{)}_2$$
$$n=2n_1+n_0=(n_1,n_0{)}_2$$
所以NTT又可改写如下：

上式成立，因为$a^N\equiv 1(modM)$，且$N=2^2$,则
$$a^{2^2{m}_1{n}_1}\equiv 1(modM)$$
然后，令$x(n_1,n_0)=x_0(n_1,n_0)$,且
$$x_1(m_0,n_0)=\sum _{n_1=0}^1{x}_0(n_1,n_0)a^{2n_1{m}_0}=x_0(0,n_0)+x_0(1,n_0)a^{2m_0}(modM)$$
则，遍历所有$m_0,n_0$的值，我们可得到：

相应的蝶形操作如为：

将$x_1(m_0,n_0)$代入到原式中有：
$$x_2(m_0,m_1)=\sum _{n_0=0}^1{x}_1(m_0,n_0)a^{(2m_1+m_0)n_0}=x_1(m_0,0)+x_1(m_0,1)a^{2m_1+m_0}(modM)$$
进一步有：
$$X(m_1,m_0)=x_2(m_0,m_1)(modM)$$
相应的蝶形操作如下图：

(2)当$N=2^3$时，同理可得：

将红色部分提取出来有：
$$x_1(m_0,n_1,n_0)=x_0(0,n_1,n_0)+x_0(1,n_1,n_0)a^{4m_0}(modM)$$
遍历$m_0,n_1,n_0$的值，我们可以得到下图的蝶形操作图：

将$x_1(m_0,n_1,n_0)$代入到原式可得：
$$X(m_2,m_1,m_0)=\sum _{n_0=0}^1(\sum _{n_1=0}^1{x}_1(m_0,n_1,n_0)a^{(2^2{m}_1+m_0)n_1})a^{(2^2{m}_2+2m_1+m_0)n_0}(modM)$$
将括号内的式子提取出来，并令：
$$x_2(m_0,m_1,n_0)=\sum _{n_1=0}^1{x}_1(m_0,n_1,n_0)a^{(2^2{m}_1+m_0)n_1}=x_1(m_0,0,n_0)+x_1(m_0,1,n_0)a^{4m_1+2m_0}(modM)$$
遍历$m_0,m_1,n_0$的值，我们可以得到下图的蝶形操作图：

再将$x_2(m_0,m_1,n_0)$代回到原式，我们可以得到：
$$x_3(m_0,m_1,m_2)=\sum _{n_0=0}^1{x}_2(m_0,m_1,n_0)a^{(2^2{m}_2+2m_1+m_0)n_0}=x_2(m_0,m_1,0)+x_2(m_0,m_1,1)a^{4m_2+2m_1+m_0}(modM)$$
同样遍历$m_0,m_1,m_2$的值，我们可以得到下图的蝶形操作图：

(3)根据数学归纳法，我们可以得到$N=2^l$的情况：
NTT中的正整数序列$m,n$可以写为二进制形式：


则$a^{mn}$可以计算为：

进而NTT可以写为下面的形式：

利用(1)、(2)中的思想，将红色内容提取出来有：

根据数学归纳法，我们可以得到$x_r(m_0,\dots ,m_{r-1},n_{l-r-1},\dots ,n_0)$的一般递推式：

最后，就可以得到：

上面内容就是NTT的一般化流程，而每一次红色部分内容的提出本质上就是蝶形操作的一轮。下面就以$r$轮为例，详细说明蝶形操作的具体流程，即上面的公式(2.7)。在上面的公式中，我们可以容易地发现$x_r$和$x_{r-1}$中除了$n_{l-r}$不同外，其他部分都一样，故为了方面描述，我们可以简写为：

公式(2.7)可以写为：

因为$m_{r-1}$和$n_{l-r}$的取值为0和1，故对应的蝶形操作为：

其中，关于$p$有下面的结论：

证明过程如下，已知$a^{2^{l-1}}=-a^0$：
[1]当$r=1$时

[2]当$r>1$时

注意：以上部分就是NTT具体数学推导流程，而INTT的思想与NTT完全一致，除了需要在最后的结果上乘上$\frac{1}{N}=\frac{1}{2^l}$。而关于NTT在实现方面是有许多改进和注意的地方，这是后续的工作内容，这里就先不叙述了。