Educational implementations

fft.cpp

typedef complex<double> cpx;
typedef vector<cpx> Poly;
typedef vector<cpx> Eval;
Eval FFT(Poly P) {
    int n = P.size();
    if (n == 1)
        return P;
    Poly P_even(n / 2), P_odd(n / 2);
    for (int j = 0; j < n / 2; j++) {
        P_even[j] = P[j * 2];    // Put all the even terms (2*0, 2*1, 2*2,...) in one polynomial
        P_odd[j] = P[j * 2 + 1]; // odd terms in the other (2*0+1,2*1+1,2*2+1,...)
    }
    Eval eval_even = FFT(P_even); // Recursively evaluate the two polynomials
    Eval eval_odd = FFT(P_odd);
    cpx root_of_unity = exp(2i * M_PI / (double)n); // Generate our n-th root of unity
    Eval result(n);
    for (int j = 0; j < n; j++) {
        cpx j_root_of_unity = pow(root_of_unity, j); // Generates \omega^j
        if (j < n / 2)
            result[j] = eval_even[j] + j_root_of_unity * eval_odd[j];
        else // As (\omega^(j + n/2))^2 = (\omega^j)^2, we're reusing them here.
            result[j] = eval_even[j - (n / 2)] + j_root_of_unity * eval_odd[j - (n / 2)];
    }
    return result;
}