Number Theoretic Transform (NTT) for ML-KEM, Step by Step

Key points

ML-KEM spends most of its time multiplying polynomials in $\mathbb{Z}_q[X]/(X^{256}+1)$ . Schoolbook multiplication needs 65,536 coefficient products; the NTT cuts this to a few thousand.
The NTT is a Fourier transform over integers mod $q$ . It works because $q$ has the right roots of unity: $\zeta = 17$ has order 256 mod 3329.
Mod 3329 there is no 512th root of unity, so ML-KEM’s NTT is incomplete. It uses 7 layers, ends with 128 degree-1 pieces, and multiplies them with a small base-case formula.
The forward transform uses Cooley–Tukey butterflies, and the inverse uses Gentleman–Sande butterflies plus a final scaling by $128^{-1} = 3303$ .
Most NTT bugs come from twiddle order, missing scaling, or modular reduction. Test against schoolbook multiplication.

Why ML-KEM needs a fast multiplication

In ML-KEM every public-key and ciphertext element is a polynomial with 256 coefficients mod $q = 3329$ . Key generation computes $t = As + e$ , where $A$ is a $k \times k$ matrix of such polynomials. That is 9 polynomial products for ML-KEM-768 before any encryption happens. These are the Module-LWE equations, and their cost is dominated by polynomial multiplication.

Multiplying two degree-255 polynomials coefficient by coefficient takes $256^2 = 65{,}536$ products. The NTT turns the multiplication into cheap pointwise work after a transform that costs about $n \log n$ operations. The same idea powers the FFT for real numbers, but everything stays exact because all arithmetic is modulo $q$ .

Negacyclic convolution: what the ring does to products

The ring is $R_q = \mathbb{Z}_q[X]/(X^n + 1)$ . Because $X^n \equiv -1$ , any term that overflows past degree $n-1$ wraps around to the bottom with its sign flipped:

c_k = \sum_{i+j=k} a_i b_j \; - \sum_{i+j=k+n} a_i b_j \pmod q

This is called a negacyclic convolution. A plain FFT computes a cyclic convolution, where the wrapped terms keep their sign. The negacyclic sign comes for free if we evaluate at roots of $X^n + 1$ rather than roots of $X^n - 1$ .

Intuition

A polynomial is fully determined by its values at enough points. Multiplying two polynomials is hard coefficient by coefficient, but easy value by value: $(a \cdot b)(r) = a(r) \cdot b(r)$ . The NTT evaluates at cleverly chosen points (the roots of $X^n+1$ mod $q$ ), you multiply the values, and the inverse NTT turns them back into coefficients. The cleverness is in choosing points so that evaluating at all of them at once reuses almost all the work.

Roots of unity mod q

An element $\omega \in \mathbb{Z}_q$ is a primitive $m$ -th root of unity if $\omega^m \equiv 1$ and no smaller positive power equals 1. Such an element exists exactly when $m$ divides $q-1$ . To split $X^n + 1$ into linear factors you need a primitive $2n$ -th root $\psi$ , because the roots of $X^n+1$ are the odd powers $\psi, \psi^3, \dots, \psi^{2n-1}$ .

ML-KEM’s modulus

For ML-KEM, $q - 1 = 3328 = 2^8 \cdot 13$ . So a primitive 256th root exists, and FIPS 203 fixes $\zeta = 17$ ( $17^{128} \equiv -1$ , $17^{256} \equiv 1 \pmod{3329}$ ). A 512th root does not exist, so the best available factorization is into 128 quadratics:

$X^{256} + 1 = \prod_{i=0}^{127} \left(X^2 - \zeta^{\,2\,\mathrm{BitRev}_7(i) + 1}\right) \pmod{3329}$
By the Chinese remainder theorem, $R_q \cong \prod_{i=0}^{127} \mathbb{Z}_q[X]/(X^2 - \gamma_i)$ . The NTT of a polynomial is its list of 128 remainders, each a degree-1 polynomial stored as two coefficients.

A worked example you can check by hand

Shrink everything: $n = 4$ , $q = 17$ . Here $\psi = 2$ is a primitive 8th root of unity (powers: 1, 2, 4, 8, 16, 15, 13, 9). Multiply

a = 1 + X, \qquad b = 1 + X^3

Schoolbook first: $ab = 1 + X + X^3 + X^4$ . Since $X^4 \equiv -1$ , the $X^4$ term cancels the constant and $c = X + X^3$ .

Full NTT: evaluate at the four roots

The roots of $X^4+1$ mod 17 are the odd powers of 2: 2, 8, 15 and 9. Evaluate and multiply:

root r	a(r)	b(r)	a(r)·b(r) mod 17	c(r) = r + r³
2	3	9	10	10
8	9	3	10	10
15	16	10	7	7
9	10	16	7	7

The pointwise products match $c$ at every root, so interpolating them back gives exactly $X + X^3$ . That is the whole idea; the rest is doing it quickly.

Incomplete NTT: stop one layer early, like ML-KEM

Pretend 17 had no 8th root, only the 4th root $\zeta = \psi^2 = 4$ ( $4^2 = 16 \equiv -1$ ). Then $X^4 + 1 = (X^2 - 4)(X^2 + 4)$ and one butterfly layer computes both remainders:

f \bmod (X^2 \mp 4) = (f_0 \pm 4 f_2) + (f_1 \pm 4 f_3)\,X

For $b = (1, 0, 0, 1)$ this gives $1 + 4X$ and $1 - 4X = 1 + 13X$ . For $a$ both remainders are $1 + X$ . Now multiply each pair modulo its own quadratic with the base-case rule, where $X^2 = \gamma$ :

(a_0 + a_1 X)(b_0 + b_1 X) \bmod (X^2 - \gamma) = (a_0 b_0 + a_1 b_1 \gamma) + (a_0 b_1 + a_1 b_0)\,X

Mod $X^2 - 4$ : $(1 + X)(1 + 4X) = (1 + 4\cdot 4) + 5X = 17 + 5X \equiv 5X$ . Check: $c = X + X^3 \equiv X + 4X = 5X$ .
Mod $X^2 + 4$ ( $\gamma = 13$ ): $(1 + X)(1 + 13X) = (1 + 13\cdot 13) + 14X \equiv 14X$ . Check: $X + 13X = 14X$ .

ML-KEM does exactly this with 7 layers instead of 1, ending in 128 quadratics $X^2 - \gamma_i$ . FIPS 203 calls the pairwise step BaseCaseMultiply (Algorithm 12) and applies it through MultiplyNTTs (Algorithm 11).

Butterflies: how the layers work

Each layer splits every factor $X^{2m} - \zeta^2$ into $(X^m - \zeta)(X^m + \zeta)$ . For a block $f = f_{\text{lo}} + X^m f_{\text{hi}}$ , the two remainders are $f_{\text{lo}} \pm \zeta f_{\text{hi}}$ . Computed in place, that is the Cooley–Tukey butterfly:

t = \zeta \cdot f[j+m], \qquad f[j+m] = f[j] - t, \qquad f[j] = f[j] + t

ML-KEM’s forward NTT (FIPS 203 Algorithm 9) runs this with block half-lengths 128, 64, …, 2. That is 7 layers of 128 butterflies, 896 multiplications in total. The inverse (Algorithm 10) runs the layers in reverse with the Gentleman–Sande butterfly, which undoes one split:

t = f[j], \qquad f[j] = t + f[j+m], \qquad f[j+m] = \zeta \cdot (f[j+m] - t)

Each inverse layer doubles the values, so the result is finally multiplied by $128^{-1} \bmod 3329 = 3303$ .

Why the zetas are in bit-reversed order

An in-place Cooley–Tukey transform leaves its outputs in bit-reversed order. Instead of shuffling the array, FIPS 203 keeps the output in that order and indexes the twiddle factors by $\zeta^{\mathrm{BitRev}_7(k)}$ , where $\mathrm{BitRev}_7$ reverses the 7-bit binary expansion of $k$ . The first entries are 1, 1729, 2580, 3289, …. The NTT output is only ever consumed by MultiplyNTTs and the inverse NTT, so the permuted order is harmless as long as all three agree. ML-KEM even samples the matrix $A$ directly in the NTT domain, so it is never transformed at all.

Reference implementation in Python

This follows FIPS 203 Algorithms 9–12 line for line. It is written for clarity, not speed or constant time, and is useful as a test oracle for an optimized implementation.

Q = 3329
F = 3303  # 128^-1 mod Q

def bitrev7(i):
    return int(f"{i:07b}"[::-1], 2)

ZETAS = [pow(17, bitrev7(i), Q) for i in range(128)]            # FIPS 203 Alg. 9/10
GAMMAS = [pow(17, 2 * bitrev7(i) + 1, Q) for i in range(128)]   # FIPS 203 Alg. 11

def ntt(f):
    f, k, length = list(f), 1, 128
    while length >= 2:
        for start in range(0, 256, 2 * length):
            zeta = ZETAS[k]; k += 1
            for j in range(start, start + length):
                t = zeta * f[j + length] % Q          # Cooley-Tukey butterfly
                f[j + length] = (f[j] - t) % Q
                f[j] = (f[j] + t) % Q
        length //= 2
    return f

def intt(f):
    f, k, length = list(f), 127, 2
    while length <= 128:
        for start in range(0, 256, 2 * length):
            zeta = ZETAS[k]; k -= 1
            for j in range(start, start + length):
                t = f[j]                              # Gentleman-Sande butterfly
                f[j] = (t + f[j + length]) % Q
                f[j + length] = zeta * (f[j + length] - t) % Q
        length *= 2
    return [x * F % Q for x in f]

def multiply_ntts(a, b):
    c = [0] * 256
    for i in range(128):
        a0, a1, b0, b1 = a[2*i], a[2*i+1], b[2*i], b[2*i+1]
        c[2*i]     = (a0 * b0 + a1 * b1 * GAMMAS[i]) % Q   # mod X^2 - gamma
        c[2*i + 1] = (a0 * b1 + a1 * b0) % Q
    return c

# product of a and b in Z_q[X]/(X^256 + 1):
# intt(multiply_ntts(ntt(a), ntt(b)))

Test it the obvious way: for random $a$ and $b$ , check that $\mathrm{intt}(\mathrm{ntt}(a)) = a$ and that the NTT product equals a schoolbook negacyclic convolution. Both hold for the code above.

Montgomery and Barrett reduction

Production code never calls % on secret data, because hardware division can be slow and variable-time. Two reductions do the job with multiplications and shifts:

Montgomery

Keep values as $aR \bmod q$ with $R = 2^{16}$ . A product $x < q \cdot 2^{15}$ is reduced to $xR^{-1} \bmod q$ using $q^{-1} \bmod 2^{16} = -3327$ . The reference Kyber code stores its zetas in Montgomery form so each butterfly needs one Montgomery multiplication.

Barrett

Approximate $\lfloor x/q \rfloor$ as $(x \cdot v) \gg 26$ with $v = \lfloor (2^{26} + \lfloor q/2 \rfloor)/q \rfloor = 20159$ , then subtract that multiple of $q$ . It is used to pull coefficients back into range between layers, where lazy additions let values grow.

With Montgomery-form data, the final $128^{-1}$ of the inverse NTT is merged with the $R$ factors into one constant. That is why the reference code uses a different number there than FIPS 203’s 3303.

Common implementation bugs

Twiddles in natural order. Using $\zeta^k$ instead of $\zeta^{\mathrm{BitRev}_7(k)}$ gives a transform that round-trips but multiplies wrongly.
Wrong $\gamma$ in the base case. The quadratic factors use odd powers, $\zeta^{2\,\mathrm{BitRev}_7(i)+1}$ , and alternate in sign between neighbours.
Missing or doubled scaling. Forgetting 3303, or applying it on top of a Montgomery constant that already includes it, scales the whole product.
Overflow in lazy reduction. Skipping reductions between layers is a standard optimization. Bound the growth per layer, or 16-bit lanes overflow silently.
Mixing domains. Adding a normal-domain polynomial to an NTT-domain one compiles fine and produces garbage. Keep the domain in the type.
Only testing round-trips. $\mathrm{intt}(\mathrm{ntt}(a)) = a$ passes for several of the bugs above. Always compare products with a schoolbook oracle and the official test vectors.

Constant-time considerations

The NTT itself has a fixed memory-access pattern that does not depend on data, which makes it easy to run in constant time. The risks are in the arithmetic around it: % or / on secrets, conditional subtractions written as branches, and compilers turning branch-free code back into branches. The KyberSlash timing issues (2023–2024) came from a division by $q$ on secret-dependent values in the compression code, not the NTT, which shows where to look. Use Barrett or Montgomery reduction everywhere, and check the compiled output with a tool such as dudect or valgrind-based checks for secret-dependent branching.

Comparison: the NTT in ML-DSA

ML-DSA (Dilithium) uses the same ring shape with $q = 8380417$ . There $q - 1 = 2^{13} \cdot 3 \cdot 11 \cdot 31$ , so a primitive 512th root exists (FIPS 204 uses $r = 1753$ ). ML-DSA therefore runs a complete 8-layer NTT down to 256 linear factors and multiplies purely pointwise, with no base case. Its coefficients need 32-bit arithmetic, while ML-KEM’s fit in 16 bits. The small modulus is the price of that compactness, and it is why ML-KEM’s transform stops one layer early.

Check Your Understanding

Q1. Why is the ML-KEM NTT “incomplete”, stopping at 128 degree-1 pieces instead of 256 single values?
Q2. In $\mathbb{Z}_q[X]/(X^4 + 1)$ , what does the term $X^3 \cdot X^2$ become?
Q3. You multiply two polynomials by calling NTT on both, multiplying pointwise and calling the inverse NTT, but the result is wrong by a constant factor. What is the most likely bug?
Q4. Why does Dilithium / ML-DSA get a full 8-layer NTT while ML-KEM does not?

3.5 The NTT, step by step