数论变换（NTT）及其变体

参考文献：

1. Chung C M M, Hwang V, Kannwischer M J, et al. NTT multiplication for NTT-unfriendly rings: New speed records for Saber and NTRU on Cortex-M4 and AVX2[J]. IACR Transactions on Cryptographic Hardware and Embedded Systems, 2021: 159-188.
2. Boemer F, Kim S, Seifu G, et al. Intel HEXL: accelerating homomorphic encryption with intel AVX512-IFMA52[C]//Proceedings of the 9th on Workshop on Encrypted Computing & Applied Homomorphic Cryptography. 2021: 57-62.

CRT

令$R$是基环 (base ring)，那么$R[x]$是一元多项式环。

如果$f,g\in R[x]$且$gcd(f,g)=1$，那么存在环同构：
$$\begin{array}{rl}\varphi :R[x]/(f(x)g(x))& \cong R[x]/(f(x))\times R[x]/(g(x))\\ \varphi (h)& =(h\mathrm{mod}f,h\mathrm{mod}g)\end{array}$$
且：
$$\begin{array}{rl}if& \\ & \varphi (h)=(h_f,h_g)\\ let& \\ & g_f^{-1}\cdot g\equiv 1\mathrm{mod}f\\ & f_g^{-1}\cdot f\equiv 1\mathrm{mod}g\\ then& \\ & h\equiv h_f\cdot g\cdot g_f^{-1}+h_g\cdot f\cdot f_g^{-1}\mathrm{mod}fg\end{array}$$

(cyclic) NTT 【模$x^{2^k}-1$】

• FFT是DFT的快速算法，对于$n$长向量，计算复杂度为$O(n\log n)$

• NTT是工作在素域上$R=F_p$的FFT，因此它不含复数运算，计算结果没有舍入误差。

• 令$n=2^k$，NTT需要寻找$x^n-1$的本原单位根：$\xi \in F_p$

• $\exists \xi \in F_p⟺n|p-1$，形如$p=2^k{p}^{\prime }+1$的素数，叫做NTT友好的素数 (NTT-friendly prime)

• 给定本原元$g\in F_p$，那么n次本原单位根为：${\xi }_n=g^{\frac{p-1}{n}}$

• 易知：$gcd(x-{\xi }^i,x-{\xi }^j)=1,i\not\equiv j$，利用CRT：
  $$\begin{array}{r}\varphi :R[x]/(x^n-1)\to R[x]/(x-1)\times R[x]/(x-\xi )\times \cdots \times R[x]/(x-{\xi }^{n-1})\end{array}$$

• 如果$f(x)={\sum }_{i=0}^n{f}_i{x}^i$，那么$f(x)\equiv f(a)\mathrm{mod}x-a$；因此，上式为：$\varphi (f)=[f(1),f(\xi ),\cdots ,f({\xi }^{n-1})]$

• 当我们只关注多项式乘积时，分解的顺序不重要。由于$n=2^k$，可利用CT蝴蝶或者GS蝴蝶实现NTT，注意比特反序。

• 由于DFT和IDFT公式相似，所以也可以用NTT计算IDFT：将$\xi$替换为${\xi }^{-1}$，并将结果再乘以$\frac{1}{n}$

• 将环元素$f(x)\in R[x]/(x^n-1)$写作$n$长的系数向量：$[f_{n-1},\cdots ,f_1,f_0]$

• 环元素乘法$h(x)=f(x)g(x)\mathrm{mod}{x}^n-1$ $⟺$ 系数向量循环卷积$h_k=f(k)⊛g(k):={\sum }_{t=0}^{n-1}{f}_{(k-t)}{g}_t$，这里$(k-t):=k-t\mathrm{mod}n$

• 对于$f,g\in R[x]/(x^n-1)$，做NTT分别得到$[f(\xi ),\cdots ,f({\xi }^{n-1}),f({\xi }^n)]$和$[g(\xi ),\cdots ,g({\xi }^{n-1}),g({\xi }^n)]$

• $f(x)g(x)⟺[f({\xi }^j)g({\xi }^j)]$，左边是多项式乘积(向量卷积)，右边是阵列乘积。

NTT变体

Negacyclic NTT【模$x^{2^{k-1}}+1,2^k\mid q-1$】

• 对于$R[x]/(x^n+1)$，其上多项式乘法是反循环卷积 (Negacyclic convolution)。

• 第一种解法 （直接使用NTT第$j=1$层的下半部分）

  • 如果$n=2^{k-1}$，寻找$2n$次本原单位根：${\xi }_{2n}\in R$；由于$(x^{2n}-1)=(x^n+1)(x^n-1)$，且$gcd(x^n+1,x^n-1)=1$，因此$R[x]/(x^{2n}-1)\cong R[x]/(x^n+1)\times R[x]/(x^n-1)$

  • $x^{2n}-1$的所有根是 { ξ 2 n , ξ 2 n 2 , ⋯   , ξ 2 n 2 n } \{\xi_{2n},\xi_{2n}^2,\cdots,\xi_{2n}^{2n}\} {ξ2n​,ξ2n2​,⋯,ξ2n2n​}，而$x^n-1$的所有根是 { ξ 2 n 2 , ξ 2 n 4 , ⋯   , ξ 2 n 2 n } \{\xi_{2n}^2,\xi_{2n}^4,\cdots,\xi_{2n}^{2n}\} {ξ2n2​,ξ2n4​,⋯,ξ2n2n​}，因此$x^n+1$的所有根是 { ξ 2 n 1 , ξ 2 n 3 , ⋯   , ξ 2 n 2 n − 1 } \{\xi_{2n}^1,\xi_{2n}^3,\cdots,\xi_{2n}^{2n-1}\} {ξ2n1​,ξ2n3​,⋯,ξ2n2n−1​}

  • $R[x]/(x^n+1)\cong R[x]/(x-\xi )\times \cdots \times R[x]/(x-{\xi }^{2n-1})$

  • 模$x^n+1$的Negacyclic NTT，可以看做模$x^{2n}-1$的(cyclic) NTT的一半：先对长度$n$的 { f i } \{f_i\} {fi​}补零到$2n$长序列，计算NTT，截取结果 { f ( ξ 1 ) , f ( ξ 3 ) , ⋯   , f ( ξ 2 n − 1 ) } \{f(\xi^1),f(\xi^3),\cdots,f(\xi^{2n-1})\} {f(ξ1),f(ξ3),⋯,f(ξ2n−1)}

  • 对于IDFT，直接用上述 Negacyclic NTT 实现即可。

• 第二种解法 (感觉更好理解)

  • 令$Y={\xi }_{2n}X$，则$R[X]/(X^n+1)\cong R[Y]/(Y^n-1)$
    $$\begin{array}{rl}\varphi (\sum _{i=0}^{n-1}{f}_i{x}^i):=& \sum _{i=0}^{n-1}{f}_i{\xi }^{-i}(\xi x{)}^i=\sum _{i=0}^{n-1}(f_i{\xi }^{-i})y^i\\ {\varphi }^{-1}(\sum _{i=0}^{n-1}{g}_i{y}^i):=& \sum _{i=0}^{n-1}{g}_i{\xi }^i({\xi }^{-1}y{)}^i=\sum _{i=0}^{n-1}(g_i{\xi }^i)x^i\end{array}$$

  • 先把$f(x)\in R[X]/(X^n+1)$映射到$g=\varphi (f)\in R[Y]/(Y^n-1)$，然后做(cyclic) NTT得到 { g ( 1 ) , g ( ξ ) , ⋯   , g ( ξ n − 1 ) } \{g(1),g(\xi),\cdots,g(\xi^{n-1})\} {g(1),g(ξ),⋯,g(ξn−1)}

  • 对于IDFT，先用NTT从 { g ( ξ j ) } \{g(\xi^j)\} {g(ξj)}得到$g\in R[Y]/(Y^n-1)$，然后输出$f={\varphi }^{-1}(g)$

Incomplete NTTs【模$x^{h\cdot 2^k}-1,2^k\mid q-1$】

• 对于$f\in R[x]/(x^n-1),n=h\cdot 2^k,2^k\mid q-1$，使用FFT对它做$k$层分解，得到$2^k$个$R[x]/(x^h-{\xi }_i)$的直积：$[f\mathrm{mod}{x}^h-{\xi }_i,\cdots ]$
• 对于$f,g\in R[x]/(x^n-1)$，做不完全NTT，分别得到$[f_1,\cdots ,f_{2^k}]$和$[g_1,\cdots ,g_{2^k}]$
• 大卷积$f(x)g(x)$ $⟺$ 若干小卷积$[f_j(x)g_j(x)]$，左右两边都是多项式乘积
• 对于固定参数$p,n$的商环$F_p[x]/(x^n-1)$，根据硬件实现的差异 (寄存器、运算器的数目，支持运算数的比特长度)，选择合适的分解层数，对于小卷积直接使用经典算法即可。

Good’s Trick【模$x^{h\cdot 2^k}-1,2^k\mid q-1$】

• 对于$R[x]/(x^n-1)$，假如$n=h\cdot 2^k$，其中$h$是奇数。我们令$x=yw$，且$y^{2^{k-1}}=-1,w^{h-1}+w^{h-2}+\cdots +w=-1$

• 做映射：
  $$\begin{array}{rl}\varphi :& R[x]/(x^{h\cdot 2^k})\to R[y,w]\\ & x^i\mapsto y^{i\mathrm{mod}{2}^k}{w}^{i\mathrm{mod}h}\end{array}$$
  令$b(y,w)=\varphi (a(x))$，则：${\mathrm{deg}}_{y}b<2^k,{\mathrm{deg}}_{w}b<h$

• 记$b(y,w)={\sum }_{i=0}^{h-1}{w}^i{b}_i(y)$，其中$b_i(y)\in R[y]/(y^{2^k}-1)$

• Good’s Permutation：用一维阵列$a[i]$记录$a(x)={\sum }_{i=0}^{h\cdot 2^k-1}{a}_i{x}^i$；用二维阵列$b[i][j]$记录$b(y,w)={\sum }_{i=0}^{h-1}{\sum }_{j=0}^{2^k-1}{b}_{ij}{w}^i{y}^j$，其中的每一行$b[i]$记录了$b_i(y)$

• $a\in R[x]/(x^{h\cdot 2^k}-1)$，我们映射到$b=\varphi (a)={\sum }_{i=0}^{h-1}{w}^i{b}_i(y)$；然后对每一行$b[i]$，并行地做$h$次大小$2^k$的NTT，得到NTT域结果：阵列$B[i][j]$

• 对于$b_1(y,w)={\sum }_{i=0}^{h-1}{w}^i{b}_i^{(1)}(y)$和$b_2(y,w)={\sum }_{i=0}^{h-1}{w}^i{b}_i^{(2)}(y)$，乘积为：
  $$b_1(y,w)b_2(y,w)=\sum _{t=0}^{h-1}\sum _{i=0}^{h-1}[b_{t-i}^{(1)}(y)b_i^{(2)}(y)]w^t\mathrm{mod}{w}^h-1$$
  且：$b^{(3)}(y)={\sum }_{i=0}^{h-1}{b}_{t-i}^{(1)}(y)b_i^{(2)}(y)⟺\forall j,b^{(3)}({\xi }^j)={\sum }_{i=0}^{h-1}{b}_{t-i}^{(1)}({\xi }^j)b_i^{(2)}({\xi }^j)$

• 多项式乘积

  • 对于$a_1,a_2\in R[x]/(x^{h\cdot 2^k}-1)$，为了计算$a_1(x)a_2(x)$，我们先得到$B_1[i][j],B_2[i][j]$
  • 阵列的第$j$列，表示定点$y={\xi }^j$处的关于$w$的模$w^h-1$的多项式
  • 对$B_1$和$B_2$的第$j$列做循环卷积 (这对应于$h=1$时的阵列乘)，也就是定点$y={\xi }^j$处的关于$w$的两个多项式乘积：$b_1(y,w)b_2(y,w)⟺[{\sum }_{i=0}^{h-1}{b}_{t-i}^{(1)}({\xi }^j)b_i^{(2)}({\xi }^j)\mathrm{mod}{w}^h-1]$
  • 最后，对结果阵列$B_3$做并行的$h$个INTT，得到$b_3=b_1{b}_2$，并输出结果$a_3={\varphi }^{-1}(b_3)=a_1{a}_2$

• 复杂度：$2hO(2^k\log 2^k)+2^{k}O(h^2)=O(2n\log \frac{n}{h})+O(nh)$；当$h\to 1$时为$O(n\log n)$，当$h\to n$时为$O(n^2)$

  也就是说，对于$n=h\ast 2^0$，对于$R[x]/(x^n-1)$上的多项式乘积 NTT 没用？奥，我们可以模$x^{n^{\prime }}-1$嘛，保证$n^{\prime }>2n$就行了。

• 为什么我觉得它和不完全NTT一模一样呢？

环嵌入【模$q\ne 2^k{p}^{\prime }+1$】

• 如果我们需要计算$Z_q[x]/(x^n-1)$上的多项式乘积，但是$n=2^k\nmid q-1$，即$q\ne 2^k{p}^{\prime }+1$，于是$Z_q$中不存在n次本原单位根。
• 对于$f,g\in Z_q[x]/(x^n-1)$，令它们的系数范围是：$f_i,g_i{\text{ mod}}^{\pm }q\in [-\frac{q}{2},\frac{q}{2})$。容易知道，多项式乘积$f(x)g(x)$的每一项的绝对值大小至多是$\frac{n{q}^2}{4}$
• 我们选取一个NTT友好的大素数$p>\frac{n{q}^2}{2}$，满足$n\mid p-1$
• 我们把$f,g\in Z_q[x]/(x^n-1)$映射到$f^{\prime },g^{\prime }\in Z_p[x]/(x^n-1)$，这里就是简单的系数嵌入：$\varphi :a\in Z_p\mapsto a\in Z_q$
• 将$f,g$嵌入到$Z_p[x]/(x^n-1)$中得到$f^{\prime },g^{\prime }$，计算多项式乘积，然后简单取模即可：$f(x)g(x)\equiv f^{\prime }(x)g^{\prime }(x)\mathrm{mod}q$
• 由于$p>\frac{n{q}^2}{2}$，因此乘积的系数不会模掉$p$，所以上述结果一定正确。由于$n\mid p-1$，因此可以利用NTT加速多项式乘法。

多模数【模$q\ne 2^k{p}^{\prime }+1$】

• 这是解决$n=2^k\nmid q-1$的另一种方案。

• 我们挑选一系列的NTT友好的小素数 { p j } \{p_j\} {pj​}，满足$n\mid p_j-1$且$P={\prod }_j{p}_j>\frac{n{q}^2}{2}$

• 我们把$f,g\in Z_q[x]/(x^n-1)$映射到$f^{\prime },g^{\prime }\in Z_P[x]/(x^n-1)$，这里也是简单的系数嵌入。计算$f^{\prime }(x)g^{\prime }(x)\mathrm{mod}q$即可，且由于$P>\frac{n{q}^2}{2}$，所以结果一定正确。

• 根据CRT，定义如下双射 (因为$gcd(p_i,p_j)=1$)：
  $$\begin{array}{r}\varphi :f^{\prime }\mathrm{mod}P\mapsto [f^{\prime }\mathrm{mod}{p}_1,\cdots ,f^{\prime }\mathrm{mod}{p}_s]\end{array}$$
  利用CRT的求解公式，当然是可以的；但有一种更高效的迭代算法，Garner’s Algorithm：

  假设有同余方程组$u\equiv u_i\mathrm{mod}{p}_i$，且满足：$|u_i|<\frac{p_i}{2},|u|<\frac{P}{2}$

  那么，
  $$\begin{array}{rl}{y}_1& =u_1\\ y_2& =y_1+((u_2-y_1)m_2\mathrm{mod}{}^{\pm }{p}_2)\cdot p_1\\ y_3& =y_2+((u_3-y_2)m_3\mathrm{mod}{}^{\pm }{p}_3)\cdot p_1{p}_2\\ & ⋮\\ y_s& =y_{s-1}+((u_s-y_{s-1})m_s\mathrm{mod}{}^{\pm }{p}_s)\cdot p_1{p}_2\cdots p_{s-1}\end{array}$$
  其中$m_i:=(p_1{p}_2\cdots p_{i-1}{)}^{-1}\mathrm{mod}{}^{\pm }{p}_i$，那么$u=y_s$

  算法结果正确，且对于小的$s$，上述迭代算法比CRT更高效。

• 计算多项式乘法

  • 将$f,g$嵌入到$Z_p[x]/(x^n-1)$中得到$f^{\prime },g^{\prime }$
  • 做映射，得到$[f_i^{\prime }(x)\mathrm{mod}{p}_i],[g_i^{\prime }(x)\mathrm{mod}{p}_i]$
  • 对每个对应的分量，分别计算多项式乘积 (利用NTT)，得到序列$[(f^{\prime }{g}^{\prime }{)}_i(x)\mathrm{mod}{p}_i]$
  • 求解方程组$(f^{\prime }{g}^{\prime })\equiv (f^{\prime }{g}^{\prime }{)}_i\mathrm{mod}{p}_i$，最后$(fg)\equiv (f^{\prime }{g}^{\prime })\mathrm{mod}q$

总结

• 用快速数论变换 ( NTT) 实现时域、频域之间的快速转换。使用DFT的卷积性质，可加速$Z_q[x]/(x^n-1)$上多项式乘积的计算。时间复杂度：$2O(n\log n)+O(n)=O(n\log n)$
• 一、对于环$Z_q[x]/(x^n-1),n\mid q-1,n=2^k$，用NTT可以直接求解。
• 二、对于环$Z_q[x]/(x^n+1),n\mid q-1,n=2^{k-1}$，用Negacyclic NTT，做映射，在$Z_q[Y]/(Y^n-1)$上利用NTT求解。
• 三、对于环$Z_q[x]/(x^n-1),n\mid q-1,n=h\cdot 2^k$，用Good’s Trick，将一元多项式映射为二元多项式 (以一元多项式为系数的一元多项式)，然后利用NTT求解。
• 四、对于环$Z_q[x]/(x^n-1),n\nmid q-1,n=2^k$，用环嵌入，找到一个NTT友好的素数$n\mid p-1$，然后用NTT求解。
• 五、对于环$Z_q[x]/(x^n-1),n\nmid q-1,n=2^k$，用多模数，找到一系列NTT友好的素数$n\mid p_i-1$，然后用NTT分别求解，并用CRT组合结果。
• 六、对于环$Z_q[x]/(x^n-1),n\mid q-1,n=h$，$h$是奇数，那么直接用Good’s Trick对多项式乘积没有帮助。我们可以寻找$n^{\prime }=h^{\prime }\times 2^k>2n$，并寻找相应的$n^{\prime }\mid q^{\prime }-1$；然后在$Z_{q^{\prime }}[x]/(x^{n^{\prime }}-1)$上计算多项式乘积，最后结果$fg=(f^{\prime }{g}^{\prime })\mathrm{mod}q,x^n-1$即可。 要保证两个n长多项式一次相乘后不会模$x^{n^{\prime }}-1$，以保证最终结果正确，所以要$n^{\prime }>2n$
• 七、对于环$Z_q[x]/(x^n-1),n\nmid q-1,n=h\cdot 2^k$，先用环嵌入或者多模数，使得$n\mid p-1$或者$n\mid p_i-1$；这样问题就转化为了若干个问题“三”；用Good’s Trick，将一元多项式映射为二元多项式，然后利用NTT求解。
• 我们可以解决“六”和“七”，也就是说，对于任意的参数$q,n$，环$Z_q[x]/(x^n\pm 1)$上多项式乘积都可以用NTT快速计算。