📐 Комплексни числа и FFT: Експертно ниво

Тази тема покрива работата с комплексни числа в състезателното програмиране и един от най-важните алгоритми в компютърните науки – Бързото преобразуване на Фурие (Fast Fourier Transform - FFT).

1. Комплексни числа (Complex Numbers)

Комплексните числа са разширение на реалните числа. Те имат вида $z = a + bi$, където: * $a$ е реална част ($ ext{Re}(z)$). * $b$ е имагинерна част ($ ext{Im}(z)$). * $i$ е имагинерната единица, за която $i^2 = -1$.

1.1. Аритметични операции

Ако имаме две числа $z_1 = a_1 + b_1 i$ и $z_2 = a_2 + b_2 i$: * Събиране: $(a_1 + a_2) + (b_1 + b_2)i$ * Изваждане: $(a_1 - a_2) + (b_1 - b_2)i$ * Умножение: $(a_1 a_2 - b_1 b_2) + (a_1 b_2 + a_2 b_1)i$ (защото $i^2 = -1$)

1.2. Геометрична интерпретация

В комплексната равнина (Argand diagram), числото $z = x + iy$ се представя като точка $(x, y)$. * Модул (дължина): $|z| = \sqrt{x^2 + y^2}$. * Аргумент (ъгъл): $\arg(z) = \theta$ е ъгълът с положителната ос $x$. * Полярен вид: $z = r(\cos \theta + i \sin \theta)$, където $r = |z|$. * Формула на Ойлер: $e^{i\theta} = \cos \theta + i \sin \theta$. Следователно $z = r e^{i\theta}$.

Умножението на комплексни числа в полярен вид е много интуитивно: $z_1 \cdot z_2 = (r_1 e^{i\theta_1}) (r_2 e^{i\theta_2}) = (r_1 r_2) e^{i(\theta_1 + \theta_2)}$. * Модулите се умножават. * Ъглите се събират.

1.3. C++ `std::complex`

C++ предоставя вграден клас <complex>.

#include <iostream>
#include <complex>
#include <cmath>

using namespace std;

typedef complex<double> cd;
const double PI = acos(-1);

int main() {
    cd a(3, 4); // 3 + 4i
    cd b(1, -2); // 1 - 2i

    cd sum = a + b;
    cd prod = a * b;

    cout << "Sum: " << sum.real() << " + " << sum.imag() << "i\n";
    cout << "Abs of a: " << abs(a) << "\n"; // 5
    cout << "Arg of a: " << arg(a) << "\n"; // в радиани

    // Създаване чрез полярни координати
    cd c = polar(1.0, PI / 2); // 1 * e^(i*pi/2) = i
    cout << "Polar (0, 1): " << c << "\n";

    return 0;
}

2. Полиноми (Polynomials)

Полином от степен $n-1$ се дефинира като: $A(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_{n-1} x^{n-1}$

Има два основни начина за представяне на полином: 1. Коефициентна форма: Списък с коефициентите $[a_0, a_1, \dots, a_{n-1}]$. * Събиране: $O(n)$. * Умножение: $O(n^2)$ (всеки с всеки). * Пресмятане на стойност (Horner): $O(n)$. 2. Точкова форма (Point-Value form): Множество от $n$ двойки ${\{(x_0, y_0), \dots, (x_{n-1}, y_{n-1})\}}$, където $y_k = A(x_k)$. Всеки полином от степен $n-1$ е уникално определен от $n$ точки. * Събиране: $O(n)$ (събират се стойностите за едни и същи $x$). * Умножение: $O(n)$ (умножават се стойностите $C(x_k) = A(x_k) \cdot B(x_k)$)!

Идеята на FFT: Искаме да умножим два полинома бързо. 1. Преобразуваме ги от коефициентна в точкова форма (DFT). 2. Умножаваме ги за $O(n)$ в точкова форма. 3. Връщаме резултата в коефициентна форма (Inverse DFT).

Ако изберем произволни $x$, преобразуването е бавно ($O(n^2)$). Но ако изберем специални точки – комплексните корени на единицата, можем да го направим за $O(n \log n)$.

3. Бързо преобразуване на Фурие (FFT)

3.1. Корени на единицата

Комплексните $n$-ти корени на единицата са решенията на уравнението $x^n = 1$. Има точно $n$ такива числа, разположени равномерно върху единичната окръжност. $\\omega_n^k = e^{i \frac{2\pi k}{n}} = \cos\left(\frac{2\pi k}{n}\right) + i \sin\left(\frac{2\pi k}{n}\right)$

Основният (примитивен) корен е $\\omega_n = e^{i \frac{2\pi}{n}}$. Останалите са негови степени: $\\omega_n^0, \\omega_n^1, \dots, \\omega_n^{n-1}$.

Важни свойства: 1. $\\omega_n^n = 1$. 2. $\\omega_n^{n/2} = -1$. 3. $\\omega_{2n}^{2k} = \\omega_n^k$.

3.2. Алгоритъмът на Cooley-Tukey

Нека $n$ е четно число. Разделяме полинома $A(x)$ на две части – четни и нечетни индекси: $A(x) = (a_0 + a_2 x^2 + \dots) + x(a_1 + a_3 x^2 + \dots)$ $A(x) = A_{even}(x^2) + x A_{odd}(x^2)$

За да пресметнем $A(x)$ в точките $x = \\omega_n^k$, забелязваме симетрия: За $0 \le k < n/2$: $x = \\omega_n^k \implies x^2 = \\omega_n^{2k} = \\omega_{n/2}^k$ $A(\omega_n^k) = A_{even}(\\omega_{n/2}^k) + \\omega_n^k A_{odd}(\\omega_{n/2}^k)$

За "втората половина" $k' = k + n/2$: $\\omega_n^{k + n/2} = -\\omega_n^k$ $A(\omega_n^{k+n/2}) = A_{even}(\\omega_{n/2}^k) - \\omega_n^k A_{odd}(\\omega_{n/2}^k)$

Така проблемът за размер $n$ се свежда до два подпроблема за размер $n/2$. Сложност: $T(n) = 2T(n/2) + O(n) \implies O(n \log n)$.

3.3. Имплементация (Рекурсивна)

#include <vector>
#include <complex>
#include <cmath>
#include <algorithm>

using namespace std;
typedef complex<double> cd;
const double PI = acos(-1);

void fft(vector<cd> & a, bool invert) {
    int n = a.size();
    if (n == 1) return;

    vector<cd> a0(n / 2), a1(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        a0[i] = a[2 * i];
        a1[i] = a[2 * i + 1];
    }

    fft(a0, invert);
    fft(a1, invert);

    double ang = 2 * PI / n * (invert ? -1 : 1);
    cd w(1), wn(cos(ang), sin(ang));

    for (int i = 0; 2 * i < n; i++) {
        a[i] = a0[i] + w * a1[i];
        a[i + n / 2] = a0[i] - w * a1[i];
        if (invert) {
            a[i] /= 2;
            a[i + n / 2] /= 2;
        }
        w *= wn;
    }
}

3.4. Итеративна имплементация (Оптимизирана)

Рекурсията може да бъде бавна. Итеративният подход използва "Bit-reversal permutation" – индексите се пренареждат според огледалния запис на битовете им.

void fft_iterative(vector<cd> & a, bool invert) {
    int n = a.size();

    // Bit-reversal permutation
    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) swap(a[i], a[j]);
    }

    // Butterfly operations
    for (int len = 2; len <= n; len <<= 1) {
        double ang = 2 * PI / len * (invert ? -1 : 1);
        cd wlen(cos(ang), sin(ang));
        for (int i = 0; i < n; i += len) {
            cd w(1);
            for (int j = 0; j < len / 2; j++) {
                cd u = a[i + j], v = a[i + j + len / 2] * w;
                a[i + j] = u + v;
                a[i + j + len / 2] = u - v;
                w *= wlen;
            }
        }
    }

    if (invert) {
        for (cd & x : a) x /= n;
    }
}

3.5. Умножение на полиноми

За да умножим два полинома $A$ и $B$: 1. Разширяваме векторите с нули до най-близката степен на 2, която е поне $\deg(A) + \deg(B) + 1$. 2. Прилагаме FFT върху $A$ и $B$. 3. Умножаваме поелементно $A[i] = A[i] \cdot B[i]$. 4. Прилагаме Inverse FFT върху резултата. 5. (Опционално) Закръгляме реалните части към най-близкото цяло число.

vector<int> multiply(vector<int> const& a, vector<int> const& b) {
    vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
    int n = 1;
    while (n < a.size() + b.size()) 
        n <<= 1;
    fa.resize(n);
    fb.resize(n);

    fft_iterative(fa, false);
    fft_iterative(fb, false);

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

    fft_iterative(fa, true);

    vector<int> result(n);
    for (int i = 0; i < n; i++)
        result[i] = round(fa[i].real());

    // Премахване на водещи нули ако е нужно
    return result;
}

4. Number Theoretic Transform (NTT)

FFT използва комплексни числа (double), което води до грешки от закръгляне. При задачи, изискващи резултат по модул, използваме NTT. Това е еквивалентът на FFT в крайни полета (Modular Arithmetic). * Вместо комплексни корени на единицата, използваме примитивни корени по модул $P$. * Изисква $P$ да е от специален вид: $P = c \cdot 2^k + 1$. * Пример: $998244353 = 119 \cdot 2^{23} + 1$. Примитивен корен $g=3$.

Предимство: Работи с цели числа (long long), няма загуба на точност.

Приложение на NTT

Всичко е същото като FFT, но: 1. Всички операции са по модул $P$. 2. $\\omega_n$ се заменя с $g^{(P-1)/n} \pmod P$. 3. $\\omega_n^{-1}$ се заменя с модулярното обратно.

const int mod = 998244353;
const int root = 3; // Примитивен корен
// Имплементацията изисква модулно степенуване и модулно обратно

5. Fast Walsh-Hadamard Transform (FWHT)

FWHT е вариант на FFT, който се използва за операции с битови маски. Ако искаме да намерим полином $C$, такъв че: $C[k] = \sum_{i \oplus j = k} A[i] B[j]$ (XOR конволюция) Или за AND/OR операции.

При XOR FWHT матрицата на трансформация е много проста: $\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$ Това позволява умножение за $O(N \log N)$ вместо $O(N^2)$.

6. Задачи за упражнение

SPOJ POLARIVAL: Оценка на стойност на полином.
SPOJ MUL: Умножение на големи числа (разглеждаме цифрите като коефициенти).
Codeforces 993E: Nikita and Order Statistics (броене на подмасиви).
CSES Problem Set: Polynomial Queries.

FFT е мощен инструмент не само за полиноми, но и за всякакви задачи, включващи комбинаторика и суми с фиксиран индекс (конволюции).

```