Быстрое Преобразование ФурьеFast Fourier Transform

Многочлен степени n-1 можно задать двумя равноценными способами: коэффициентами (a0, a1, ..., a_{n-1}) — привычная запись — или значениями в n различных точках (x0,y0), ..., (x_{n-1}, y_{n-1}) — потому что n точек однозначно определяют многочлен степени A degree-(n-1) polynomial can be specified two equivalent ways: by its coefficients (a0, a1, ..., a_{n-1}) — the familiar form — or by its values at n distinct points (x0,y0), ..., (x_{n-1}, y_{n-1}) — because n points uniquely determine a polynomial of degree

Трюк БПФ: выбрать не произвольные n точек, а специальные — n-е комплексные корни из единицы (n степень двойки), которые обладают симметрией ω^(k+n/2) = −ω^k. Эта симметрия позволяет вычислить значения многочлена P(x) степени n-1 во всех n корнях из единицы через значения двух многочленов вдвое меньшей степени — P_even (составленного из коэффициентов при чётных степенях x) и P_odd (из коэффициентов при нечётных степенях) — в квадратах этих же корней, которые сами являются корнями из единицы степени n/2. Отсюда рекурсия: P(x) = P_even(x²) + x·P_odd(x²). Вычислив P_even и P_odd рекурсивно в n/2 точках каждый (это n/2-е корни из единицы), можно получить все n значений P за O(n) дополнительной работы благодаря симметрии ω^(k+n/2) = −ω^k: P(ω^k) = P_even(ω^(2k)) + ω^k·P_odd(ω^(2k)), и P(ω^(k+n/2)) = P_even(ω^(2k)) − ω^k·P_odd(ω^(2k)) — та же пара значений P_even/P_odd используется дважды с разным знаком. Рекуррентное соотношение T(n) = 2T(n/2) + O(n) по Master theorem (2.1.2) даёт O(n log n) для evaluation. Обратное преобразование (interpolation) — это, что не очевидно сразу, применение практически того же алгоритма БПФ с корнями ω^(-1) вместо ω и делением на n в конце — тоже O(n log n). Итог: O(n log n) + O(n) + O(n log n) = O(n log n) вместо O(n²) для всей свёртки.The FFT's trick: instead of arbitrary points, choose special ones — the n-th complex roots of unity (n a power of two), which have the symmetry ω^(k+n/2) = −ω^k. This symmetry lets you evaluate a degree-(n-1) polynomial P(x) at all n roots of unity via the values of two half-degree polynomials — P_even (built from the even-indexed coefficients) and P_odd (from the odd-indexed ones) — at the squares of those same roots, which are themselves the (n/2)-th roots of unity. Hence the recursion: P(x) = P_even(x²) + x·P_odd(x²). Having computed P_even and P_odd recursively at n/2 points each (the (n/2)-th roots of unity), all n values of P can be obtained with O(n) extra work thanks to the symmetry ω^(k+n/2) = −ω^k: P(ω^k) = P_even(ω^(2k)) + ω^k·P_odd(ω^(2k)), and P(ω^(k+n/2)) = P_even(ω^(2k)) − ω^k·P_odd(ω^(2k)) — the same pair of P_even/P_odd values is reused twice with opposite signs. The recurrence T(n) = 2T(n/2) + O(n), by the Master theorem (2.1.2), gives O(n log n) for evaluation. The inverse transform (interpolation) — not obvious at first — turns out to be essentially the same FFT algorithm with ω^(-1) instead of ω and a division by n at the end — also O(n log n). The total: O(n log n) + O(n) + O(n log n) = O(n log n) for the entire convolution, instead of O(n²).

// evaluation: рекурсивный БПФ (Кули-Тьюки), n — степень двойки// evaluation: recursive FFT (Cooley-Tukey), n a power of two
std::vector<std::complex<double>> fft(std::vector<std::complex<double>> a, bool invert) {
    int n = a.size();
    if (n == 1) return a;

    std::vector<std::complex<double>> a_even(n/2), a_odd(n/2);
    for (int i = 0; i < n/2; i++) { a_even[i] = a[2*i]; a_odd[i] = a[2*i+1]; }

    auto y_even = fft(a_even, invert);         // P_even в n/2-х корнях// P_even at the n/2-th roots
    auto y_odd  = fft(a_odd, invert);          // P_odd в n/2-х корнях// P_odd at the n/2-th roots

    double ang = 2 * M_PI / n * (invert ? -1 : 1);
    std::complex<double> w(1), wn(cos(ang), sin(ang));
    std::vector<std::complex<double>> y(n);
    for (int k = 0; k < n/2; k++) {
        auto t = w * y_odd[k];
        y[k]       = y_even[k] + t;      // P(ω^k)// P(ω^k)
        y[k + n/2] = y_even[k] - t;      // P(ω^(k+n/2)) = P_even - ω^k·P_odd// P(ω^(k+n/2)) = P_even - ω^k·P_odd
        w *= wn;
    }
    return y;
}
// умножение многочленов: fft(a), fft(b), поточечно перемножить, обратный fft, / n// polynomial multiply: fft(a), fft(b), multiply pointwise, inverse fft, / n

Полное преобразование Фурье требует комплексной арифметики, которую сложно честно показать в простой SVG-анимации, не создавая ложного впечатления точности. Демонстрация ниже намеренно ограничена тем, что можно показать без искажений: разбиение многочлена на чётную/нечётную части рекурсивно вниз до базового случая (структура дерева рекурсии Кули—Тьюки) — именно та часть алгоритма, которая даёт множитель log n. Шаг интерполяции (обратное преобразование) на диаграмме обозначен свёрнуто как единый блок, а не прослеживается арифметически — это сделано намеренно, чтобы не создавать иллюзию полной точности вычислений с комплексными числами.A full Fourier transform requires complex arithmetic, which is hard to show honestly in a simple SVG animation without creating a false impression of precision. The demo below is deliberately limited to what can be shown without distortion: splitting the polynomial into even/odd parts recursively down to the base case (the Cooley–Tukey recursion-tree structure) — exactly the part of the algorithm that produces the log n factor. The interpolation step (the inverse transform) is deliberately shown collapsed as a single summarized block on the diagram, not traced arithmetically — this is intentional, to avoid creating an illusion of full precision in complex-number computation.