Fast Fourier Transform (FFT)

Polynomial multiplication combines coefficient arrays through convolution. The direct nested-loop method is quadratic. FFT evaluates polynomials at carefully structured complex roots of unity, multiplies values pointwise, and transforms back in near-linearithmic time. This advanced lesson derives the split and provides an iterative complex implementation for small bounded integer coefficients, with explicit rounding caveats and an introduction to exact modular NTT.

Prerequisites: Complex numbers (real/imaginary pairs), roots of unity, polynomials, divide-and-conquer and modular fields.

Build the idea

1. 1. Define coefficient convolution

   For A(x)=Σaᵢxⁱ and B(x)=Σbⱼxʲ, coefficient cₖ of their product is Σᵢ₊ⱼ₌ₖaᵢbⱼ. Arrays of lengths n,m produce n+m−1 coefficients, including zeros. For [1,2,3] and [4,5], the result is [4,13,22,15]. Index is polynomial degree, so input is lowest degree first.

   $$c_k=\sum _{i+j=k}{a}_i{b}_j$$

2. 2. Evaluate at roots of unity

   Choose a power-of-two N≥n+m−1 and an Nth complex root ω=exp(2πi/N). The DFT records A(ωʲ) for j=0..N−1. Pointwise products equal (AB)(ωʲ). Zero-padding prevents coefficients beyond degree N−1 from wrapping into smaller indices as cyclic convolution. An inverse DFT recovers coefficients using the opposite exponent sign and division by N.

3. 3. Derive the butterfly split

   Separate A(x)=E(x²)+xO(x²) into even- and odd-index coefficients. Squaring Nth roots gives N/2th roots, so recursively evaluate E and O once each. Because ω^(j+N/2)=−ωʲ, paired values are E(ω²ʲ)±ωʲO(ω²ʲ). This is the butterfly operation. Two half transforms plus O(N) combination work gives O(N log N).

   $$A\left(x\right)=E\left(x\mathrm{²}\right)+xO\left(x\mathrm{²}\right)$$

4. 4. Read the iterative implementation

   Bit-reversal permutation puts coefficients in the order needed for bottom-up butterflies. Lengths 2,4,8,… combine already transformed halves with successive root powers. The inverse uses the conjugate-angle direction and divides every entry by N. This implementation uses positive forward angles and negative inverse angles consistently; the opposite convention also works if both ends agree.

5. 5. Bound rounding and know the NTT alternative

   The sample transforms back near [4,13,22,15], then rounds. Input lengths≤4096 and |coefficients|≤100 keep exact answers modest, but complex arithmetic is still numerical. A visible noninteger/imaginary residual triggers NUMERICAL WARNING; this heuristic is not an exactness certificate and does not detect every possible error. For exact modular convolution, NTT uses a prime p with a suitable power-of-two divisor of p−1 and a primitive root; supported N must divide p−1. It avoids floating rounding but returns coefficients modulo p; recovering larger exact integers may require several moduli and CRT.

This is a complete program. Compile with C++17 and enter the input described below.

// Key idea: 5. Bound rounding and know the NTT alternative The sample transforms back near
// [4,13,22,15], then rounds.
#include<iostream>
#include<vector>
#include<complex>
#include<cmath>
#include<algorithm>
using Complex=std::complex<long double>;
void fft(std::vector<Complex>&a,bool inverse){
    int n=(int)a.size();
    for(int i=1,j=0;i<n;++i){
        int bit=n>>1;
        for(;j&bit;bit>>=1){
            j^=bit;
        }
        j^=bit;
        if(i<j){
            std::swap(a[i],a[j]);
        }
    }
    const long double PI=std::acos(-1.0L);
    for(int length=2;length<=n;length<<=1){
        long double angle=2*PI/length*(inverse?-1:1);
        Complex step(std::cos(angle),std::sin(angle));
        for(int start=0;start<n;start+=length){
            Complex root(1,0);
            for(int j=0;j<length/2;++j){
                Complex even=a[start+j],odd=a[start+j+length/2]*root;
                a[start+j]=even+odd;
                a[start+j+length/2]=even-odd;
                root*=step;
            }
        }
    }
    if(inverse){
        for(auto&value:a){
            value/=n;
        }
    }
}
int main(){
    int n,m;
    std::cin>>n>>m;
    int size=1;
    while(size<n+m-1){
        size<<=1;
    }
    std::vector<Complex>a(size),b(size);
    for(int i=0;i<n;++i){
        int x;
        std::cin>>x;
        a[i]=Complex(x,0);
    }
    for(int i=0;i<m;++i){
        int x;
        std::cin>>x;
        b[i]=Complex(x,0);
    }
    fft(a,false);
    fft(b,false);
    for(int i=0;i<size;++i){
        a[i]*=b[i];
    }
    fft(a,true);
    for(int i=0;i<n+m-1;++i){
        if(std::abs(a[i].imag())>1e-4L||std::abs(a[i].real()-std::round(a[i].real()))>1e-4L){
            std::cout<<"NUMERICAL WARNING\n";
            return 0;
        }
    }
    for(int i=0;i<n+m-1;++i){
        std::cout<<std::llround(a[i].real())<<(i+1==n+m-1?'\n':' ');
    }
}

Complexity analysis

Let N be the smallest power of two at least n+m−1. O(N log N) time for two forward FFTs and one inverse, plus O(N) pointwise work; O(N) auxiliary memory. This is O((n+m) log(n+m)) up to power-of-two padding. Numerical precision and coefficient growth impose limits beyond this operation count.

Example input and output

Read n m (1..4096), then n and m integer coefficients, each in −100..100, ordered by increasing degree. Print n+m−1 rounded product coefficients, or NUMERICAL WARNING if the residual heuristic detects trouble. Arithmetic is numerical despite integer input/output.

Input

3 2
1 2 3
4 5

Output

4 13 22 15

(1+2x+3x²)(4+5x)=4+13x+22x²+15x³.

Practical applications

Use for polynomial products, pair-sum frequency convolution and correlation reductions when numeric accuracy is sufficient. Direct convolution is simpler and often faster for short arrays. NTT is preferable for exact modular tasks with a supported prime/length, and multi-modulus reconstruction needs verified coefficient bounds.

Check your understanding

Why zero-pad to at least n+m−1?

Put it into practice

Compare small random convolutions against the direct double loop.