【NTT(Number Theoretic Transforms)（一）】

提示：Kyber中的NTT(Number Theoretic Transforms)

NTT(Number Theoretic Transforms)（一）

前言

数论变换（number-theoretic transform, NTT）是离散傅里叶变换（DFT）在数论基础上的实现；快速数论变换（fast number-theoretic transform, FNTT）是 快速傅里叶变换（FFT）在数论基础上的实现。针对具有代数结构的格密码算法，其中的多项式乘法为主要的核心算子，目前多数具有代数结构的格密码算法均通过NTT进行多项式乘法加速。本文以NIST公布的标准算法Kyber为例，以此介绍NTT算法的实现原理。

一、中国剩余定理（CRT）

Kyber中的多项式乘法均取自多项式环$R_q={\mathbb{Z}}_q/(x^n+1)$，其中第三轮的Kyber中$q=3329=2^8+1$，$n=256$。由于$256||(q-1)$，所以有限域${\mathbb{Z}}_q$中存在256次本原单位根，但不存在512次本原单位根，因此$(x^n+1)modq$可分解为128个二次多项式。令$\zeta =17$为该256次本原单位根，因此$\zeta ,{\zeta }^3,{\zeta }^5,\cdots ,{\zeta }^{255}$便为${\mathbb{Z}}_q$中所有的256次单位根，固有如下式子：
$$x^n+1=(x^2-\zeta )(x^2-{\zeta }^3)\cdots (x^2-{\zeta }^{255})modq$$于是由中国剩余定理有：$${\mathbb{Z}}_q/(x^{256}+1)={\mathbb{Z}}_q/\prod _{i=0}^{127}(x^2-{\zeta }^{2i+1})\cong {\mathbb{Z}}_q/(x^2-\zeta )⨂\cdots ⨂Z_q/(x^2-{\zeta }^{255})$$ 具体的分解过程可通过平方差公式得到更详细的分解过程，首先${\mathbb{Z}}_q/(x^{256}+1)$，由于${\zeta }^{256}=1$，所以${\zeta }^{128}=-1$，于是有${\mathbb{Z}}_q/(x^{256}+1)={\mathbb{Z}}_q/(x^{256}-{\zeta }^{128})$，由平方差公式便有
$${\mathbb{Z}}_q/(x^{256}-{\zeta }^{128})\cong {\mathbb{Z}}_q/(x^2-{\zeta }^{64})\otimes {\mathbb{Z}}_q/(x^2-{\zeta }^{192})$$其中${\mathbb{Z}}_q/(x^2-{\zeta }^{192})$是由$x^2+{\zeta }^{64}=x^2-{\zeta }^{128}\cdot {\zeta }^{64}$得来。由此继续向下利用平方差公式进行分解，直到最后的
$${\mathbb{Z}}_q/(x^2-{\zeta }^{2brv(0)+1})\otimes {\mathbb{Z}}_q/(x^2-{\zeta }^{2brv(1)+1})\otimes \cdots \otimes {\mathbb{Z}}_q/(x^2-{\zeta }^{2brv(127)+1})$$ 便不可再分。其中$brv(i)$表示为$i$的$\log n-1$位二进制表示比特翻转，即Kyber文档中的$b{r}_7(i)$。比特翻转的顺序结果是由于每次对第二个加的式子使用平方差公式时均需要乘上${\zeta }^{128}$使其变成满足平方差公式的形式。详见下图（图片来源于知乎）。
故Kyber中多项式可表示为：
$$(fmod(x^2-\zeta ),fmod(x^2-{\zeta }^{2brv(1)+1}),\cdots ,fmod(x^2-{\zeta }^{2brv(127)+1}))$$
所以f(x)的NTT表示为：
$$NTT(f)=\hat{f}(x)={\hat{f}}_0+{\hat{f}}_{1}x+\cdots +{\hat{f}}_{255}{x}^{255}$$
其中：
$${\hat{f}}_{2i}=\sum _{j=0}^{127}{f}_{2j}\cdot {\zeta }^{(2brv(i)+1)j};$$ $${\hat{f}}_{2i+1}=\sum _{j=0}^{127}{f}_{2j+1}\cdot {\zeta }^{(2brv(i)+1)j}.$$也可表示为：
$$NTT(f)=\hat{f}(x)=({\hat{f}}_0+{\hat{f}}_{1}x,{\hat{f}}_2+{\hat{f}}_{3}x,\cdots ,{\hat{f}}_{254}+{\hat{f}}_{255}x).$$
于是$h(x)=f(x)\cdot g(x)$利用NTT计算便为：
$$\hat{h}(x)=\hat{f}(x)\cdot \hat{g}(x)=({\hat{h}}_0+{\hat{h}}_{1}x,{\hat{h}}_2+{\hat{h}}_{3}x,\cdots ,{\hat{h}}_{254}+{\hat{h}}_{255}x).$$其中：
$${\hat{h}}_{2i}+{\hat{h}}_{2i+1}x=({\hat{f}}_{2i}+{\hat{f}}_{2i+1}x)\cdot ({\hat{g}}_{2i}+{\hat{g}}_{2i+1}x)mod{x}^2-{\zeta }^{2i+1}$$故
$${\hat{h}}_{2i}={\hat{f}}_{2i}\cdot {\hat{g}}_{2i}+{\hat{f}}_{2i+1}\cdot {\hat{g}}_{2i+1}\cdot {\zeta }^{2i+1};$$ $${\hat{h}}_{2i+1}={\hat{f}}_{2i}\cdot {\hat{g}}_{2i+1}+{\hat{f}}_{2i+1}\cdot {\hat{g}}_{2i}.$$那么如何快速的计算${\hat{f}}_{2i}$和${\hat{f}}_{2i+1}$以及其逆变换呢？这可以从快速傅里叶变换（FFT）中得到启发。

二、离散傅里叶变换（DFT）

$n-1$次多项式有两种表述方式，一种便是由$n$个系数直接写出，另一种便是利用插值多项式原理，$n$个不同变量的函数值唯一决定了一个多项式。即$f(x)$可以表示为$f_0+f_{1}x+\cdots +f_{n-1}{x}^{n-1}$，也可表述为满足$n$个函数值对$(x_0,f(x_0)),\cdots ,(x_{n-1},f(x_{n-1}))$的$n-1$多项式，其中$x_i\ne x_j$,$i\ne j$。通过这两种等价表述，$h(x)=f(x)\cdot g(x)$可通过计算$(x_0,f(x_0)\cdot g(x_0)),\cdots ,(x_{(}n-1),f(x_{n-1})\cdot g(x_{n-1}))$得到。于是可寻找$n$个不同的自变量$x_0,x_1,\cdots ,x_{n-1}$，分别计算$f(x_i)$和$g(x_i)$，再分别进行相乘即可得到$n$个$h(x)$的不同点对，再利用插值多项式的计算方法便可得到最终的乘积。
在离散傅里叶变换（DFT）中，上述$n$个不同的自变量取$n$次单位根即可，即$w_n^k=e^{2\pi i\cdot k/n}$，其中$k=0,1,2,\cdots ,n-1$；NTT便取有限域中的$n$个本原单位根。于是离散傅里叶变换即为计算$f(e^{2\pi i\cdot k/n})$，所以$f(x)$的离散傅里叶变换公式为：
$$F_k=\sum _{j=0}^{n-1}{f}_j\cdot (e^{2\pi i\cdot k/n}{)}^j=\sum _{j=0}^{n-1}{f}_j\cdot e^{2kj\pi i/n}$$.逆变换便是由$(F_0,F_1,\cdots ,F_{n-1})$计算$(f_0,f_1,\cdots ,f_{n-1})$，上述式子利用矩阵的形式表示为：
$$\left(\begin{array}{c}{F}_0\\ F_1\\ ⋮\\ F_{n-1}\end{array}\right)=\left(\begin{array}{cccc}1& 1& \cdots & 1\\ 1& e^{2\pi i\cdot 1/n}& \cdots & (e^{2\pi i\cdot 1/n}{)}^{n-1}\\ ⋮& ⋮& ⋮& ⋮\\ 1& e^{2\pi i\cdot (n-1)/n}& \cdots & (e^{2\pi i\cdot (n-1)/n}{)}^{n-1}\end{array}\right)\cdot \left(\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_{n-1}\end{array}\right)$$ 所以有：$$\left(\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_{n-1}\end{array}\right)={\left(\begin{array}{cccc}1& 1& \cdots & 1\\ 1& e^{2\pi i\cdot 1/n}& \cdots & (e^{2\pi i\cdot 1/n}{)}^{n-1}\\ ⋮& ⋮& ⋮& ⋮\\ 1& e^{2\pi i\cdot (n-1)/n}& \cdots & (e^{2\pi i\cdot (n-1)/n}{)}^{n-1}\end{array}\right)}^{-1}\cdot \left(\begin{array}{c}{F}_0\\ F_1\\ ⋮\\ F_{n-1}\end{array}\right)$$ 因为由单位根构成的系数矩阵为一个范德蒙德矩阵，所以该矩阵的逆容易求出，即：$$\left(\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_{n-1}\end{array}\right)=\frac{1}{n}\left(\begin{array}{cccc}1& 1& \cdots & 1\\ 1& e^{-2\pi i\cdot 1/n}& \cdots & (e^{-2\pi i\cdot 1/n}{)}^{n-1}\\ ⋮& ⋮& ⋮& ⋮\\ 1& e^{-2\pi i\cdot (n-1)/n}& \cdots & (e^{-2\pi i\cdot (n-1)/n}{)}^{n-1}\end{array}\right)\cdot \left(\begin{array}{c}{F}_0\\ F_1\\ ⋮\\ F_{n-1}\end{array}\right)$$ 故离散傅里叶逆变换为：$$f_j=\frac{1}{n}\sum _{k=0}^{n-1}{F}_k\cdot (e^{-2\pi i\cdot j/n}{)}^k=\frac{1}{n}\sum _{k=0}^{n-1}{F}_k\cdot e^{-2kj\pi i/n}$$ 类似地，也可得到NTT及其INTT的计算公式。此时可以分为两种情形，${\mathbb{Z}}_q$中存在$2n$次本原单位根（Dilithium）和只存在$n$次本原单位根（Kyber）。一般讨论n为2的方幂的情形，如不满足则对多项式的高次系数填充0即可。
对于第一种情形，由于${\mathbb{Z}}_q$中存在$2n$次本原单位根${\zeta }_{2n}$，于是$x^n+1$可以完全分解为$n$个1次多项式，$({\zeta }_{2n},{\zeta }_{2n}^3,\cdots ,{\zeta }_{2n}^{2n-1})\in {\mathbb{Z}}_q$即为$n$个$n$次本原单位根，因此可以通过上述离散傅里叶变换的方式进行NTT变换，只需将其中的$n$次单位根变为本原单位根即可。即$${\hat{f}}_k=\sum _{k=0}^{n-1}{f}_j\cdot ({\zeta }_{2n}^{2k+1}{)}^j=\sum _{k=0}^{n-1}{f}_j\cdot {\zeta }_{2n}^{(2k+1)j};$$ $$f_j=\frac{1}{n}\sum _{k=0}^{n-1}{\hat{f}}_k\cdot ({\zeta }_{2n}^{-(2j+1)}{)}^k=\frac{1}{n}\sum _{k=0}^{n-1}{\hat{f}}_k\cdot {\zeta }_{2n}^{-(2j+1)k}.$$ 也可通过第一章中中国剩余定理的原理得到。但与第一章中不同之处在于此时可使用平方差公式分解至一次多项式，即：$${\mathbb{Z}}_q/(x^n+1)\cong {\mathbb{Z}}_q/(x-\zeta )\otimes {\mathbb{Z}}_q/(x-{\zeta }^3)\cdots \otimes {\mathbb{Z}}_q/(x-{\zeta }^{2n-1}).$$ 将${\zeta }^{2k+1}$分别代入$f(x)$便得到了在${\mathbb{Z}}_q/(x-{\zeta }^{2k+1})$处的模数，即NTT变换后的NTT多项式系数，与上述${\hat{f}}_k$一致。
对于NTT，当有限域${\mathbb{Z}}_q$中包含$2n$次本原单位根时，$({\zeta }_{2n},{\zeta }_{2n}^3,\cdots ,{\zeta }_{2n}^{2n-1})\in {\mathbb{Z}}_q$，$(f({\zeta }_{2n}),f({\zeta }_{2n}^3),\cdots ,f({\zeta }_{2n}^{2n-1}))\in {\mathbb{Z}}_q$，所也可以如此进行计算。但若不存在$2n$次本原单位根，只存在$n$次本原单位根时，该方法便无法奏效，因为${\mathbb{Z}}_q$中不存在这样的函数对，所以无法以此方式来表示多项式。
如若为第二种情形则稍微复杂一些，此时通过第一章的中国剩余定理的方式更好理解NTT。NTT逆变换也可通过矩阵的方式得到，此时NTT变化对应的矩阵变为：$$\left(\begin{array}{c}{\hat{f}}_0\\ {\hat{f}}_2\\ ⋮\\ {\hat{f}}_{n-2}\\ {\hat{f}}_1\\ ⋮\\ {\hat{f}}_{n-1}\end{array}\right)=\left(\begin{array}{cc}A& 0\\ 0& A\end{array}\right)\cdot \left(\begin{array}{c}{f}_0\\ f_2\\ ⋮\\ f_{n-2}\\ f_1\\ ⋮\\ f_{n-1}\end{array}\right)$$ 其中：$$A={\left(\begin{array}{cccc}1& {\zeta }_n^{(2b{r}_{\log n-1}(0)+1)\times 1}& \cdots & {\zeta }_n^{(2b{r}_{\log n-1}(0)+1)\times \frac{n-2}{2}}\\ 1& {\zeta }_n^{(2b{r}_{\log n-1}(1)+1)\times 1}& \cdots & {\zeta }_n^{(2b{r}_{\log n-1}(0)+1)\times \frac{n-2}{2}}\\ ⋮& ⋮& ⋮& ⋮\\ 1& {\zeta }_n^{(2b{r}_{\log n-1}(\frac{n-2}{2})+1)\times 1}& \cdots & {\zeta }_n^{(2b{r}_{\log n-1}(\frac{n-2}{2})+1)\times \frac{n-2}{2}}\end{array}\right)}_{\frac{n}{2}\times \frac{n}{2}}$$ 由分块矩阵的求逆公式有$${\left(\begin{array}{cc}A& 0\\ 0& A\end{array}\right)}^{-1}=\left(\begin{array}{cc}{A}^{-1}& 0\\ 0& A^{-1}\end{array}\right)$$所以上述逆变换变为了求A的逆，其中A是一个范德蒙矩阵，所以逆矩阵为：$$A^{-1}={\left(\begin{array}{cccc}1& {\zeta }_n^{-(2b{r}_{\log n-1}(0)+1)\times 1}& \cdots & {\zeta }_n^{-(2b{r}_{\log n-1}(0)+1)\times \frac{n-2}{2}}\\ 1& {\zeta }_n^{-(2b{r}_{\log n-1}(1)+1)\times 1}& \cdots & {\zeta }_n^{-(2b{r}_{\log n-1}(0)+1)\times \frac{n-2}{2}}\\ ⋮& ⋮& ⋮& ⋮\\ 1& {\zeta }_n^{-(2b{r}_{\log n-1}(\frac{n-2}{2})+1)\times 1}& \cdots & {\zeta }_n^{-(2b{r}_{\log n-1}(\frac{n-2}{2})+1)\times \frac{n-2}{2}}\end{array}\right)}_{\frac{n}{2}\times \frac{n}{2}}$$ 所以此时NTT逆变换为:$$f_{2i}=\frac{2}{n}\sum _{j=0}^{n/2-1}{\hat{f}}_{2j}\cdot {\zeta }_n^{-(2b{r}_{\log n-1}(i)+1)j};$$ $$f_{2i+1}=\frac{2}{n}\sum _{j=0}^{n/2-1}{\hat{f}}_{2j+1}\cdot {\zeta }_n^{-(2b{r}_{\log n-1}(i)+1)j}.$$

三、快速离散傅里叶变换(FFT)

接下来介绍如何快速计算这些函数值。首先注意到单位根有如下性质($n$为2的方幂)：$$w_n^{2k}=e^{2\pi i\cdot 2k/n}=e^{2\pi i\cdot k/(n/2)}=w_{n/2}^k;$$ $$w_n^{k+n/2}=e^{2\pi i\cdot (k+n/2)/n}=e^{2\pi i\cdot k/n+\pi i}=-e^{2\pi i\cdot k/n}=-w_n^k.$$ 同理，有限域${\mathbb{Z}}_q$中的$n$次本原单位根（$n$整除$q-1$）也有该性质，令$\xi$为${\mathbb{Z}}_q$乘法群中的生成元，$n$次单位根即为${\zeta }_n={\xi }^{(q-1)/n}$,于是有：$${\zeta }_n^{2k}={\xi }^{\frac{q-1}{n}\times 2k}={\xi }^{\frac{q-1}{n/2}\times k}={\zeta }_{n/2}^k;$$ $${\zeta }_n^{k+n/2}={\xi }^{\frac{q-1}{n}\times (k+n/2)}={\xi }^{\frac{q-1}{n}\times k}\cdot {\xi }^{(q-1)/2}=-{\zeta }_n^k.$$接下来利用上述性质对$f(x)$进行处理。$$\begin{gathered} f(x)=f_0+f_{1}x+f_2{x}^2+\cdots +f_{n-2}{x}^{n-2}+f_{n-1}{x}^{n-1} \\ =f_0+f_2{x}^2+\cdots +f_{n-2}{x}^{n-2}+x(f_1+f_3{x}^2+\cdots +f_{n-1}{x}^{n-2}):=E_1(x^2)+x{O}_1(x^2) \end{gathered}$$ 将$f(x)$按系数分为奇偶两部分，分别记为$E_1(x^2)$和$O_1(x^2)$，其中$$E_1(X)=f_0+f_{2}X+\cdots f_{2i}{X}^i+\cdots +f_{(}n-2)X^{(n-2)/2};$$ $$O_1(X)=f_1+f_{3}X+\cdots f_{2i+1}{X}^i+\cdots +f_{(}n-1)X^{(n-2)/2}.$$ 因此计算一个$n-1$次多项式的值变为计算两个$(n-2)/2$次多项式的函数值，当$k<n/2$有$f(w_n^k)=E_1(w_n^{2k})+w_n^k\times O_1(w_n^{2k})=E_1(w_{n/2}^k)+w_n^k\times O_1(w_{n/2}^k)$；当$k^{\prime }=k+n/2$有$f(w_n^{k^{\prime }})=E_1(w_n^{2k^{\prime }})+w_n^{k^{\prime }}\times O_1(w_n^{2k^{\prime }})=E_1(w_{n/2}^k)-w_n^k\times O_1(w_{n/2}^k)$。所以计算$n$个函数值变为了计算$n/2$个函数值，且复杂度降低一半。然后继续往下分解函数，共$or{d}_2(n)$次后不可再分解，因为此时的单位根不能继续往下分解，且复杂度不断降低一半，故两个$n-1$次多项式的乘法计算复杂度由原来的$O(n^2)$变为$O(n\log n)$。
以上便是快速傅里叶变换（FFT）的计算思路，以下以一个7次多项式为例，即模多项式为$x^8+1$，$f(x)=f_0+f_{1}x+f_2{x}^2+\cdots +f_7{x}^7$，模数$q$满足$8|(q-1)$。同样先考存在$2n$次本原单位根，即16次本原单位根，此时NTT可仿照FFT的计算思路。
接下来计算$f({\zeta }_{16}^k)$，$(k=1,3,5,7,9,11,13,15)$，由于$$\begin{gathered} f(x)=f_0+f_{1}x+f_2{x}^2+f_3{x}^3+f_4{x}^4+f_5{x}^5+f_6{x}^6+f_7{x}^7 \\ =f_0+f_2{x}^2+f_4{x}^4+f_6{x}^6+x(f_1+f_3{x}^2+f_5{x}^4+f_7{x}^6) \\ =f_0+f_4{x}^4+x^2(f_2+f_6{x}^4)+x(f_1+f_5{x}^4+x^2(f_3+f_7{x}^4)) \end{gathered}$$ 分别计算$f_0+f_4{x}^4$，$f_2+f_6{x}^4$，$f_1+f_5{x}^4$，$f_3+f_7{x}^4$在${\zeta }_1{6}^k$时的取值，由于${\zeta }_{16}^{4}k={\zeta }_4^k$，因此只需计算上面4个多项式在${\zeta }_4^1$和$\zeta {\zeta }_4^3$处的取值即可，又因为${\zeta }_4^1=-{\zeta }_4^3$，因此只需计算其中一个即可。即$$\begin{gathered} f_0^{\prime }=f_0+f_4\cdot {\zeta }_4^1,f_4^{\prime }=f_0-f_4\cdot {\zeta }_4^1; \\ {f}_2^{\prime }=f_2+f_6\cdot {\zeta }_4^1,f_6^{\prime }=f_2-f_6\cdot {\zeta }_4^1; \\ {f}_1^{\prime }=f_1+f_5\cdot {\zeta }_4^1,f_5^{\prime }=f_1-f_5\cdot {\zeta }_4^1; \\ {f}_3^{\prime }=f_3+f_7\cdot {\zeta }_4^1,f_7^{\prime }=f_3-f_7\cdot {\zeta }_4^1. \end{gathered}$$ 再往上一层，即$f_0+f_4{x}^4+x^2(f_2+f_6{x}^4)$和$f_1+f_5{x}^4+x^2(f_3+f_7{x}^4)$。此处需要在原来的基础上乘$x^2$，即${\zeta }_{16}^{2k}={\zeta }_8^k$，又因为${\zeta }_8^1=-{\zeta }_8^5$，${\zeta }_8^3=-{\zeta }_8^7$，而$x^4$处，即$(x^2{)}^2$又为${\zeta }_4^1=-{\zeta }_4^3$，所以总的有8个值，即：$$\begin{gathered} f_0^{\prime \prime }=f_0^{\prime }+{\zeta }_8^1\cdot f_2^{\prime },f_2^{\prime \prime }=f_0^{\prime }-{\zeta }_8^1\cdot f_2^{\prime }; \\ {f}_4^{\prime \prime }=f_4^{\prime }+{\zeta }_8^3\cdot f_6^{\prime },f_6^{\prime \prime }=f_4^{\prime }-{\zeta }_8^3\cdot f_6^{\prime }; \\ {f}_1^{\prime \prime }=f_1^{\prime }+{\zeta }_8^1\cdot f_3^{\prime },f_3^{\prime \prime }=f_1^{\prime }-{\zeta }_8^1\cdot f_3^{\prime }; \\ {f}_5^{\prime \prime }=f_5^{\prime }+{\zeta }_8^3\cdot f_7^{\prime },f_7^{\prime \prime }=f_5^{\prime }-{\zeta }_8^3\cdot f_7^{\prime }; \end{gathered}$$ 最后，计算$f_0+f_4{x}^4+x^2(f_2+f_6{x}^4)+x(f_1+f_5{x}^4+x^2(f_3+f_7{x}^4))$，即$f(x)$。在上述的基础上此处出现了$x$，即${\zeta }_{16}^k$，因为${\zeta }_{16}^1=-{\zeta }_{16}^9$，${\zeta }_{16}^3=-{\zeta }_{16}^{11}$，${\zeta }_{16}^5=-{\zeta }_{16}^{13}$，${\zeta }_{16}^7=-{\zeta }_{16}^{15}$。而平方后的值均已经计算了，所以这8个值为：$$\begin{gathered} F_0=f({\zeta }_{16}^1)=f_0^{\prime \prime }+{\zeta }_{16}^1\cdot f_1^{\prime \prime }; \\ {F}_4=f({\zeta }_{16}^9)=f_0^{\prime \prime }-{\zeta }_{16}^1\cdot f_1^{\prime \prime }; \\ {F}_2=f({\zeta }_{16}^5)=f_2^{\prime \prime }+{\zeta }_{16}^5\cdot f_3^{\prime \prime }; \\ {F}_6=f({\zeta }_{16}^{13})=f_2^{\prime \prime }-{\zeta }_{16}^5\cdot f_3^{\prime \prime }; \\ {F}_1=f({\zeta }_{16}^3)=f_4^{\prime \prime }+{\zeta }_{16}^3\cdot f_5^{\prime \prime }; \\ {F}_5=f({\zeta }_{16}^{11})=f_4^{\prime \prime }-{\zeta }_{16}^3\cdot f_5^{\prime \prime }; \\ {F}_3=f({\zeta }_{16}^7)=f_6^{\prime \prime }+{\zeta }_{16}^7\cdot f_7^{\prime \prime }; \\ {F}_7=f({\zeta }_{16}^{15})=f_6^{\prime \prime }-{\zeta }_{16}^7\cdot f_7^{\prime \prime }; \end{gathered}$$接下来将在下一篇文章中介绍NTT中与FFT类似的快速计算方法以及相应伪代码。欢迎大家交流指正。