Number Theoretic Transform, in depth: roots of unity mod p, a hand-traced example, tested code, three-prime CRT and the negacyclic NTT in ML-KEM

Roots of unity in a finite field

The FFT works because the complex number w = e^(2πi/n) has three properties: w^n = 1, no smaller positive power of w equals 1, and w^(n/2) = −1. The even/odd split, the butterfly and the inverse transform use nothing else. So any number system with an element having those properties supports the same algorithm.

In the integers mod a prime p, the nonzero elements form a cyclic group of order p − 1. There is a primitive root g whose powers run through every nonzero residue. An element of order exactly n exists if and only if n divides p − 1, and it is w = g^((p−1)/n). Because p is prime, the only square roots of 1 are 1 and −1, so w^(n/2), which squares to 1 but is not 1, must be −1. All three properties hold.

Since the radix-2 algorithm needs n to be a power of two, you want primes where p − 1 has a large power of two as a factor: p = c · 2^k + 1. The exponent k caps the transform length at 2^k. The inverse transform also needs n^−1 mod p, which exists because n is a power of two and p is odd. Computing inverses and powers mod p relies on fast exponentiation, covered in Modular Exponentiation.

Every row was checked: each p is prime, has the stated form, and has 3 as its smallest primitive root. 998244353 is the usual default because it is just under 2^30, so the product of two residues fits in a 64-bit integer with room to spare.

A worked example by hand

Take a tiny field so you can check every number: p = 17, n = 4. Since 16 = 4 · 4, a fourth root of unity exists, and w = 4 works: 4^2 = 16 = −1 and 4^4 = 1. The four evaluation points are w^0..w^3 = 1, 4, 16, 13. Multiply f(x) = 1 + 2x by g(x) = 3 + 4x. The product has degree 2, which needs three coefficients, so length 4 is enough and the inputs are padded to [1, 2, 0, 0] and [3, 4, 0, 0].

Now invert. The inverse transform is the forward transform with w^−1 = 13 (since 4 · 13 = 52 = 1 mod 17), followed by multiplication by n^−1 = 4^−1 = 13. For the constant term: 13 · (4 + 1 + 1 + 6) = 156 = 3 mod 17. The other three come out as 10, 8 and 0, so the product is 3 + 10x + 8x^2, which is exactly (1 + 2x)(3 + 4x). No rounding, and the trailing 0 confirms that the padding was large enough.

That last check is the most important lesson in the example. The NTT computes the product modulo p, and the cyclic transform computes it modulo x^n − 1. If the true product has degree n or more, the high terms wrap around into the low ones; if a true coefficient is p or larger, it wraps around mod p. Both wraps are silent.

Iterative implementation

The iterative version first permutes the array into bit-reversed index order, then runs log2 n stages. Stage s combines blocks of length 2^s using the root of unity of that order. The code below was checked against schoolbook multiplication on random inputs.

MOD, ROOT = 998244353, 3

def ntt(a, invert=False, mod=MOD):
    n = len(a)                       # must be a power of two dividing mod - 1
    j = 0
    for i in range(1, n):            # bit-reversal permutation
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j |= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    length = 2
    while length <= n:
        w = pow(ROOT, (mod - 1) // length, mod)   # primitive length-th root
        if invert:
            w = pow(w, mod - 2, mod)
        half = length // 2
        for start in range(0, n, length):
            t = 1
            for k in range(start, start + half):
                u, v = a[k], a[k + half] * t % mod
                a[k] = (u + v) % mod
                a[k + half] = (u - v) % mod
                t = t * w % mod
        length <<= 1
    if invert:
        n_inv = pow(n, mod - 2, mod)
        for i in range(n):
            a[i] = a[i] * n_inv % mod

def multiply(f, g, mod=MOD):
    need = len(f) + len(g) - 1
    n = 1
    while n < need:
        n <<= 1
    if (mod - 1) % n:
        raise ValueError("transform length does not divide mod - 1")
    fa = f + [0] * (n - len(f))
    fb = g + [0] * (n - len(g))
    ntt(fa, mod=mod); ntt(fb, mod=mod)
    for i in range(n):
        fa[i] = fa[i] * fb[i] % mod
    ntt(fa, invert=True, mod=mod)
    return fa[:need]

Three details matter. The guard on (mod - 1) % n turns a silent wrong answer into an error when the transform is longer than the prime supports. The twiddle t is recomputed per block here for clarity; production code precomputes a table of powers once per length. And (u - v) % mod relies on Python returning a non-negative remainder; in C or Rust, write u - v + mod before reducing, or a negative intermediate leaks through.

Any modulus: three primes and Garner

Often the answer is wanted modulo a prime that is not NTT-friendly, such as 10^9 + 7 (its p − 1 is divisible only by 2^1), or you want the exact integer product of big numbers. The standard fix is to compute the convolution modulo several NTT primes and reconstruct each coefficient with the Chinese remainder theorem.

First bound the true coefficients. If both inputs have entries below m and length at most L, each coefficient of the product is a sum of at most L products, so it is below L · m^2. With m = 10^9 + 7 and L = 2^20 that is about 2^80. The three primes 998244353, 167772161 and 469762049 multiply to about 2^86, so the residues determine each coefficient uniquely. Reconstruct with Garner's method, which never forms a number larger than the moduli product:

P1, P2, P3 = 998244353, 167772161, 469762049

def multiply_any_mod(f, g, m):
    r1 = multiply([x % P1 for x in f], [x % P1 for x in g], P1)
    r2 = multiply([x % P2 for x in f], [x % P2 for x in g], P2)
    r3 = multiply([x % P3 for x in f], [x % P3 for x in g], P3)
    inv_p1 = pow(P1, -1, P2)
    inv_p1p2 = pow(P1 * P2 % P3, -1, P3)
    out = []
    for a, b, c in zip(r1, r2, r3):
        t1 = (b - a) * inv_p1 % P2
        x12 = a + P1 * t1                      # exact value mod P1*P2
        t2 = (c - x12) % P3 * inv_p1p2 % P3
        out.append((x12 + P1 * P2 % m * t2) % m)
    return out

This costs three forward-and-inverse pairs instead of one, so roughly three times the work. Two primes suffice when the coefficient bound is below their product, which is the first thing to check before paying for the third. The extended Euclid and CRT background is in Number Theory Algorithms.

The negacyclic NTT in lattice cryptography

Lattice schemes such as ML-KEM (Kyber) and ML-DSA (Dilithium) work in the ring of polynomials modulo x^256 + 1, not x^256 − 1. Multiplying in that ring is a negacyclic convolution: terms that wrap past degree 255 come back with their sign flipped. The trick is a 2n-th root of unity ψ with ψ^2 = w. Multiply input coefficient i by ψ^i, run an ordinary length-n cyclic NTT, multiply pointwise, inverse, then multiply output coefficient i by ψ^−i. No zero padding is needed, which halves the transform length compared with the cyclic approach. Real implementations merge the ψ twists into the butterfly twiddles.

ML-DSA uses q = 8380417 = 2^23 − 2^13 + 1. Since 512 divides q − 1, a primitive 512th root exists (1753 is the one in the specification; 1753^256 = −1 mod q), and the full negacyclic NTT runs down to single coefficients. ML-KEM uses q = 3329, where q − 1 = 2^8 · 13. A primitive 256th root exists (17, with 17^128 = −1 mod 3329), but no 512th root does, so x^256 + 1 cannot split into linear factors. ML-KEM therefore stops after seven layers, leaving 128 polynomials of degree less than 2, and multiplies those pairs directly. If you implement ML-KEM, that incomplete NTT is specified exactly and the pointwise step is a small degree-1 multiplication, not a scalar product.

In cryptography two extra rules apply. Arithmetic must run in constant time, so no data-dependent branches in reductions. And since the transform touches secret polynomials, use the vetted reference or a reviewed library rather than your own.

Making it fast

The work is about (n/2) log2 n butterflies, each one modular multiplication plus an add and a subtract. Speed comes from making that multiplication cheap and the memory access regular.

• Avoid division. A hardware % is slow. Montgomery multiplication keeps values in a scaled form so reduction becomes multiplies and a shift; Barrett reduction precomputes an approximate reciprocal. When one operand is a fixed twiddle, Shoup's trick precomputes a quotient per twiddle and is cheaper still.
• Lazy reduction. With p below 2^30, values can be kept in the range 0 to 4p, which still fits a 32-bit word, so most conditional subtractions are skipped and normalised at the end.
• Precompute twiddles in the order the loop consumes them, so the inner loop reads a contiguous array.
• Vectorise and batch. AVX2 and NEON implementations process several butterflies per instruction. On GPUs, homomorphic encryption libraries run thousands of independent NTTs at once, because a single small transform cannot fill the device.
• Know the crossover. For short inputs, schoolbook or Karatsuba multiplication wins. Measure the cutoff on your machine instead of assuming one.

Failure modes

• Length beyond the prime. Asking 998244353 for a transform of 2^24 has no valid root; without a guard the code computes garbage. Check that n divides p − 1.
• Too little padding. Sizing n to max(len f, len g) instead of len f + len g − 1 makes the product wrap cyclically. The output looks plausible and is wrong.
• Coefficient overflow mod p. If true coefficients can exceed p, results are silently reduced. Bound them first, then pick one, two or three primes.
• Wrong root. Using a generator that is not a primitive root, or reusing w for the inverse without inverting it, gives a transform that does not invert. A round-trip test catches this instantly.
• Integer overflow in C. The product of two 30-bit residues needs 64 bits; with a 62-bit prime it needs 128 bits or Montgomery form.
• Forgetting n^−1. The inverse without scaling returns every coefficient multiplied by n.

Trade-offs

What to do next

1. Run the code above and test it against schoolbook multiplication on random inputs of random lengths, including lengths 1 and exact powers of two.
2. Add a round-trip test: the inverse NTT of the NTT must return the input unchanged.
3. Write down the coefficient bound for your workload and choose one, two or three primes from it.
4. Replace % with Montgomery or Barrett reduction, precompute twiddle tables, and benchmark against Karatsuba to find your cutoff.
5. For cryptographic use, read the ML-KEM and ML-DSA NTT definitions, run their known-answer tests, and use a reviewed constant-time library.
6. Continue with Advanced Number Theory for primitive roots and group order.