Recursive NTT for Kyber

Question

since the NTT is a variant of the DFT respectively FFT, it seems possible to implement a recursive algorithm for the NTT of Kyber with a divide-and-conquer concept. This is the reason for my question here, I would be interested in how to recursively realize the NTT from the current Kyber specification which is defined as follows:

$$ NTT(f) = \hat{f} = (\hat{f}_0+\hat{f}_1 X, \hat{f}_2 + \hat{f}_3X,...,\hat{f}_{254}+\hat{f}_{255}X)$$ with $$\hat{f}_{2i} = \sum_{j=0}^{127} f_{2j} \xi^{(2 br_7(i) +1) j}$$ $$\hat{f}_{2i+1} = \sum_{j=0}^{127} f_{2j+1} \xi^{(2 br_7(i) +1) j}$$

So the question is: How do we obtain from this formulas a recursive NTT.

---

Ideas/Reasoning

For a recursive implementation of the formula for $\hat{f}_{2i}$, I have considered the following:

$$\hat{f}_{2i} = \sum_{j=0}^{127} f_{2j} \cdot \xi^{(2 br_{7}(i)+1)j} = \sum_{j=0}^{63} f_{4j} \xi^{(2br_{7}(i)+1)2j} + \sum_{j=0}^{63} f_{4j+2} \xi^{(2br_{7}(i)+1)(2j+1)} $$ $$= \sum_{j=0}^{63} f_{4j} \xi^{(2br_{7}(i)+1)2j} + \xi^{2br_{7}(i)+1}\sum_{j=0}^{63} f_{4j+2} \xi^{(2br_{7}(i)+1)2j}$$

also applies due to the periodicity:

$$\hat{f}_{2i + N/2} = \sum_{j=0}^{63} f_{4j} \xi^{(2br_{7}(i)+1)2j} - \xi^{2br_{7}(i)+1}\sum_{j=0}^{63} f_{4j+2} \xi^{(2br_{7}(i)+1)2j}$$

So now we have two separate sums that also halve the input size in each step. This is very similar to a divide-and-conquer approach. Overall, this is very similar to the FFT.

Inspired by the FFT algorithm, a recursive version could look something like this (I think):

$\textbf{recursive_ntt}(a, \xi):$
$\quad n \leftarrow \textbf{length}(a)$
$\quad \textbf{if} \ n = 1$
$\quad \quad \textbf{return} \ a$
$\quad a^\textbf{even} \leftarrow (a_0,a_2,...,a_{n-2})$
$\quad a^\textbf{odd} \leftarrow (a_1,a_3,...,a_{n-1})$
$\quad y^\textbf{even} \leftarrow \textbf{recursive_ntt}(a^\textbf{even},\xi^2)$
$\quad y^\textbf{odd} \leftarrow \textbf{recursive_ntt}(a^\textbf{odd},\xi^2)$

$\quad \textbf{for} \ k \leftarrow 0 \ \textbf{to} \ n/2 -1 \ \textbf{do}$
$\quad \quad y_k = y^\textbf{even}_k + \xi^{2 \cdot br_7(i)+1} \cdot y^\textbf{odd}_k$
$\quad \quad y_{k + n/2} = y^\textbf{even}_k - \xi^{2 \cdot br_7(i)+1} \cdot y^\textbf{odd}_k$

$\quad \textbf{return} \ y$

Here I would like to note that the input $a$ is an array that contains the even indices and is already reduced to half size. The idea is $f$ consists of all values. For $\hat{f}_{2i}$ we only need the even coefficients $f_{2i}$. If we split the even coefficients further into even and "odd", this corresponds roughly to what formulas 3 and 4 express. So if $f = (f_0, f_1, f_2, f_3, f_4,...f_{255})$ then $a = (f'_0=f_0, f'_1 = f_2, f'_2 = f_4, f'_3=f_6, ..., f'_{127}=f_{254})$, so $f_{4j}$ becomes $f'_{2j'}$ and $f_{4j+2}$ becomes $f'_{2j'+1}$, where $j'=0,1,...,64$.

---

Question(s)

1. It is still an interesting and open question how we could define a recursive formalization of the NTT formulas of the specification. I made a try, but I'am not sure.

2. Is the bit-reverse operation in a recursive call needed? This is still an open question. The FFT for itself does not need a bit-reversal in it's recursive description, even so the recursive NTT? Is it possible to ommit the bit-reverse operation in the pseudo-code from above? Do we need the $br_7(i)$ operation in a recursive implementation or is it only required if we work iterative in-place?

3. I have implemented $\hat{f}_{2i}$ as it is given in the specification and also recursively according to my description. The outputs are completely different. But if I leave out the bit-reverse in both implementations, the output is the same, how can that be?

---

Answer 1

I don't fully understand your confusion, but the following might help. The NTT for Kyber is not the standard NTT (which has a standard recursive algorithm, by analogy with the FFT, as you mentioned). Call this standard NTT $\mathcal{F}_{\mathsf{standard}}$, and Kyber's NTT $\mathcal{F}_{\mathsf{Kyber}}$. These two are related via

$$\mathcal{F}_{\mathsf{Kyber}}(x) = \sigma^{-1}(\mathcal{F}_{\mathsf{standard}}(x))$$

where $\sigma$ is a permutation (involution even, e.g. $\sigma(\sigma(x)) = x$) often called "Bit reversal".

So, if you want to write $\mathcal{F}_{\mathsf{Kyber}}(x)$ recursively, you could do the "obvious" recursive algorithm for $\mathcal{F}_{\mathsf{standard}}$, and then post-processes it by the "bit reversal" permutation. This requires additional computation ($O(n)$ random accesses to compute the permutation) which is asymptotically smaller than the $O(n\log n)$ cost of either $\mathcal{F}$, e.g. asymptotically it doesn't matter, so likely is good enough for you.

In other words, if your only issue with verifying your algorithm to Kyber's spec is that the coefficients $\widehat{f}_i$ are in the wrong order, this is expected. Kyber uses this (mathematically less elegant) order as it works better with AVX implementations, see page 5 of this for a quick comment stating that.