FFT-trick及Culey-Tukey蝴蝶变换原理及实现

FFT-trick计算多项式环乘法

定理：

${\mathbb{Z}}_q/(a(x)\cdot b(x))\cong {\mathbb{Z}}_q/(a(x))\times {\mathbb{Z}}_q/(b(x))$,其中$a(x),b(x)$互素

$\varphi :f(x)\to (f(x)moda(x),f(x)modb(x))$

• ${\mathbb{Z}}_q[x]/(x^n+1)={\mathbb{Z}}_q[x]/(x^n-{\omega }^n)\cong {\mathbb{Z}}_q[x]/(x^{\frac{n}{2}}-{\omega }^{\frac{n}{2}})\times {\mathbb{Z}}_q[x]/(x^{\frac{n}{2}}+{\omega }^{\frac{n}{2}})$
• $f={\sum }_{i=0}^{n-1}{f}_i{x}^i\to (fmod(x^{\frac{n}{2}}-{\omega }^{\frac{n}{2}}),fmod(x^{\frac{n}{2}}+{\omega }^{\frac{n}{2}}))$

当$n=4$的时候则可以得到以下分解

快速数论变换NTT的结果都可以以这种方式计算出做完NTT变换后的每项系数，也就是对应这个二叉树叶节点的值（都为常数）。其实求每个分支的值的方法也非常简单。例如我们要求第二层也就是根节点的左右儿子的值，左边：只需要将$x^2={\omega }^2$带入根节点的式子中即可得左儿子的值。

这里解释一下为什么根节点的右儿子可以分解成${\mathbf{Z}}_q[x]/\left(x-{\omega }^3\right)$和${\mathbf{Z}}_q[x]/\left(x+{\omega }^3\right)$的形式，是因为这里$\omega$是模$q$的本原单位根，满足${\omega }^n=-1$的性质，因此$x^2+{\omega }^2$又可以表示为$x^2-{\omega }^4\times {\omega }^2$也就是$x^2-{\omega }^6$

Cooley-Tukey蝴蝶变换

刚刚分解已经讲过怎么去求分支的值，我们可以写一下公式分别看看根的左结点$f_l$和根的右结点$f_r$分别是什么

• $f_l=fmod(x^2-{\omega }^2)=f_{l_0}+f_{l_1}x=(f_0+f_2{\omega }^2)+(f_1+f_3{w}^2)x$
• $f_r=fmod(x^2+{\omega }^2)=f_{r_0}+f_{r_1}x=(f_0-f_2{\omega }^2)+(f_1-f_3{w}^2)x$

有没有发现规律，这两个结点其实每项系数都是对称的$f_{l_0}$对应$(f_0+f_2{\omega }^2)$,$f_{r_0}$对应$(f_0-f_2{\omega }^2)$。
那么我们就可以画出这样的图
$\left[\begin{array}{l}{f}_{l_0}\\ f_{r_0}\end{array}\right]\leftarrow \left[\begin{array}{l}{f}_0+t\\ f_0-t\end{array}\right]$

形如蝴蝶，因此称作蝴蝶变换。其实对于上图树中的每一个节点的多项式环分解成左右两棵子树的过程都是该节点多项式环中的间隔$\frac{l}{2}$(这里$l$为当前待分解多项式环的度数)距离的两组系数做一个Cooley-Tukey蝴蝶变换。
完整来看，这样一颗二叉树（下图1）就对应这样的蝴蝶变换（下图2）

从上图的蝴蝶变换我们可以看出，往往下层的某项系数是由上层的某两项系数决定的，而这两个系数是对称的，只要确定其中一个的位置也就能确定另外一个的位置。例如第一层的第0项系数$f_0$和$f_2$决定了下一层的零项系数（左右都是，不过是对称的）。那么我们就可以固定其中一个即可。

Gentleman-Sande 蝴蝶变换

如果我们已知分解后的左右节点系数值，如何求得根节点的值？这里就引入Gentleman-Sande蝴蝶变换。我们先看下面一个例子（对应树中根节点${\mathbb{Z}}_q/(x^4-w^4)$分解到左子树${\mathbb{Z}}_q/(x^2-w^2)$和右子树${\mathbb{Z}}_q/(x^2+w^2)$过程的逆过程）：

在这个例子中我们已知$f_l$和$f_r$的系数$f_{l0},f_{l1},f_{r0},f_{r1}$，如何求得$f(x)$的系数？
其实也很简单，因为我们求$f_l$和$f_r$的时候，是将$x^2=w^2$和$x^2=-w^2$代入到$f(x)$中，因此我们可以得到$f_0+f_{1}x+f_2\cdot w^2+(f_3\cdot w^2)x=f_{l0}+f_{l1}x$和，$f_0+f_{1}x-f_2\cdot w^2-(f_3\cdot w^2)x=f_{l0}+f_{l1}x$那么可得下列式子:
$$\begin{gathered} f_0+f_2\cdot w^2=f_{l0}(1) \\ {f}_{1}x+f_3\cdot w^2=f_{l1}(2) \\ {f}_0-f_2\cdot w^2=f_{r0}(3) \\ {f}_{1}x-f_3\cdot w^2=f_{r1}(4) \end{gathered}$$
四个式子四个未知数可以求得每个系数的值为
$$\begin{gathered} f_0=\frac{1}{2}(f_{l0}+f_{r_0}) \\ {f}_2=\frac{1}{2}{w}^{-2}(f_{l_0}-f_{r_0}) \\ {f}_1=\frac{1}{2}(f_{l1}+f_{r_1}) \\ {f}_3=\frac{1}{2}{w}^{-2}(f_{l_1}-f_{r_1}) \end{gathered}$$
因此GS变换如下图

如何减少乘法？
可以从一开始就不乘$\frac{1}{2}$，则每层都是原来真实系数的两倍，最终得到系数应该是真实系数的$2^{logn}=n$倍。因此在逆NTT的时候需要乘$\frac{1}{n}$，但这样减少了一半的乘法。

FFT-trick及Cooley-Tukey代码实现

由上面分析可以知道，我们只需要遍历每一层的一半系数就能确定下一层的系数。因此我们可以设置变量$j$用来遍历该层一半的系数，变量$start$确定$j$的开始位置，以及变量$len$确定$j$的一个遍历范围。这里要注意，每层结点都是2的该层的幂次方-1，每个结点也就相当于一个多项式环也是要进行相应分解的，因此这里$start$就对应每个结点多项式环的第一个系数。对于第一层$len=\frac{n}{2}$代表的就是树中根节点的分解过程，对于第二层两个start，第一个start对应的就是树中第二层左儿子的分解，第二个start对应的树中第二层右儿子的分解，依次类推。

那么我们就可以得到如下伪代码：

Forward NTT Alogrithm:
	k = 1
	for len = n/2, len >  0, len = len/2 do
		for start = 0 ; start < n; start = j+len do
			for j = start; j < start + len; j++ do
			 	t = w[k] * f[j+len]
			 	f[j+len] = f[j] - t
			 	f[j] = f[j] + t
			 end for
			 k++
		end for
	end for
			 	

对应C代码就是

void simple_ntt(int16_t r[N])
{
    int len, start, j, k;
    int16_t t, omega;

    k = 1;
    for (len = N / 2; len >= 1; len >>= 1)
    {
        for (start = 0; start < N; start = j + len)
        {
            omega = omegas[k++];
            for (j = start; j < start + len; ++j)
            {
                t = (omega * r[j + len]) % Q;
                r[j + len] = (r[j] - t) % Q;
                r[j] = (r[j] + t) % Q;
            }
        }
    }
}

这里的omega是该层对应的本原单位根

优化NTT过程

因为前面的计算结果都是要回到${\mathbb{Z}}_q$域上的，因此要不断取模，C语言中取模运算(%)用到了除法，是比较耗时且效率低的。有两种约减算法可以改善这种模运算的耗时性。分别是

• Barret约减
• Montgomery约减

Barret约减

计算$r=amodq$，其中$a$为16位有符号整数，$-2^{15}\le r<q$。标准的约减需要一次整数除法，但是可以通过倒数形式表示为乘法。
$r=a\mathrm{m}\mathrm{o}\mathrm{d}q=a-q\cdot \left[\frac{a}{q}\right]=a-q\cdot \lfloor a\cdot q^{-1}\rfloor$
C语言实现

int16_t barrett_reduce(int16_t a) {
  int16_t t;
  const int16_t v = ((1U << 26) + KYBER_Q/2)/KYBER_Q;

  t  = (int32_t)v*a >> 26;
  t *= KYBER_Q;
  return a - t;
}

Montgomery约减

首先会选择一个基数$\beta$，使得$\beta$与$q$互素，且$\beta >q$，一般会取2的整数倍，当满足$0\le a<\beta a$时，将计算$r=amodq$转换为先计算$r^{{}^{\prime }}=a{\beta }^{-1}modq$再计算$r=r^{{}^{\prime }}\beta modq$.计算$r^{{}^{\prime }}$用到的算法就是蒙哥马利约减算法。
C语言代码如下

int16_t montgomery_reduce(int32_t a)
{
  int32_t t;
  int16_t u;

  u = a*QINV;
  t = (int32_t)u*KYBER_Q;
  t = a - t;
  t >>= 16;
  return t;
}