NTT

本文总结了与NTT相关的容易混淆的问题：DIT、DIF；原址计算、流水线；正常顺序、比特逆序；正包卷积、负包卷积。

---

数论变换（NTT）是一种用来计算有限域上卷积的快速算法。由于多项式乘法等效于计算系数的卷积，NTT常用于加速计算多项式乘法。

1. NTT定义与使用条件

定义

数论变换正变换的公式为：

\[F[k] = \sum\limits_{n=0}^{N-1} f[n]\omega_N^{nk}\mod q,\quad k=0,1,2,...,N-1\]

数论变换逆变换(INTT)的公式为：

\[f[n] = N^{-1}\sum\limits_{k=0}^{N-1} F[k]\omega_N^{-nk}\mod q,\quad n=0,1,2,...,N-1\]

符号解释：

• q为有限域的模数，有限域内所有计算均要模q
• N为进行数论变换的数字的个数，也称为NTT的点数
• $f[n]$为变换前的第n个数，$F[k]$为变换后的第k个数
• $\omega_N$为模q有限域的N阶单位根，阶数N用下标表示，需要满足两个条件：
  • $\omega_N^N=1\mod q$
  • $\omega_N^k\neq1\mod q,\quad k=1,2,…,N-1$
• $N^{-1}$可以理解为有限域内的$1/N$，满足：$N\cdot N^{-1}=1\mod q$

NTT的使用条件：

• 模数q需要是质数
• $N$必须是$q-1$的因数

用NTT来计算卷积的方法：

• 输入：N个数的数组a，N个数的数组b
• 计算：数组a与数组b的卷积
• 输出：2N个数的卷积结果数组c
• 计算方法：
  • $A=NTT_{2N}(零填充(a))$，零填充指在前面填充N个零，使数组长度变为2N
  • $B=NTT_{2N}(零填充(b))$
  • $C=A\odot B$，即$C[k] = A[k]\cdot B[k],\quad k=0,1,2,…,2N-1$
  • $c=INTT_{2N}(C)$
• 注：当卷积结果需要模多项式时，有时可通过正负包卷积加速，见后文。

2. DIT与DIF

为了公式简洁，以下计算中不再标记模q，所有计算默认在模q有限域内。

由于单位根满足以下的这些性质，NTT也可以像FFT用分治思想加速：

• 对称性：$\omega_N^{k+N/2}=-\omega_N^k$
• 周期性：$\omega_N^{k+N} = \omega_N^k$
• 放缩性：$\omega_{N/m}^{k/m} = \omega_N^k$

Radix-2 DIT NTT

基2的DIT将输入（对应DFT的时域）按奇偶分为2组，将$F[k] = \sum\limits_{n=0}^{N-1} f[n]\omega_N^{nk}$分解为：

\[F[k]=NTT_{N/2}(f[2i],k)+\omega_N^kNTT_{N/2}(f[2i+1],k)\]

由于单位根的对称性与周期性，对于$k’=k+N/2>N/2$的$F[k’]$可以简化计算为：

\[\begin{eqnarray} F[k+N/2] &=& NTT_{N/2}(f[2i],k+N/2&)&+&\omega_N^{k+N/2}&&NTT_{N/2}(f[2i+1],&k+N/2&)\\ &=& NTT_{N/2}(f[2i],k&)&-&\omega_N^k&&NTT_{N/2}(f[2i+1],&k&) \end{eqnarray}\]

当N=8时的示意图：

将此分解方式不断重复，最终得到：

Radix-2 DIF NTT

基2的DIF将输出（对应DFT的频域）按奇偶分为2组

利用放缩性与周期性将偶数输出变形为

\[F[2k]=NTT_{N/2}(f[i]+f[i+N/2],k)\]

利用放缩性、周期性与对称性将奇数输出变形为

\[F[2k+1]=NTT_{N/2}((f[i]-f[i+N/2])\omega_N^i,k)\]

当N=8时的示意图：

将此分解方式不断重复，最终得到：

蝶形运算

由以上DIT和DIF的最终分解图可以发现，NTT的运算最终被分解为一个个相似的运算单元，这些运算单元被称为蝶形运算，对应的单位根的幂$\omega_N^k$也称为蝶形因子或旋转因子。

• DIT NTT对应的蝶形运算形式为：$a\pm\omega\cdot b$

• DIF NTT对应的蝶形运算形式为：$a+b,\ \ (a-b)\cdot\omega$

INTT

与NTT几乎相同，只需做如下两个修改：

• 将所有蝶形因子，指数部分取负。例如将$\omega_8^1$变为$\omega_8^{-1}$，实际计算时，可利用周期性和对称性变换：$\omega_8^{-1}=\omega_8^7=-\omega_8^3$
• 在计算结束后，所有结果乘上$N^{-1}$，常称为后处理

也有学者提出可将后处理合并在蝶形运算中，参见论文¹

3. 原址计算(In-place)

由于NTT的点数N有时非常大，这时中间运算结果的数量非常大。为了节省存储空间，一种常用的方法是原址计算(In-place calculation)，即每次蝶形运算的结果存入蝶形运算输入变量的相同地址（覆盖原值），这样就不需要额外的存储空间来存储中间结果。

原址计算可参考上小节的蝶形图，每行的输入、中间结果、输出存在同一个地址，计算时必须保证左侧的蝶形运算先运算，以避免数据非法覆盖。

原址计算存在以下几个缺点：对存储器的读写次数很多、蝶形运算读写的地址较复杂、存在输入或输入比特逆序问题。

正常顺序与比特逆序

观察第2节中DIT和DIF NTT的蝶形图，可以发现输入或输出总有一个不是正常的顺序。

我们把0, 1, 2, …, N-1这样的顺序称为正常顺序（normal order ）。

将正常顺序的二进制写出，例如N=8时：000, 001, 010, 011, 100, 101, 110, 111

如果将二进制的index全部前后调转，会变为：000, 100, 010, 110, 001, 101, 011, 111

对应的顺序0, 4, 2, 6, 1, 5, 3, 7就称为比特逆序（bit-reversed order）。

原址计算的比特逆序问题

原址计算中的正常顺序和比特逆序，其实是指的数据在存储器中存储的顺序，即地址顺序与index对应的顺序是正常顺序还是比特逆序。并非指的参与运算的数的index，因为那是算法已经确定了的。

原址计算的特性是经过一次NTT/INTT，地址会被比特逆序。按照是DIT还是DIF，以及输入数据在存储器中是按正常顺序还是比特逆序存储的，可以将原址计算的NTT分为4类：$NTT_{no\to br}^{DIT}$, $NTT_{br\to no}^{DIT}$, $NTT_{no\to br}^{DIF}$, $NTT_{br\to no}^{DIF}$

其中，DIT还是DIF会影响蝶形运算的结构与蝶形因子，存储顺序则主要影响存取数的地址计算以及蝶形因子的index计算。其中$NTT_{br\to no}^{DIT}$及$NTT_{no\to br}^{DIF}$的蝶形因子index计算较简单（另两个index计算需要比特逆序），因此更为常用。

流水线

与原址计算不同，另一种硬件常用的计算结构为流水线结构。流水线结构在FFT领域中已得到充分的研究和应用，而在NTT中使用的较少。其特点是，只从存储器读一次输入数据、写一次输出数据，中间结果存储在大量的移位寄存器（延时单元）中，一般使用多个蝶形单元，每个单元负责一层的蝶形运算。大致可以分为三类：

• 单路延时反馈（SDF）架构
• 多路延时换向（MDC）架构
• 单路延时换向（SDC）架构

例图（来自论文²中图1）：

流水线结构由于只读写一次存储器，因此不存在逆序的问题，只需要按照所需index对应的地址去读写就行了。SDF与SDC的乘法器利用率是无法达到100%的，而MDC可以通过少许的修改（每周期输入多个数据）达到100%的乘法器利用率。

流水线结构的一个主要缺点是，需要使用大量的移位寄存器用作延时单元，论文³中研究了使用FPGA的BRAM替代移位寄存器的方法。

正包卷积与负包卷积

第1节中介绍了，使用NTT来加速两个长度为N的序列的卷积，需要进行零扩展及2N点的NTT/INTT。当我们计算的是多项式的乘法（即系数的卷积），且多项式在有限域$\mathbb{Z}_q[x]/g(x)$上时，最终的计算结果需要模多项式$g(x)$。当模多项式$g(x)$为两种特殊多项式时，可以用正负包卷积来加速。

当模多项式为$g(x)=x^N-1$时，可以用正包卷积/循环卷积（positive wrapped convolution/cyclic convolution）来加速，直接将原来的2N点NTT/INTT换为N点NTT/INTT即可：

\[c=INTT_N(NTT_N(a)\oplus NTT_N(b))\]

当模多项式为$g(x)=x^N+1$时，可以用负包卷积（negative wrapped convolution）来加速，但是情况相较正包卷积会复杂一些。在NTT之前，需要对输入多项式进行一个预处理；在INTT之后，需要对其结果进行后处理：

• 预处理：

  \[\hat a[k] = a[k]\cdot\omega_{2N}^k\\ \hat b[k] = b[k]\cdot\omega_{2N}^k\]

• NTT、逐点乘法（point-wise multiplication）、INTT：

  \[\hat c=INTT_N(NTT_N(\hat a)\oplus NTT_N(\hat b))\]

• 后处理：

  \[c[k]=\hat c[k]\cdot\omega_{2N}^{-k}\]

注意：负包卷积的前后处理中，用到了2N阶单位根。所以必须满足2N是q-1的因数，才能使用负包卷积。

另外，论文¹中给出了将负包卷积的前后处理合并进NTT/INTT的蝶形运算中的方法。

NTT的代数理解

• 多项式环$\mathbb Z_q[x]/(x^N-1)$上的乘法叫循环卷积，因为$c_k=\Sigma_{i=0}^ka_ib_{k-i}+\Sigma_{i=k+1}^{N-1}a_ib_{k+N-i}\bmod q$.

  • 利用NTT加速乘法，实际是利用了环同构：

    $\mathbb Z_q[x]/(x^N-1)\cong\mathbb Z_q[x]/(x^{\frac N2}-1)\times\mathbb Z_q[x]/(x^{\frac N2}+1)$.

    再利用$x^{\frac N2}+1=x^{\frac N2}-\omega_N^{\frac N2}$进一步分解，最终分解到N个一次多项式的商环上

    $\mathbb Z_q[x]/(x^N-1)\cong\Pi_{i=0}^{N-1}\mathbb Z_q[x]/(x-\omega_N^{br(i)})$.

    一次多项式的商环上乘法就是整数模乘，没有多项式操作

  • 此时只要求$\mathbb Z/q\mathbb Z$存在N阶原根，即乘法群$\mathbb Z_q^*$中阶为$N$的元素

    • 只有$q$是素数时，$\mathbb Z_q$才是有限域，所以要求$q$是素数
      • 质数存在原根$r$，即$r$的阶为$\phi(q)=q-1$⁴
      • 因此乘法群是循环群，原根$r$是其生成元

    • $q$是素数时，乘法群$\mathbb Z_q^\star=\mathbb Z_q\setminus\left{0\right}$的阶为$

      \mathbb Z_q^\star

      =q-1$

    • 元素的阶必是群的阶的因子，因此存在$\omega_N$的条件是$q-1$是$N$的倍数
      • 若已知原根$r$，可选择$\omega_N=r^{\frac{q-1}{N}}$

• 多项式环$\mathbb Z_q[x]/(x^N+1)$上的乘法叫负包卷积，因为$c_k=\Sigma_{i=0}^ka_ib_{k-i}-\Sigma_{i=k+1}^{N-1}a_ib_{k+N-i}\bmod q$.

  • 这里同样利用环同构，并利用$x^N+1=x^N-\omega_N^{N/2}$:

    $\mathbb Z_q[x]/(x^N+1)=\mathbb Z_q[x]/(x^N-\omega_N^{N/2})\cong\mathbb Z_q[x]/(x^{\frac N2}-\omega_N^{N/4})\times\mathbb Z_q[x]/(x^{\frac N2}+\omega_N^{N/4})$.

    但是仅利用$\omega_N$只能分解到2次多项式的商环上

    如果$\mathbb Z_q$上存在2N阶原根$\psi_{2N}$，2次多项式可以进一步分解，最终分解到：

    $\mathbb Z_q[x]/(x^N+1)\cong\Pi_{i=0}^{N-1}\mathbb Z_q[x]/(x-\psi_{2N}^{2br(i)+1})$.

  • 此时要求$q-1$是$2N$的倍数

• $N^{-1}$是乘法群$\mathbb Z_q^*$中$N$的逆，可通过扩展欧几里得算法求