Offensive Security Researcher & Quantum Cryptography Analyst

Polynomials and the FFT: A Comprehensive Guide to Fast Fourier Transform

Abdul Wahab Junaid
Polynomials and the FFT A Comprehensive Guide

I want to explain how the Fast Fourier Transform (FFT) lets me multiply polynomials — and, more generally, transform between coefficient and value representations of a polynomial — far faster than the naive approach. Multiplying two polynomials of degree n the obvious way costs me $O(n^2)$ operations, but with FFT I can do it in $O(n \log n)$. I care about this algorithm because it sits at the intersection of algebra, signal processing, and algorithm design, and it shows up everywhere from audio processing to big-integer multiplication.

History and Background

I credit the modern, widely-used formulation of the FFT to James Cooley and John Tukey, who published their algorithm in 1965. However, I should note that similar ideas were discovered much earlier — Carl Friedrich Gauss described an equivalent method around 1805, long before the electronic computer existed, though his work went largely unnoticed for over a century. The Cooley-Tukey rediscovery in the context of digital computing transformed fields such as signal processing, telecommunications, and numerical analysis, and it is often cited as one of the most important algorithms of the 20th century.

Problem Statement

I frame the core problem as: given two polynomials A(x) and B(x) of degree less than n, each represented by their coefficients, compute their product C(x) = A(x) · B(x), a polynomial of degree less than 2n, as efficiently as possible. More generally, I want an efficient way to convert a polynomial between its coefficient representation and its value representation (its values at a specific set of points), since multiplication is easy in the value representation but hard in the coefficient representation.

Core Concepts

How It Works

I break polynomial multiplication via FFT into three phases:

  1. Evaluation (DFT): I evaluate both input polynomials A and B at the n-th roots of unity, using the FFT algorithm, which does this in $O(n \log n)$ time instead of the naive $O(n^2)$.
  2. Pointwise multiplication: I multiply the evaluated values together pointwise: C(x_k) = A(x_k) · B(x_k) for each root of unity x_k. This step costs $O(n)$.
  3. Interpolation (Inverse DFT): I convert the value representation of C back into coefficient form using the inverse FFT, again in $O(n \log n)$ time.

Before I start, I pad both polynomials with zero coefficients so their combined degree fits within a power of two n large enough to hold the product (n ≥ 2 · degree + 1), since the FFT’s divide-and-conquer structure works most naturally on power-of-two sizes.

Working Principle

I rely on the divide-and-conquer trick at the heart of FFT: I split a polynomial’s coefficients into even-indexed and odd-indexed subsets, forming two smaller polynomials of half the size. I recursively compute the DFT of each half, then combine the results using the special algebraic structure of roots of unity — in particular, the fact that $\omega_n^{k + n/2} = -\omega_n^k$, which lets me reuse each recursive result twice (once with a positive sign, once with a negative sign) instead of recomputing it. This halving-and-combining is what turns an $O(n^2)$ computation into an $O(n \log n)$ one.

Mathematical Foundation

I define the n-th principal root of unity as:

$$ \omega_n = e^{\frac{2\pi i}{n}} $$

The DFT of a coefficient vector $(a_0, \dots, a_{n-1})$ is the vector $(y_0, \dots, y_{n-1})$ where:

$$ y_k = \sum_{j=0}^{n-1} a_j , \omega_n^{jk}, \quad k = 0, 1, \dots, n-1 $$

The Cooley-Tukey recursive decomposition splits A(x) into even and odd terms:

$$ A(x) = A_{\text{even}}(x^2) + x \cdot A_{\text{odd}}(x^2) $$

which lets me write, for $k

$$ A(\omega_n^k) = A_{\text{even}}(\omega_{n/2}^k) + \omega_n^k \cdot A_{\text{odd}}(\omega_{n/2}^k) $$

$$ A(\omega_n^{k + n/2}) = A_{\text{even}}(\omega_{n/2}^k) – \omega_n^k \cdot A_{\text{odd}}(\omega_{n/2}^k) $$

The inverse DFT recovers the coefficients using the conjugate root of unity and a normalization factor:

$$ a_j = \frac{1}{n} \sum_{k=0}^{n-1} y_k , \omega_n^{-jk} $$

Diagrams

flowchart TD
    A["Coefficient vectors A, B (size n)"] --> B["FFT: evaluate A at roots of unity"]
    A --> C["FFT: evaluate B at roots of unity"]
    B --> D["Pointwise multiply: C(x_k) = A(x_k) * B(x_k)"]
    C --> D
    D --> E["Inverse FFT: interpolate back to coefficients"]
    E --> F["Coefficient vector of product polynomial C"]

Pseudocode

FFT(a)                      // a is a coefficient array of length n (power of 2)
    n = length(a)
    if n == 1
        return a
    w_n = e^(2*pi*i / n)
    w = 1
    a_even = FFT(a[0, 2, 4, ...])
    a_odd  = FFT(a[1, 3, 5, ...])
    y = array of size n
    for k = 0 to n/2 - 1
        y[k]         = a_even[k] + w * a_odd[k]
        y[k + n/2]   = a_even[k] - w * a_odd[k]
        w = w * w_n
    return y

MULTIPLY-POLYNOMIALS(A, B)
    n = smallest power of 2 >= 2 * max(len(A), len(B))
    pad A and B with zeros to length n
    yA = FFT(A)
    yB = FFT(B)
    yC[k] = yA[k] * yB[k] for all k
    C = INVERSE-FFT(yC)     // divide every entry by n
    return C

Step-by-Step Example

I multiply A(x) = 1 + 2x (coefficients [1, 2]) by B(x) = 3 + 4x (coefficients [3, 4]).

By hand, I know the answer should be:

$$ (1 + 2x)(3 + 4x) = 3 + 4x + 6x + 8x^2 = 3 + 10x + 8x^2 $$

so the coefficient vector should be [3, 10, 8].

Using FFT conceptually: I pad both to length n = 4 (the next power of 2 that fits a degree-2 result): A = [1, 2, 0, 0], B = [3, 4, 0, 0]. I evaluate both at the 4th roots of unity {1, i, -1, -i}, multiply the evaluations pointwise, and apply the inverse FFT. After interpolation, I recover the coefficient vector [3, 10, 8, 0], which matches my hand computation once I drop the trailing zero.

Time Complexity

Space Complexity

I need $O(n)$ auxiliary space to hold the recursive arrays created during the FFT (for the even/odd splits) and $O(n)$ for the output arrays, so the overall space complexity is $O(n)$. An iterative, in-place FFT implementation (using bit-reversal permutation) can reduce constant factors but keeps the same asymptotic space bound.

Correctness Analysis

I justify correctness in two parts. First, the DFT/inverse-DFT pair is mathematically exact: I can prove using orthogonality of the roots of unity that applying the inverse DFT formula to the DFT values exactly recovers the original coefficients (subject to floating-point rounding in practice). Second, the recursive even/odd decomposition I use in the FFT is an exact algebraic identity — I am not approximating the DFT, I am just computing it in a smarter order — so the recursive FFT computes precisely the same values as the direct $O(n^2)$ DFT formula, just faster.

Advantages

Disadvantages

Applications

Implementation in C

#include stdio.h>
#include math.h>
#include complex.h>

/* I implement a recursive, in-place-style FFT over complex numbers.
   n must be a power of 2. sign = -1 for forward FFT, +1 for inverse FFT. */
void fft(double complex *a, int n, int sign) {
    if (n == 1) {
        return;
    }

    /* I split into even and odd indexed coefficients */
    double complex even[n / 2], odd[n / 2];
    for (int i = 0; i  n / 2; i++) {
        even[i] = a[2 * i];
        odd[i]  = a[2 * i + 1];
    }

    /* I recursively transform each half */
    fft(even, n / 2, sign);
    fft(odd, n / 2, sign);

    /* I combine the two halves using roots of unity */
    for (int k = 0; k  n / 2; k++) {
        double complex w = cexp(sign * 2.0 * I * M_PI * k / n);
        double complex t = w * odd[k];
        a[k]         = even[k] + t;
        a[k + n / 2] = even[k] - t;
    }
}

/* I multiply two real-coefficient polynomials using FFT-based convolution */
void multiplyPolynomials(double *A, int lenA, double *B, int lenB, double *result) {
    int resultSize = lenA + lenB - 1;
    int n = 1;
    while (n  resultSize) n  1;   /* I round up to the next power of 2 */

    double complex fa[n], fb[n];
    for (int i = 0; i  n; i++) {
        fa[i] = (i  lenA) ? A[i] : 0.0;
        fb[i] = (i  lenB) ? B[i] : 0.0;
    }

    fft(fa, n, -1);
    fft(fb, n, -1);

    for (int i = 0; i  n; i++) {
        fa[i] = fa[i] * fb[i];
    }

    fft(fa, n, 1); /* inverse FFT */

    for (int i = 0; i  resultSize; i++) {
        result[i] = creal(fa[i]) / n;  /* I normalize and take the real part */
    }
}

int main(void) {
    double A[] = {1, 2};   /* 1 + 2x */
    double B[] = {3, 4};   /* 3 + 4x */
    int resultSize = 2 + 2 - 1;
    double result[resultSize];

    multiplyPolynomials(A, 2, B, 2, result);

    printf("Product coefficients: ");
    for (int i = 0; i  resultSize; i++) {
        printf("%.0f ", result[i]);
    }
    printf("\n");

    return 0;
}

I wrote fft as a straightforward recursive implementation for clarity rather than maximum performance; a production-grade version would use an iterative, in-place implementation with bit-reversal permutation to avoid the overhead of recursive array splitting.

Sample Input and Output

Input: A(x) = 1 + 2x, B(x) = 3 + 4x

Output:

Product coefficients: 3 10 8

which matches 3 + 10x + 8x^2 exactly, as I computed by hand earlier.

Optimization Techniques

Common Mistakes

Further Reading

Exit mobile version