参考文献：

1. [CT65] Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics of computation, 1965, 19(90): 297-301.
2. [Shoup95] Shoup V. A new polynomial factorization algorithm and its implementation[J]. Journal of Symbolic Computation, 1995, 20(4): 363-397.
3. [CHKKS18] Cheon J H, Han K, Kim A, et al. Bootstrapping for approximate homomorphic encryption[C]//Advances in Cryptology–EUROCRYPT 2018: 37th Annual International Conference on the Theory and Applications of Cryptographic Techniques, Tel Aviv, Israel, April 29-May 3, 2018 Proceedings, Part I 37. Springer International Publishing, 2018: 360-384.
4. [CHH18] Cheon J H, Han K, Hhan M. Faster homomorphic discrete fourier transforms and improved fhe bootstrapping[J]. Cryptology ePrint Archive, 2018.
5. Nussbaumer Transform 以及 Amortized FHEW bootstrapping
6. Chimera：混合的 RLWE-FHE 方案
7. 快速乘法技巧：Karatsuba, Toom, Good, Schonhage, Strassen, Nussbaumer
8. Paterson-Stockmeyer 多项式求值算法

Baby-Step Giant-Step

Shoup95

文章 [Shoup95] 研究并实现了 BSGS factoring method，用于将单变元多项式分解为不可约因子。其中使用了 CRT 和 FFT 来表示多项式（比 GHS12 的 Doube-CRT 更早） ，并且实现了多项式的快速乘法、除法、逆元、平方、GCD 等等运算。

多项式分解可以分为三步，

1. square-free factorization：将多项式分解为 $f={\prod }_i{f}_i$，其中的 $f_i$ 是平方自由的
2. distinct-degree factorization：将平方自由多项式分解为 $f_i={\prod }_j{f}_{i,j}$，其中的 $f_{i,j}$ 是一些度数都为 $j$ 的不可约因式的乘积
3. equal-degree factorization：将不可约因式都是相同度数的平方自由多项式 $f_{i,j}$，分解为这些不可约多项式

主要步骤集中在 step 2，[Shoup95] 观察到事实：对于任意的非负整数 $a,b\in {\mathbb{Z}}^{+}$，多项式 $h_{a,b}(x)=x^{p^a}-x^{p^b}\in GF(p)[x]$ 以所有的满足 $\mathrm{deg}f|(a-b)$ 的不可约多项式 $f$ 为因式。

对于 $\mathrm{deg}f\le n$ 的平方自由多项式，它的真因子度数不超过 $n/2$，我们令 $f_d,1\le d\le n$ 是它的全部 $d$ 次不可约因子的乘积。我们可以枚举全部的 $1\le a-b\le n$，计算出 $h_{a,b}(x)$，再计算 $\gcd (h_{a,b},f)$ 从而获得这些 $f_d$

[Shoup95] 使用 BSGS 算法来计算这些 $h_{a,b}$，设置真因子的度数上界 $B$，将它分为 $B=l\cdot m$，Baby-Step 就是 { i : 1 ≤ i ≤ l } \{i:1 \le i \le l\} {i:1≤i≤l}，Giant-Step 就是 { l ⋅ j : 1 ≤ j ≤ m } \{l \cdot j:1 \le j \le m\} {l⋅j:1≤j≤m}，

然而，如果简单地直接计算 $h_i,H_j$，上述算法依旧是不实用的。[Shoup95] 以迭代的方式计算它们：$h_{i+1}=h_i(h_1)(\mathrm{mod}f)$，$H_{j+1}=H_j(H_1)(\mathrm{mod}f)$，现在的问题是如何快速计算这些 modular-composition，形如 $g(h)(\mathrm{mod}f)$

[Shoup95] 依旧采取 BSGS 算法（类似于 [PS73] 多项式求值算法），选取参数 $t\approx \sqrt{n}$，预计算 $h^i(\mathrm{mod}f),0\le i\le t$ 表格，那么：
$$g(x)=\sum _{j=0}^{n/t}{g}_j(x)\cdot y^j,y=x^t,\mathrm{deg}{g}_j<t$$
于是，直接使用预计算表的内容，简单计算加法（以及数乘），
$$g_j(h)=\sum _{i=0}^t{g}_{j,i}\cdot h^i(\mathrm{mod}f)$$
接着，采取 Horner 法则，计算出
$$g(x)=((g_{n/t}\cdot h^t+\cdots )\cdot h^t+g_1)\cdot h^t+g_0$$
其中的多项式运算都是以 Double-CRT 方式计算的，总的复杂度为 $O(n^{2.5}+n\log n\log \log n\log p)$

CT65

[CT65] 给出了递归形式的 DFT 分解，其实也是可以视为一种 BSGS 版本的 FFT 算法。DFT 公式为：
$$A_j:=\sum _{k=0}^{N-1}{a}_k\cdot {\zeta }^{jk}$$
采取 BSGS 算法，分解为 $N=N_1\cdot N_2$，设置索引
$$\begin{array}{rl}j& :=N_1\cdot j_1+j_0,j_0\in [N_1],j_1\in [N_2]\\ k& :=N_2\cdot k_1+k_0,k_0\in [N_2],k_1\in [N_1]\end{array}$$
那么有
$$\begin{array}{rl}{A}_{j_1,j_0}& :=\sum _{k_0\in [N_2]}\sum _{k_1\in [N_1]}{a}_{k_1,k_0}\cdot {\zeta }^{jk}\\ & =\sum _{k_0\in [N_2]}\left(\sum _{k_1\in [N_1]}{a}_{k_1,k_0}\cdot {\zeta }^{N_{2}j{k}_1}\right)\cdot {\zeta }^{j{k}_0}\\ & =\sum _{k_0\in [N_2]}\left(\sum _{k_1\in [N_1]}{a}_{k_1,k_0}\cdot {\zeta }^{N_{2}j{k}_1}\right)\cdot {\zeta }^{j_0{k}_0}\cdot {\zeta }^{N_1{j}_1{k}_0}\end{array}$$
于是，将 $a_N$ 按照行主序，排列为形状 $N_1\times N_2$ 的矩阵 $a_{N_1\times N_2}$，

1. 对于每一个 $k_0$，利用形状 $N_1\times N_1$ 的矩阵
   W 1 : = { ζ j 0 k 1 } j 0 , k 1 W_1:=\{\zeta^{j_0k_1}\}_{j_0,k_1} W1​:={ζj0​k1​}j0​,k1​​
   计算长度为 $N_1$ 的各个列矢 $a_{k_0}$ 的 NTT 变换（单位根为 { ζ N 1 j 0 , j 0 ∈ [ N 1 ] } \{\zeta_{N_1}^{j_0},j_0 \in [N_1]\} {ζN1​j0​​,j0​∈[N1​]}），得到形状 $N_1\times N_2$ 的矩阵
   $$W_1\times a_{N_1\times N_2}={\left\{A_{j_0,k_0}^{\prime }:=\sum _{k_1\in [N_1]}{a}_{k_1,k_0}\cdot {\zeta }^{N_{2}j{k}_1}\right\}}_{j_0,k_0}$$

2. 利用形状 $N_1\times N_2$ 的矩阵
   W 2 : = { ζ j 0 k 0 } j 0 , k 0 W_2 := \{\zeta^{j_0k_0}\}_{j_0,k_0} W2​:={ζj0​k0​}j0​,k0​​
   对它做 Hadamard 乘积，扭曲矩阵 $A^{\prime }$ 使得接下来的运算是标准 NTT（否则就需要恰当的扭曲后续 NTT 采用的单位根），此时的结果是形状 $N_1\times N_2$ 的矩阵
   $$W_2\odot A_{N_1\times N_2}^{\prime }={\left\{A_{j_0,k_0}^{\prime \prime }:={\zeta }^{j_0{k}_0}\cdot \sum _{k_1\in [N_1]}{a}_{k_1,k_0}\cdot {\zeta }^{N_{2}j{k}_1}\right\}}_{j_0,k_0}$$

3. 对于每一个 $j_1$，利用形状 $N_2\times N_2$ 的矩阵
   W 3 : = { ζ N 2 j 1 k 0 } j 1 , k 0 W_3 := \{\zeta_{N_2}^{j_1k_0}\}_{j_1,k_0} W3​:={ζN2​j1​k0​​}j1​,k0​​
   计算长度为 $N_2$ 的各个行矢 $A_{j_0}^{\prime \prime }$ 的 NTT 变换（单位根为 { ζ N 2 j 1 , j 1 ∈ [ N 2 ] } \{\zeta_{N_2}^{j_1},j_1 \in [N_2]\} {ζN2​j1​​,j1​∈[N2​]}），得到形状 $N_2\times N_1$ 的矩阵
   W 3 × ( A N 1 × N 2 ′ ′ ) T = { A j 1 , j 0 } j 1 , j 0 W_3 \times (A_{N_1 \times N_2}'')^T = \{A_{j_1,j_0}\}_{j_1,j_0} W3​×(AN1​×N2​′′​)T={Aj1​,j0​​}j1​,j0​​

对于形状 $N_2\times N_1$ 的矩阵 $A_{N_2\times N_1}$，按照行主序读取为 $A_N=NTT(a_N)$

总之，$a_N,A_N$ 都按照行主序排列为矩阵（形状不同），那么有：
$$A_{N_2\times N_1}=W_3\times (W_2\odot (W_1\times a_{N_1\times N_2}){)}^T$$
其实，这个过程可以用 Nussbaumer Transform 的环同态表示
$$\mathbb{F}[x]/(x^N-1)\cong (\mathbb{F}[y]/(y^{N_1}-1))[x]/(x^{N_2}-y)\cong (\mathbb{F}[y]/(y^{N_1}-1))[z]/(z^{N_2}-1)$$

CHKKS18

[CHKKS18] 利用 SIMD 技术的 Hadamard 和 Rotate 运算，实现同态的线性运算，它也采取了 BSGS 技巧。 我们默认 index 都是自动 $(\mathrm{mod}n)$ 的，基本符号：

• 对于任意的线性变换 $M\in {\mathbb{C}}^{n\times n}$，简记 $dia{g}_i(M)=[M_{0,i},M_{1,i+1},\cdots ,M_{n,i+n}]$ 是第 $i\in \mathbb{Z}$ 条对角线（可以是负数 $-i$，就是第 $n-i$ 条对角线）
• 对于任意的矢量 $v\in {\mathbb{C}}^n$，简记 $ro{t}_i(v)=[v_i,v_{i+1},\cdots ,v_{i+n-1}]$ 是循环左移 $i\in \mathbb{Z}$ 距离（可以是负数 $-i$，循环右移 $i$ 距离）

采取 BSGS 算法，分解 $n=l\times k$，线性变换可以表示为：

最优化时选取 $k\approx \sqrt{n}$，计算复杂度为：$O(\sqrt{n})$ 次 Rotate 运算（关于 $v$ 的密文），$O(n)$ 次 Hadamard 运算。对于公开的固定矩阵 $M$，其中的 $ro{t}_{-ki}(dia{g}_{ki+j}(M))$ 是预计算的常数多项式（用 InvDFT 编码），它不需要 CKKS 密文下的 Rotate 运算。

[CHKKS18] 将上述的线性变换转换到 slot-packing CKKS 下同态计算，并利用它来实现 coeff-to-slot，从而批处理 CKKS 自举。用到的线性变换是 DFT 和 InvDFT，[CHKKS18] 将它们视为通用的线性变换，利用这个同态矩阵乘法来实现。

不过，对于公开的线性变换，直接使用 TFHE 提出的那种 Functional Key-Switch，效率要高得多。对于秘密的线性变换，TFHE 也提出了根据 $M$ 来构造 KS-Key for M，从而支持 Private Functional Key-Switch。但是如果没有提供这种特殊的 KS-Key for M，而是将 $M$ 加密为一般的 CKKS 密文，那么就只能使用上面的同态矩阵乘法，利用 Rotate 和 Hadamard 慢慢计算。

Faster Homomorphic DFT

[CHH18] 观察到 DFT 矩阵拥有稀疏分解（也就是蝴蝶算法），因此对于这种特殊的线性变换，可以比 [CHKKS18] 的通用矩阵乘法，复杂度降低一个 $n$ 因子。它可以应用在 [CHKKS18] 的 CKKS 批自举上，计算速度提高了数百倍。

Sparse-Diagonal matrix Factorization

[CHH18] 它说根据 [CT65] 的递归式 FFT，可以将 DFT 矩阵做如下的稀疏分解：

继续迭代地分解前者，最终可以获得：

容易看出， d i a g i ( D 2 i ( n ) ) ≠ 0 ⃗ ⟺ k ∈ { 0 , ± n 2 i } diag_i(D_{2^i}^{(n)}) \neq \vec 0 \iff k\in \{0,\pm \dfrac{n}{2^i}\} diagi​(D2i(n)​)=0 ⟺k∈{0,±2in​}，只有 3 条对角线是非零的，因此采取斜线乘法来计算是十分高效的，
D 2 i ( n ) ⋅ v = ∑ k ∈ { 0 , ± n / 2 i } d i a g k ( D 2 i ( n ) ) ⊙ r o t k ( v ) D_{2^i}^{(n)} \cdot v = \sum_{k\in \{0,\pm n/2^i\}} diag_k(D_{2^i}^{(n)}) \odot rot_{k}(v) D2i(n)​⋅v=k∈{0,±n/2i}∑​diagk​(D2i(n)​)⊙rotk​(v)
算法为：

注意到 $ro{t}_0(v)=v$ 不必计算，对于特殊情况 $i=1$ 根据 $ro{t}_{n/2^i}(v)=ro{t}_{-n/2^i}(v)$ 可节约计算。最终，$DFT\cdot v={\prod }_{i=0}^{{\log }_{2}n}{D}_{2^i}^{(n)}\cdot v$ 的复杂度为 $O({\log }_{2}n)$

对于逆变换，由于 DFT 的逆矩阵恰好是其厄米，

因此很明显 CT 蝴蝶和上述的 GS 蝴蝶一样，也是稀疏对角的，从而也存在类似的快速矩阵乘法。

Radix-r

不过虽然上述的操作次数很小，但是计算深度为 ${\log }_2(n)$，需要多层 Rotate 和 CMult 串行，可能导致噪声控制的问题。我们可以合并某些连续的 $k$ 个矩阵，那么深度就降低为 ${\log }_{r}n,r=2^k$，代价是计算复杂度的上升。

根据对角阵的乘法性质：两个对角阵的乘积，依旧是对角阵，
$$dia{g}_i(a)\cdot dia{g}_j(b)=dia{g}_{i+j}(a\odot ro{t}_i(b))$$
可以证明，连续 $k$ 个矩阵的合并
$$D_{k,s}=D_{2^{s+k}}^{(n)}\cdots D_{2^{s+2}}^{(n)}\cdot D_{2^{s+1}}^{(n)}$$
的非零对角线的索引为
$$e_1\cdot \frac{n}{2^{s+1}}+e_2\cdot \frac{n}{2^{s+2}}+\cdots +e_t\cdot \frac{n}{2^{s+k}}$$
其中 e i ∈ { 0 , ± 1 } e_i \in \{0,\pm1\} ei​∈{0,±1}，易知这些 index 都是 $\frac{n}{2^{s+k}}$ 的倍数，绝对值上界为 $\frac{(2^k-1)n}{2^{s+k}}$，这些 index 的个数至多为 $2^{k+1}-1$

此时的 DFT 的复杂度为：$O(r{\log }_{r}n)$ 次 Rotate 和 Hadamard，深度为 $O({\log }_{r}n)$

Hybrid method

由于上述分解的非零 index 呈现算术级数（都是 $\frac{n}{2^{s+k}}$ 的倍数），因此可以采取 BSGS 技巧，提取出公共的计算部分，进一步提高计算效率。

最优化设置 $k_2\approx \sqrt{t}$，此时 DFT 的复杂度为：依旧 $O(r{\log }_{r}n)$ 次 Hadamard，但是只需 $O(\sqrt{r}{\log }_{r}n)$ 次 Ratate（但 $r$ 是常数），深度也还是 $O({\log }_{r}n)$

Result

同态 DFT 的执行时间：秒的量级

CKKS 自举时间：2 分钟延迟，均摊 4 毫秒