Fast Polynomial Multiplication, in depth: choosing between schoolbook, Karatsuba, FFT and NTT, measured float precision, and Newton inversion and division

The menu of methods

Every fast method trades multiplications for additions and bookkeeping. They differ in asymptotic cost, in the constant factor that decides real speed at your sizes, and in whether the result is exact.

Two facts drive the choice. The asymptotics hide large constants: an FFT pays for padding, complex arithmetic and memory traffic. And a floating-point FFT returns approximations that round correctly only while the error stays below one half. Production libraries such as GMP therefore keep several of these methods and switch between them at size thresholds tuned per CPU.

How to choose

Karatsuba on coefficient lists

Karatsuba splits each polynomial at the halfway point h, so a = a0 + x^h a1 and b = b0 + x^h b1. The product needs a0 b0, a1 b1 and the middle term a0 b1 + a1 b0. The trick is that the middle term equals (a0 + a1)(b0 + b1) - a0 b0 - a1 b1, so three recursive products replace four. Three half-size products per level gives n^(log2 3). The code below works on coefficient lists over any ring where you can add, subtract and multiply, and falls back to schoolbook below a cutoff.

def schoolbook(a, b):
    out = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        if x:
            for j, y in enumerate(b):
                out[i + j] += x * y
    return out


def _add(a, b):
    n = max(len(a), len(b))
    return [(a[i] if i < len(a) else 0) + (b[i] if i < len(b) else 0) for i in range(n)]


def karatsuba(a, b, cutoff=32):
    if not a or not b:
        return []
    if min(len(a), len(b)) <= cutoff:
        return schoolbook(a, b)
    h = max(len(a), len(b)) // 2
    a0, a1 = a[:h], a[h:]          # a = a0 + x^h * a1
    b0, b1 = b[:h], b[h:]
    z0 = karatsuba(a0, b0, cutoff)
    z2 = karatsuba(a1, b1, cutoff)
    z1 = karatsuba(_add(a0, a1), _add(b0, b1), cutoff)
    out = [0] * (len(a) + len(b) - 1)
    for i, v in enumerate(z0):
        out[i] += v
        out[i + h] -= v
    for i, v in enumerate(z2):
        out[i + 2 * h] += v
        out[i + h] -= v
    for i, v in enumerate(z1):
        out[i + h] += v
    return out

Below the cutoff, slicing and adding lists costs more than the saved multiplication, so time both methods on your coefficient type and pick the crossover.

Evaluation and interpolation with the FFT

The FFT route rests on one observation. A polynomial of degree d is fixed by its values at d + 1 points, and at any point the value of a product is the product of the values. So: evaluate both inputs at enough points, multiply the values pairwise, and interpolate back. Done naively, evaluation and interpolation each cost n^2. The FFT makes both cost n log n by choosing the points to be the complex roots of unity, whose symmetries let each evaluation reuse half of the work of the previous one.

Worked example. Take a(x) = 1 + 2x + 3x^2 and b(x) = 4 + 5x. The product has degree 3, so four points suffice: the fourth roots of unity 1, i, -1 and -i.

Interpolating those four values gives 4 + 13x + 22x^2 + 15x^3, the same as schoolbook (1*4 = 4; 1*5 + 2*4 = 13; 2*5 + 3*4 = 22; 3*5 = 15). Check one row: the product evaluated at 1 is 4 + 13 + 22 + 15 = 54. In code, NumPy's real-input transforms do the work:

import numpy as np


def fft_mul(a, b):
    """Exact only while every true output coefficient stays well below 2**53."""
    n = len(a) + len(b) - 1
    size = 1 << (n - 1).bit_length()
    fa = np.fft.rfft(np.asarray(a, dtype=np.float64), size)
    fb = np.fft.rfft(np.asarray(b, dtype=np.float64), size)
    return np.rint(np.fft.irfft(fa * fb, size)[:n]).astype(np.int64)


def fft_mul_mod(a, b, m, bits=15):
    """Product mod m for coefficients < 2**30: split each into 15-bit halves."""
    mask = (1 << bits) - 1
    a = np.asarray(a, dtype=np.int64)
    b = np.asarray(b, dtype=np.int64)
    a_lo, a_hi = a & mask, a >> bits
    b_lo, b_hi = b & mask, b >> bits
    lo = fft_mul(a_lo, b_lo) % m
    mid = (fft_mul(a_lo, b_hi) % m + fft_mul(a_hi, b_lo) % m) % m
    hi = fft_mul(a_hi, b_hi) % m
    s = (1 << bits) % m
    return (lo + mid * s % m + hi * (s * s % m) % m) % m

Padding to a power of two at least n + m - 1 is essential. A shorter transform computes a cyclic convolution, and the high coefficients wrap around onto the low ones.

How far can floating point go?

Doubles carry 53 bits of mantissa, and FFT rounding error grows with transform size and input magnitude. Rather than quote a rule of thumb, we measured it. Two random polynomials of length n with coefficients uniform in [0, B), multiplied with fft_mul above (NumPy, pocketfft backend), compared against an exact big-integer product:

The error crossed one half once the largest coefficient passed roughly 2^49 to 2^50, well before the 2^53 representation limit. Random inputs like these reach about a quarter of n*B^2; worst-case inputs reach n*B^2 itself, so budget against that. These are single random trials on one machine. Adversarial inputs, such as all coefficients at the maximum, are worse, so leave a margin and assert it. fft_mul_mod multiplies modulo a prime like 10^9 + 7 by splitting each coefficient into 15-bit halves, doing four products well inside the safe zone and recombining modulo m.

Exact products with the NTT

If the answer is wanted modulo a prime p, you can avoid floating point entirely. Work in the integers mod p, where an element of multiplicative order n plays the role of the complex root of unity. That needs n to divide p - 1. The prime 998244353 equals 119 * 2^23 + 1, so it supports power-of-two transforms up to 2^23, with 3 as a primitive root. The butterfly is the same; only the arithmetic changes.

MOD, G = 998244353, 3          # 998244353 = 119 * 2**23 + 1, 3 is a primitive root


def ntt(a, invert=False):
    n = len(a)
    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_len = pow(G, (MOD - 1) // length, MOD)
        if invert:
            w_len = pow(w_len, MOD - 2, MOD)
        for start in range(0, n, length):
            w = 1
            half = length // 2
            for k in range(start, start + half):
                u, v = a[k], a[k + half] * w % MOD
                a[k], a[k + half] = (u + v) % MOD, (u - v) % MOD
                w = w * w_len % MOD
        length <<= 1
    if invert:
        n_inv = pow(n, MOD - 2, MOD)
        for i in range(n):
            a[i] = a[i] * n_inv % MOD
    return a


def ntt_mul(a, b):
    n = len(a) + len(b) - 1
    size = 1 << (n - 1).bit_length()
    fa = ntt(list(a) + [0] * (size - len(a)))
    fb = ntt(list(b) + [0] * (size - len(b)))
    return ntt([x * y % MOD for x, y in zip(fa, fb)], invert=True)[:n]

The NTT is exact but only modulo p. To recover true integer coefficients larger than p, run the product under two or three NTT-friendly primes and combine with the Chinese remainder theorem. The NTT article works through that three-prime construction, and Number Theory Algorithms covers the modular inverses it relies on. On our test machine, two 10,000-term products took 3.8 s with pure-Python Karatsuba, 1.1 s with this NTT and 7 ms with NumPy's compiled FFT; most of that gap is interpreter overhead, not the algorithm.

Beyond multiplication: Newton inversion and division

Fast multiplication is the engine for other polynomial operations. The most useful is the power series inverse: given a with a[0] nonzero, find b with a*b = 1 mod x^n. Newton's iteration doubles the number of correct terms each step: if b is correct mod x^k, then b(2 - a*b) is correct mod x^2k. The cost is a geometric series of multiplications, so inversion costs a constant times one multiplication of size n.

Division follows by reversal. Write A = Q*B + R with deg R below deg B. Reversing coefficient order turns this into rev(A) = rev(Q)*rev(B) mod x^k, where k is the number of quotient terms. So the quotient is one series inverse and one multiplication, and the remainder is A - Q*B. Long division costs n*m; this costs O(n log n). The code below reuses ntt_mul and MOD from the NTT block.

def poly_inverse(a, n, mul=ntt_mul):
    """b with a*b = 1 mod x^n, coefficients mod MOD; needs a[0] != 0."""
    b = [pow(a[0], MOD - 2, MOD)]
    k = 1
    while k < n:
        k *= 2
        ab = mul(a[:k], b)[:k]                    # a*b mod x^k
        corr = [(-v) % MOD for v in ab]
        corr[0] = (corr[0] + 2) % MOD             # 2 - a*b
        b = mul(b, corr)[:k]                      # b <- b*(2 - a*b)
    return b[:n]


def poly_divmod(a, b):
    """Quotient and remainder of a / b over Z_MOD, deg(b) <= deg(a)."""
    n, m = len(a), len(b)
    if m > n:
        return [0], a
    k = n - m + 1                                 # number of quotient terms
    rev_q = ntt_mul(a[::-1][:k], poly_inverse(b[::-1], k))[:k]
    q = rev_q[::-1]
    qb = ntt_mul(q, b)
    r = [(x - y) % MOD for x, y in zip(a[:m - 1], qb[:m - 1])]
    return q, r

Every code block here was tested against schoolbook on random inputs, including a*inverse(a) = 1 mod x^n and quotient*divisor + remainder = dividend. Keep such a test next to any fast multiplication you ship.

Failure modes

• Cyclic wrap-around. Transform size below n + m - 1 silently folds the high coefficients. Assert the size.
• Silent float rounding. Results look plausible and are off by one in a few coefficients. Check max |x - rint(x)| on every product in debug builds and fail if it exceeds about 0.2.
• Wrong root for the length. An NTT of length 2^24 mod 998244353 has no suitable root and returns garbage. Check that n divides p - 1.
• Unbalanced operands. A 10-term by 100,000-term product via a full FFT wastes most of the work; use schoolbook or split the long operand into blocks.
• Untuned cutoffs. A cutoff copied from another language or library can make Karatsuba slower than schoolbook. Benchmark on your target.

Trade-offs

What to do next

1. Write down your coefficient bound B, typical lengths and whether you need exact integers, results mod p, or approximations.
2. Compute n*B^2 and compare it with the measured table; if it is above about 2^48, plan on splitting or an NTT.
3. Implement schoolbook first and keep it as the test oracle.
4. Add Karatsuba or an FFT, then benchmark to find the crossover on your hardware.
5. Add a randomized test against schoolbook covering odd lengths, length 1 and unbalanced operands.
6. If you need division, square roots or many evaluations, build Newton inversion on top of your multiplication rather than writing long division.
7. Read the Divide and Conquer article for the recurrence analysis behind each cost in the table.