Polynomials & the FFT

By the end of this lesson you will be able to explain why polynomials have two "faces" (coefficients and point-values), why multiplying is easy in one face and hard in the other, how the complex roots of unity let us switch faces in Θ(n lg n) time, and you will be able to trace recursiveFft, bitReverseCopy and iterativeFft line by line, run the inverse transform, and multiply polynomials and huge integers with them. Nothing is assumed except the earlier lessons: Growth of functions (Θ notation) and Divide and conquer (recurrences, master theorem). Number display: every complex number on this page is shown rounded to 3 decimals (the program keeps full double precision), and a tiny part below 0.0005 is shown as 0.

Complex numbers in one minute

The FFT lives on the complex numbers, so we meet them first. A complex number is written z = a + bi, where a and b are ordinary (real) numbers and the special symbol i satisfies i·i = −1. We call a the real part and b the imaginary part. Because a number now has two parts, we draw it as a point (a, b) on a flat sheet, the complex plane: a is how far right, b is how far up.

Think of a treasure map: "walk a paces east, then b paces north". A complex number is exactly that instruction. Adding two instructions means adding the east parts and the north parts. Multiplying by a number of length 1 is a turn: it spins the arrow around the origin without changing its length. That turning behaviour is the whole trick of this chapter.

Small example: i has (a, b) = (0, 1), the point straight up. Multiply i by i: i·i = −1, the point (−1, 0), straight left. Each multiplication by i turned the arrow by a quarter turn (90°), and two quarter turns make a half turn. Three more words we need:

Magnitude |z| = √(a² + b²): the arrow's length. Numbers with |z| = 1 lie on the unit circle (the circle of radius 1 around the origin).
Conjugate z̄ = a − bi: the mirror image across the horizontal axis. For numbers on the unit circle, z̄ = z⁻¹ (multiplying z by z̄ gives |z|² = 1). This is why the inverse transform later uses conjugate roots.
Euler's formula e^{iθ} = cos θ + i sin θ: the point on the unit circle at angle θ (in radians; a full turn is 2π). Multiplying e^{iα}·e^{iβ} = e^{i(α+β)}: angles add.

Pitfall. Do not read "imaginary" as "fake". These numbers are just points on a plane with a rotation rule; every value on this page is computed with two ordinary decimal numbers (re, im). Also: the notation e^{2πi/n} is not the number 2.718 raised to a mysterious power, it is simply the unit-circle point at angle 2π/n, that is cos(2π/n) + i sin(2π/n).

Under the hood. In Dart (and JavaScript on this page) a complex number is a pair of 64-bit doubles. Multiplication (a+bi)(c+di) = (ac − bd) + (ad + bc)i costs four real multiplications; the FFT does a few multiplications per butterfly, so the real cost is a small constant times the complex count. Rounding errors of about 10⁻¹⁶ per operation are why the final answers are rounded (see the interview question on precision).

Complex number = arrow on a plane. Multiply by a unit-circle number = rotate. Conjugate = mirror. Those three ideas are all you need for the rest of the lesson.

Coefficient form: Horner's rule, addition, multiplication

A polynomial is a recipe: "take 3 spoons of the constant, 2 spoons of x, 1 spoon of x²". The recipe (the list of spoons) is the coefficient representation. Reading the recipe is easy; cooking with it (evaluating or multiplying) needs work.

A polynomial is A(x) = a₀ + a₁x + … + a_{n−1}x^{n−1}: a sum of powers of the variable x, each with a number in front (its coefficient). It has degree-bound n (every exponent is strictly below n; the true degree, the biggest exponent with a nonzero coefficient, may be smaller). Its coefficient representation is the vector a = (a₀, …, a_{n−1}); here index j holds the coefficient of x^j, so the maths and the code use the same indices (no shift).

Evaluate means "plug in a number x₀ and compute A(x₀)". Computing every power x₀^j separately would cost about n²/2 multiplications. Horner's rule regroups the sum so each step needs one multiplication and one addition:

A(x₀) = a₀ + x₀(a₁ + x₀(a₂ + … + x₀(a_{n−2} + x₀·a_{n−1}) …))

Work from the innermost bracket outwards: start with r = 0, then for j = n−1 down to 0 do r = a_j + x₀·r. After the last step r = A(x₀). That is Θ(n) time. Watch it run on the small example A(x) = 3 + 2x + x² at x₀ = 2 (expected value 3 + 4 + 4 = 11):

Pitfall. Walking the coefficients in the wrong direction (from a₀ upward) computes a different polynomial (the "reversed" one, x^{n−1}·A(1/x)). Horner must start at the highest coefficient. Another trap: a coefficient of 0 must still be visited (it just contributes r = 0 + x₀·r), otherwise the powers shift.

Analysis. The loop body runs exactly n times (one multiplication, one addition each; the very first multiplication is by r = 0), so Θ(n) time, O(1) extra space. Horner is even optimal in the number of multiplications for a general polynomial (a classical theorem, not proved here). Addition of two polynomials of equal degree-bound n is c_j = a_j + b_j, Θ(n). Multiplication is harder: the product C = A·B has coefficients c_j = Σ_{k=0}^{j} a_k·b_{j−k} (the convolution c = a ⊗ b) because x^k·x^{j−k} = x^j, so every a_k meets every b_m exactly once: Θ(n²) products, and degree-bound m + n − 1 for the result.

Now the grade-school multiplication itself: one frame per output coefficient, showing exactly which pairs a_k·b_{j−k} land on x^j.

Coefficient form: evaluate one point Θ(n) (Horner), add Θ(n), multiply Θ(n²). That last Θ(n²) is the bottleneck the whole chapter attacks.

Dart for these three operations (0-indexed, identical to the maths, so there is no index shift):

Input size → what's feasible. n ≤ 10⁶ coefficients → Horner's Θ(n) evaluation is instant; the Θ(n²) grade-school product is fine up to about n = 2000 (4·10⁶ steps) and hopeless at n = 10⁵ (10¹⁰ steps) → use the FFT multiplication section below.

Point-value form: evaluate, add and multiply point-wise

Picture a curve drawn on graph paper. Instead of the recipe, describe it by pinning it down at several places: "at x=0 it is 3, at x=1 it is 6, at x=2 it is 11". A point-value representation of A is a set of n pairs {(x₀,y₀), …, (x_{n−1},y_{n−1})} with all x_k distinct and y_k = A(x_k). Multiplying two curves means multiplying their heights point by point, which is trivially cheap.

Turning coefficients into point-values is evaluation (Horner at n points: Θ(n²)). The reverse, interpolation, comes next. The uniqueness theorem says n distinct points determine exactly one polynomial of degree-bound n (the Vandermonde matrix of the points is invertible, its determinant ∏_{j<k}(x_k − x_j) is nonzero). In this form:

Add: if C = A + B then C(x_k) = A(x_k) + B(x_k). Θ(n), provided A and B are evaluated at the same points.
Multiply: if C = A·B then C(x_k) = A(x_k)·B(x_k). Θ(n). The catch: C has degree-bound m + n − 1, so we need at least m + n − 1 points (for two length-n inputs, about 2n).

The next players evaluate A and B at chosen points, add and multiply them point-wise, then interpolate back (using Lagrange, explained in the interpolation section) and compare with the coefficient-form answer. The second example deliberately uses too few points, so you see the failure the theorem predicts.

Pitfall. (1) Using only max(m, n) points for a product: the interpolation still succeeds, but it returns a different polynomial of degree-bound N (the one with that many points), not the true product. (2) Evaluating A and B at different point sets: the pointwise product then multiplies unrelated heights. (3) Repeated x-values: two pairs (x, y₁), (x, y₂) with y₁ ≠ y₂ cannot come from any polynomial, and equal points make the Vandermonde determinant 0.

Why not simply evaluate at 2n points and multiply? Because evaluation costs Θ(n²) with Horner at arbitrary points and interpolation costs Θ(n²) with Lagrange: the conversions eat the savings. The Θ(n) pointwise multiply is worthless unless both conversions become cheap. That is the job of the FFT: choose the points as the complex roots of unity so that both conversions cost Θ(n lg n).

Point-value form: add Θ(n), multiply Θ(n), but needs ≥ m + n − 1 points for the product and Θ(n²) to convert with arbitrary points. Pipeline to remember: coefficients → evaluate → pointwise multiply → interpolate → coefficients.

Input size → what's feasible. evaluating at n arbitrary points is Θ(n²): fine for n ≤ 2000, too slow at n = 10⁵; point-wise add/multiply are Θ(n) once the values exist, so the conversion is the bottleneck.

Interpolation and Lagrange's formula

"Connect the dots": given n dots (with no two directly above each other), there is exactly one smooth polynomial curve of degree-bound n that passes through all of them. Lagrange's trick is to build it from n switches: switch k is a polynomial L_k that equals 1 at dot k and 0 at every other dot. Then y₀·L₀ + y₁·L₁ + … passes through every dot, because at dot k all switches except k are off and switch k is on with weight y_k.

A(x) = Σ_k y_k · L_k(x), L_k(x) = ∏_{j≠k} (x − x_j) / (x_k − x_j)

The numerator N_k(x) = ∏_{j≠k}(x − x_j) is a polynomial that is zero at every other point (one factor vanishes). Dividing by the number d_k = N_k(x_k) = ∏_{j≠k}(x_k − x_j) rescales it to equal exactly 1 at x_k. Computing all coefficients this way is Θ(n²) if you are smart: multiply out P(x) = ∏_j (x − x_j) once (Θ(n²)), then get every N_k(x) = P(x)/(x − x_k) by synthetic division (dividing a polynomial by a linear factor in Θ(n)), evaluate d_k = N_k(x_k) with Horner (Θ(n)), and add (y_k/d_k)·N_k into a running total (Θ(n)). Building each N_k from scratch would instead cost Θ(n²) per point, Θ(n³) overall. The player follows the smart version, one basis polynomial per step.

Pitfall. Two points with the same x make d_k = 0: division by zero (the page rejects that input). And Lagrange is numerically delicate: with points that are close together, or many points, tiny rounding errors get magnified . Interpolation at the roots of unity is the well-behaved special case.

Uniqueness theorem in one paragraph. If two polynomials of degree-bound n agreed at n distinct points, their difference would be a polynomial of degree-bound n with n roots; a nonzero polynomial of degree d has at most d roots, so the difference is the zero polynomial. That is the same statement as "the Vandermonde matrix is invertible". Solving the linear system by LU (the matrix-algorithms lesson) would cost O(n³); Lagrange brings it to Θ(n²); the FFT will bring the special case to Θ(n lg n).

Interpolation: n distinct points ⇒ a unique polynomial. Lagrange builds it in Θ(n²). Evaluation and interpolation are inverse conversions between the two faces, each Θ(n²) for arbitrary points.

Dart for Lagrange's basis polynomial L_k and for the uniqueness theorem (Vandermonde determinant, and "n distinct points force equality"). Indices k, j = 0..n−1 are the Dart indices: no shift.

Input size → what's feasible. Lagrange is Θ(n²): fine for n ≤ 1000 well-spread points; at n = 10⁵ it is 10¹⁰ steps (too slow) and numerically unstable → interpolate at the roots of unity with the inverse DFT instead.

Complex roots of unity

Picture a clock face with n evenly spaced marks, where "hour 0" is at 3 o'clock and time runs counter-clockwise. Multiplying by the mark ω = e^{2πi/n} means "advance one mark". After n advances you are back at 1. The n marks are the complex n-th roots of unity: the n numbers z with zⁿ = 1. They are ω_nᵏ = e^{2πik/n} for k = 0..n−1, all on the unit circle. ω_n = e^{2πi/n} is the principal root.

Small example: n = 4. Then ω₄ = e^{2πi/4} = cos 90° + i sin 90° = i, and the four roots are 1, i, −1, −i (the compass points). Check: i⁴ = 1. Here is n = 8 next to what happens when you square its roots:

Three facts make the FFT work. Each has a player below that walks through it on a circle, with a custom-input version:

Cancellation lemma: ω_{dn}^{dk} = ω_n^k (a d-times finer clock, taking d-times bigger steps, lands on the same mark).
Halving lemma: for even n > 0, squaring the n roots gives the n/2 roots of order n/2, each exactly twice: (ω_n^{k+n/2})² = ω_n^{2k+n} = ω_n^{2k} = ω_{n/2}^k.
Summation lemma: for k not divisible by n, Σ_{j=0}^{n−1} (ω_n^k)^j = ((ω_n^k)ⁿ − 1)/(ω_n^k − 1) = 0 (the marks balance around the circle). Also ω_n^{n/2} = −1.

Cancellation lemma

Read the players: the left circle has the n coarse marks, the right circle the dn fine marks. Mark dk on the fine circle sits at exactly the same angle as mark k on the coarse one.

Halving lemma

The left circle shows ω_n^k (blue) and its opposite partner ω_n^{k+n/2} = −ω_n^k (orange). Squaring doubles an angle, so opposite marks (180° apart) land on the same mark of the right circle (n/2 marks, green).

Summation lemma

Draw the powers r⁰, r¹, r², … of r = ω_n^k as arrows placed head to tail. When the arrows form a closed regular polygon, they add up to 0.

Pitfall. ω_n^k depends on k only through k mod n (ω_n^n = 1), and the summation lemma needs "k not a multiple of n": for k = 0 (or k = n) the arrows all point the same way and the sum is n, not 0. Also do not confuse ω_n (a single number, the principal root) with "the n-th roots of unity" (all n numbers ω_n⁰..ω_n^{n−1}). The FFT needs n to be a power of 2 so the halving lemma can be applied repeatedly.

Proofs in one line each. Cancellation: ω_{dn}^{dk} = (e^{2πi/(dn)})^{dk} = e^{2πik/n} = ω_n^k. Halving: (ω_n^{k+n/2})² = ω_n^{2k+n} = ω_n^{2k}·ω_n^n = ω_n^{2k} = ω_{n/2}^k, so k and k + n/2 give the same square; k = 0..n/2−1 gives n/2 distinct values. Summation: with r = ω_n^k ≠ 1 the geometric series is (rⁿ − 1)/(r − 1) and rⁿ = (ω_nⁿ)^k = 1. The summation lemma is what makes the DFT matrix invertible (inverse-matrix theorem).

Squaring collapses n roots to n/2 (halving) and ω_n^{n/2} = −1 pairs every mark with its opposite. Those two facts are the entire reason the FFT can split a size-n problem into two size-n/2 problems.

Dart for a root of unity (angle from k mod n, so big exponents stay accurate):

Dart that checks the three lemmas numerically, and the facts about n-th roots of unity and degree-bounds used above (a degree-bound-n polynomial has at most n roots; a product needs la + lb − 1 coefficients, padded to a power of 2). No index shift: ω_n^k uses k = 0..n−1 in the maths and the code.

Input size → what's feasible. one root of unity costs O(1) (a cos and a sin of the angle of k mod n); all n roots cost Θ(n), instant for n ≤ 10⁶.

The DFT: evaluating at the roots of unity

If you choose to "photograph" the polynomial at the n marks of the clock face, the list of n photo values is its spectrum. The DFT is just that list: the polynomial's values at the n roots of unity.

The DFT (discrete Fourier transform) of a = (a₀,…,a_{n−1}) is y_k = A(ω_nᵏ) = Σ_j a_j ω_n^{kj}, for k = 0..n−1. As a matrix-vector product y = V_n·a, where V_n has entry ω_n^{kj} in row k, column j (a Vandermonde matrix at the roots of unity). Small example, n = 4 and a = (1, 2, 3, 4) (the polynomial 1 + 2x + 3x² + 4x³), computed live by this page:

Straight from the definition this is n sums of n terms: Θ(n²). The FFT (fast Fourier transform) computes exactly the same y in Θ(n lg n) when n is a power of 2 (we assume that throughout).

Pitfall. "FFT" and "DFT" are not two different transforms: the DFT is what is computed, the FFT is how (a fast algorithm). Also the sign convention: this page uses e^{+2πi/n} for the forward transform; many libraries use e^{−2πi/n}. The answers differ by a conjugation/reversal, so match conventions when comparing with other software.

Why the roots? Any n distinct points would give a valid point-value form. The roots are special because their squares collapse (halving lemma) and their powers cancel (summation lemma); those two properties yield a divide-and-conquer evaluation and a beautiful inverse.

DFT = the n values A(ω_n⁰), …, A(ω_n^{n−1}). Definition: Θ(n²). FFT: Θ(n lg n).

Dart for the DFT definition and the inverse DFT matrix. The matrix V_n has entry (k, j) = ω_n^kj with 0-based k, j: the same indices as Dart.

Input size → what's feasible. the definition costs Θ(n²): fine for n ≤ 2000, but n = 2¹⁸ would be 6.9·10¹⁰ steps (too slow) → use the FFT, Θ(n lg n).

recursiveFft(a)

Sorting a deck by splitting it into two halves works because each half is a smaller copy of the same job. The FFT does the same for evaluation. Split the coefficients into even-indexed and odd-indexed ones: A(x) = A⁰(x²) + x·A¹(x²). Evaluating A at n roots of unity needs A⁰ and A¹ only at the squares of those roots — and by the halving lemma, there are only n/2 distinct squares, i.e. the n/2-th roots of unity. So one size-n job becomes two size-n/2 jobs of the same shape, plus n cheap "butterfly" combinations.

Small example first, by hand, n = 2 and a = (a₀, a₁): A(x) = a₀ + a₁x, roots 1 and −1, so y = (a₀ + a₁, a₀ − a₁). Here the even half is (a₀) and the odd half is (a₁), each already its own DFT, and the combine loop computes exactly y₀ = a₀ + 1·a₁ and y₁ = a₀ − 1·a₁. Bigger sizes repeat this pattern on top of each other.

Below, y0 is the DFT of the even half a0 and y1 is the DFT of the odd half a1. The stack shows every active call (a call is one running copy of the procedure, waiting for its children); the circle on the right shows the roots of unity of the running call, with ω_n^k (blue) and its partner −ω_n^k = ω_n^k+n/2 (orange).

The recursion tree and its cost (analysis visual). Each box holds the indices of a that one call receives; every call splits its indices by parity. Lines 10–14 of a call of size s do s/2 butterflies, so every level does n/2 butterflies in total, and there are lg n levels of non-trivial calls:

Pitfall. (1) Forgetting that ω must restart at 1 inside every call: line 5 runs once per call, not once globally. (2) Forgetting that the second output uses the same product ω·y1[k] with a minus sign (u − t), which is valid only because ω_n^{k+n/2} = −ω_n^k. (3) Calling it with n not a power of 2: the parity split then leaves halves of unequal size and the halving lemma does not apply.

Running time. Lines 6–7 and the loop of lines 10–14 do Θ(n) work; two recursive calls have size n/2. So T(n) = 2T(n/2) + Θ(n) = Θ(n lg n) (case 2 of the master theorem, divide-and-conquer lesson). The recursion tree above has lg n + 1 levels, each doing Θ(n) total work. Correctness: lines 11–13 use w = ω_n^k: y_k = A⁰(ω_n^{2k}) + ω_n^k A¹(ω_n^{2k}) = A(ω_n^k) by A(x) = A⁰(x²) + xA¹(x²), and y_{k+n/2} = A⁰(ω_n^{2k}) − ω_n^k A¹(ω_n^{2k}) = A(ω_n^{k+n/2}) because ω_n^{k+n/2} = −ω_n^k and (−ω_n^k)² = ω_n^{2k}. Line 14 keeps w = ω_n^k as k grows (running product instead of recomputing). Invariant: at the start of iteration k, ω = ω_n^k, and y[0], y[1] hold the correct half-size DFTs.

Even/odd split + halving lemma ⇒ T(n) = 2T(n/2) + Θ(n) = Θ(n lg n). One product ω·y1[k] serves two outputs.

Dart implementation (the pseudocode's sign argument lets the same code also compute the inverse's transform later; sign = 1 is this page's convention). Index shift: none, the pseudocode and Dart both count from 0:

The Θ(n lg n) claim in numbers: count the butterflies of iterativeFft for n = 2, 4, …, 2¹⁶, compare with the recurrence T(n) = 2T(n/2) + n/2 and with n².

Input size → what's feasible. n a power of 2 with n ≤ 2¹⁸ → about 2.4·10⁶ butterflies, instant; the recursion is only lg n = 18 levels deep but allocates new lists at every level, so for n ≥ 2²⁰ prefer iterativeFft (the iterative section below).

The inverse DFT (interpolation at the roots of unity)

If the DFT is "recipe → measurements", the inverse is "measurements → recipe". Because the measurement points are the balanced marks on the circle, the matrix that does the forward job (the Vandermonde matrix V_n with entry ω_n^{kj}) has a beautiful inverse: the same matrix with ω replaced by its conjugate ω⁻¹ (mirror image across the horizontal axis) and everything divided by n.

Inverse-matrix theorem: (V_n⁻¹)_{jk} = ω_n^{−kj}/n. Proof idea: entry (j,j′) of V_n⁻¹V_n is Σ_k ω_n^{k(j′−j)}/n, which equals 1 when j′ = j and 0 otherwise by the summation lemma (since |j′−j| < n, j′−j is not divisible by n unless zero). So a_j = (1/n) Σ_k y_k ω_n^{−kj}. Compare with y_k = Σ_j a_j ω_n^{kj}: run the same FFT with the roles of a and y swapped, ω_n replaced by ω_n⁻¹ = e^{−2πi/n} (the conjugate, since |ω| = 1), and divide by n. Hence DFT⁻¹ also costs Θ(n lg n). In the code below we reuse the fast iterative transform (next sections) with sign = −1.

Small example: n = 2. y = (a₀ + a₁, a₀ − a₁). Applying the same butterfly with ω = 1 gives (y₀ + y₁, y₀ − y₁) = (2a₀, 2a₁); dividing by n = 2 recovers (a₀, a₁). The players compute each y′_j = Σ_k y_k ω_n^{−kj} term by term (that is what the sign = −1 butterflies produce; the page checks the result against iterativeFft with the conjugate roots), then divide by n:

Pitfall. The two classic bugs: (1) forgetting the division by n — the output is exactly n times too large (the last frame of each player shows this); (2) using ω instead of ω⁻¹ — then the output is the index-reversed sequence (a₀, a_{n−1}, a_{n−2}, …, a₁) times n instead of a. Applying the forward FFT twice gives n·(a₀, a_{n−1}, …, a₁).

Why it works. V_n·V̄_n = n·I because the (j, j′) entry is Σ_k ω_n^{k(j−j′)} = n if j = j′ and 0 otherwise (summation lemma). So V̄_n/n is the inverse. Another way to see the same fact: the conjugation identity DFT⁻¹(y) = conj(DFT(conj(y)))/n (interview question) turns any forward FFT into an inverse one. Running time Θ(n lg n), same as the forward transform.

Inverse DFT = forward FFT with ω → ω⁻¹ (conjugate roots), then divide by n. Same Θ(n lg n).

Input size → what's feasible. n a power of 2 with n ≤ 2¹⁸ → one FFT plus n divisions ≈ 2.4·10⁶ butterflies; pad shorter inputs with zeros up to the next power of 2.

Polynomial multiplication in Θ(n lg n): evaluate → pointwise multiply → interpolate

To multiply two big numbers by hand you can switch to logarithms, add, and switch back. The FFT pipeline is the same idea: translate both polynomials to the "point-value language" (where multiplication is one product per point), multiply there, and translate the answer back.

Convolution theorem: a ⊗ b = DFT⁻¹_{2n}(DFT_{2n}(a) · DFT_{2n}(b)), where the vectors are padded with zeros to length 2n (more exactly: to a power of 2 that is at least the product's degree-bound) and "·" is the pointwise product. The padding is not optional: with too few points the product's high coefficients wrap around and add onto the low ones (a cyclic convolution). Steps: (1) pad, (2) two forward FFTs, (3) n pointwise multiplications, (4) one inverse FFT. Three transforms, so Θ(n lg n). Worked example computed live by this page — A(x) = 1 + 2x, B(x) = 3 + x, padded to n = 4:

Now step through the three phases yourself: evaluate both polynomials at the n roots of unity, multiply pointwise, then interpolate with the inverse DFT. The page's own FFT computes every value, and each caption cross-checks the result against grade-school (naive) multiplication.

Analysis visuals. (a) What goes wrong without padding: the same product (1 + x)² computed with n = 2 and n = 4. (b) The cost comparison: grade-school multiplies n² pairs; the FFT route does three transforms of (n/2) lg n butterflies plus n pointwise products and n divisions.

Pitfall. (1) Padding to the input length instead of the output length (m + n − 1): coefficients wrap around (the table above). (2) Comparing floating-point results without rounding: c_j comes back as 6.99999999 or 7.0000001; round to the nearest integer if the true answer is integral. (3) Very large coefficients: rounding errors grow with n and with the size of the numbers, so beyond roughly 10¹⁵ in the largest intermediate sum the rounding is no longer trustworthy (see the precision question).

Big integers. An integer is a polynomial evaluated at 10: 4728 = 8 + 2·10 + 7·10² + 4·10³. So multiplying two integers means multiplying the two digit polynomials and then carrying, because a coefficient may now exceed 9. That gives Θ(n lg n) digit operations instead of Θ(n²) for the grade-school method (and beats Karatsuba's Θ(n^{1.585}) for huge n; see divide-and-conquer lesson).

coefficients (padded to n = power of 2 ≥ m + n − 1) → FFT Θ(n lg n) → point-values → pointwise multiply Θ(n) → inverse FFT Θ(n lg n) → product coefficients. Total Θ(n lg n) instead of Θ(n²).

Dart, including big-integer multiplication: write each integer's decimal digits as the coefficients of a polynomial (least-significant digit first). Then x = P(10), y = Q(10), and x·y = (PQ)(10): multiply the polynomials, round each coefficient to an integer, then carry to get decimal digits again.

The convolution theorem as a numeric check (Dart). The transform length n decides whether the product wraps around.

Input size → what's feasible. polynomials with ≤ 10⁵ coefficients each → padded n = 2¹⁸ and ≈ 7·10⁶ butterflies (instant) against 10¹⁰ steps for the grade-school method; below ~50–100 coefficients the plain Θ(n²) loop is just as fast.

bitReverseCopy(a, A)

Imagine a library shelf that you keep sorting into "odd-numbered / even-numbered" halves again and again. After lg n rounds, each volume sits at the position you get by writing its old index in binary and reading the bits backwards. Instead of splitting recursively, the iterative FFT jumps straight to that final shelf order and only does the merging.

Small example first: n = 4. Indices 0,1,2,3 are 00, 01, 10, 11 in two bits; reversed they read 00, 10, 01, 11 = 0, 2, 1, 3. So a₁ goes to slot 2 and a₂ goes to slot 1; a₀ and a₃ stay. (A bit is one binary digit; see D03.)

Look at the recursion of recursiveFft: the first split sends indices ending in bit 0 left and ending in bit 1 right; the next split looks at the second-lowest bit, and so on. So the leaf that holds a_k is at position rev(k), where rev(k) is the lg n-bit binary expansion of k reversed. bitReverseCopy places a_k into A[rev(k)] (lines 2–3 of the panel: for k from 0 to n − 1, A[rev(k)] ← a[k]). The players show a as symbols a₀…, with binary indices as pointers.

Analysis visual. The table lists k, its lg n-bit binary, the reversed bits and the resulting slot for n = 8; the leaf order of the recursion tree above (0, 4, 2, 6, 1, 5, 3, 7) is exactly "which a_k lands in slot 0, 1, 2, …".

Pitfall. rev must use exactly lg n bits, including leading zeros: 1 reversed in 3 bits is 100 = 4, but reversed in 4 bits it is 1000 = 8. Reversing "the bits of the number" without fixing the width gives the wrong slot. Also, copying a into A in place without care overwrites values that have not moved yet; the simplest version copies into a separate array A (or swaps k with rev(k) only when k < rev(k)).

Running time. The loop runs n times; computing rev(k) by shifting bit by bit costs O(lg n), so O(n lg n) total (a lookup table, or a reverse-binary counter as in the amortized analysis lesson, makes it Θ(n)). rev is an involution (rev(rev(k)) = k), so the copy is a permutation and can even be done in place by swapping k with rev(k) when k < rev(k).

bitReverseCopy puts each a_k at slot rev(k) so that iterativeFft can merge neighbours bottom-up; O(n lg n) (Θ(n) with a table).

Input size → what's feasible. n ≤ 2²⁰ → about 2·10⁷ bit-by-bit steps, instant; a reverse-binary counter or table makes it Θ(n).

iterativeFft(a)

A tournament bracket: in round 1 neighbours pair up, in round 2 winners-of-pairs pair up, and so on. iterativeFft does this bottom-up on the bit-reversed array: stage s merges groups of size m/2 = 2^{s−1} into groups of size m = 2^s, and each merge is n/m groups × m/2 butterflies. A butterfly takes two cells u and ω·t (a "twiddle factor" ω, a root of unity, multiplied in) and writes u + ωt and u − ωt back into the same two cells — no extra array needed. Drawn as a wiring diagram the two crossing wires look like a butterfly.

Small example: n = 4, a = (a₀,a₁,a₂,a₃). After bitReverseCopy, A = (a₀,a₂,a₁,a₃). Stage 1 (m = 2, ω = 1) combines neighbours: (a₀+a₂, a₀−a₂, a₁+a₃, a₁−a₃). Stage 2 (m = 4) combines those two blocks with twiddles ω₄⁰ = 1 and ω₄¹ = i. That is the whole algorithm; for larger n there are just more stages.

The butterfly network (analysis visual) for n = 8: wires run left to right, inputs on the left already in bit-reversed order (label a_k means "a_k sits here"), lg n = 3 stages, each with n/2 = 4 butterflies. The numbers on the wires are the twiddle exponents (ω_m^j).

In the players below the current group (lines 6–7) is orange, the two cells overwritten by the latest butterfly (lines 8–13) are red, and a live copy of this network under the array turns green for every butterfly already executed and red for the one just executed. The circle shows ω_m and the twiddle factor ω = ω_m^j in use.

Pitfall. (1) Skipping bitReverseCopy: the butterflies then combine the wrong neighbours and produce garbage that still looks plausible. (2) Restarting ω at 1 for every k (line 7) is required; carrying it over between groups gives wrong twiddles. (3) Reading A[k+j+m/2] after overwriting A[k+j]: the code saves u first (line 10) so both outputs use the old values.

Running time. The innermost body (lines 8–13) runs Σ_{s=1}^{lg n} (n/2^s)·2^{s−1} = Σ n/2 = (n/2) lg n times, so iterativeFft is Θ(n lg n) plus bitReverseCopy. For n = 8 that is 12 butterflies, for n = 1024 it is 5120. Because each stage's n/2 butterflies touch disjoint pairs, they can all run in parallel: a circuit of depth Θ(lg n). Invariant: at the start of stage s, every aligned block of size 2^{s−1} in A holds the DFT of the subsequence of a that the recursion would have assigned to that block.

iterativeFft = bitReverseCopy + lg n stages of n/2 in-place butterflies = Θ(n lg n) time, Θ(n) space, parallel depth Θ(lg n).

Input size → what's feasible. n a power of 2 with n ≤ 2¹⁸ runs in about 0.1 s; in double precision integer products stay exact while coefficients are ≤ ~100 and lengths ≤ 10⁵ (see the multiplication section), beyond that use an NTT (question 23).

Quiz

Interview questions

Cheat sheet

All procedures here are deterministic and data-independent, so best = average = worst case.

Procedure	Time (best = avg = worst)	Space	In place?	Use when
Horner evaluation	Θ(n)	O(1)	yes	One point, coefficient form
Coefficient add	Θ(n)	O(n)	can be	Same degree-bound; either form works
Naive polynomial multiply (convolution)	Θ(n²)	O(n)	—	Small n (below roughly 50–100)
Point-value add / multiply	Θ(n)	O(n)	can be	After conversion; needs ≥ m + n − 1 shared points for a product
Evaluation at n arbitrary points (Horner each)	Θ(n²)	O(n)	—	Arbitrary points
Lagrange interpolation (all coefficients)	Θ(n²)	O(n)	—	Arbitrary distinct points
Naive DFT (definition)	Θ(n²)	O(n)	—	Reference / any n
recursiveFft	Θ(n lg n)	Θ(n lg n) allocated in total, Θ(n) live per level	no	Clarity; n a power of 2
bitReverseCopy	Θ(n lg n) (Θ(n) with table)	Θ(n)	can swap in place	Prepares iterativeFft
iterativeFft	Θ(n lg n) — (n/2) lg n butterflies	Θ(n)	yes (after the copy)	Production, parallel circuit depth Θ(lg n)
Inverse DFT	Θ(n lg n)	Θ(n)	yes	ω → ω⁻¹, divide by n
Polynomial / big-integer multiply via FFT	Θ(n lg n)	Θ(n)	no	Large n; round results, watch precision (use NTT for exactness)