Перейти до основного вмісту

Швидке перетворення Фур'є

У цій статті ми розглянемо алгоритм, який дозволяє перемножити два многочлени довжини nn за час O(nlog⁡n)O(n \log n), що краще за тривіальне множення, яке потребує O(n2)O(n^2) часу. Очевидно, що множення двох довгих чисел також можна звести до множення многочленів, тому й два довгих числа можна перемножити за час O(nlog⁡n)O(n \log n) (де nn — кількість цифр у числах).

Відкриття швидкого перетворення Фур'є (ШПФ, англ. FFT) приписують Кулі та Тьюкі, які опублікували алгоритм у 1965 році. Але насправді FFT відкривали неодноразово й раніше, проте до появи сучасних комп'ютерів його важливість не усвідомлювали. Деякі дослідники приписують відкриття FFT Рунге та Кеніґу в 1924 році. Та насправді Гаусс розробив такий метод іще в 1805 році, але ніколи його не публікував.

Зауважте, що наведений тут алгоритм FFT працює за час O(nlog⁡n)O(n \log n), але він не годиться для множення довільно великих многочленів із довільно великими коефіцієнтами або для множення довільно великих цілих чисел. Він легко справляється з многочленами розміру 10510^5 з невеликими коефіцієнтами або з множенням двох чисел розміру 10610^6, чого зазвичай достатньо для розв'язання задач зі змагального програмування. За межами множення чисел із 10610^6 бітами діапазону та точності чисел з рухомою комою, що використовуються під час обчислень, не вистачить, щоб дати точні кінцеві результати, хоча існують складніші варіації, здатні виконувати множення довільно великих многочленів чи цілих чисел. Наприклад, у 1971 році Шенгаґе та Штрассен розробили варіацію для множення довільно великих чисел, яка рекурсивно застосовує FFT у кільцевих структурах і працює за O(nlog⁡nlog⁡log⁡n)O(n \log n \log \log n). А нещодавно (у 2019 році) Гарві та ван дер Говен опублікували алгоритм, який працює за справжній O(nlog⁡n)O(n \log n).

Коли підходить цей алгоритм?
  • Чи задача зводиться до множення двох многочленів (або згортки двох послідовностей), де наївне O(n2)O(n^2) завелике?
  • Чи розмір даних у межах можливостей FFT (многочлени до ∼105\sim 10^5 з невеликими коефіцієнтами, або числа до ∼106\sim 10^6 цифр), де точності чисел з рухомою комою вистачає?
  • Якщо потрібен точний результат за модулем — (чи варто застосувати теоретико-числове перетворення (NTT), описане нижче, або множення з довільним модулем)?

Дискретне перетворення Фур'є​

Нехай маємо многочлен степеня n−1n - 1:

A(x)=a0x0+a1x1+⋯+an−1xn−1A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}

Без втрати загальності припустимо, що nn — кількість коефіцієнтів — є степенем 22. Якщо nn не є степенем 22, то ми просто додаємо відсутні члени aixia_i x^i і встановлюємо коефіцієнти aia_i рівними 00.

Теорія комплексних чисел каже нам, що рівняння xn=1x^n = 1 має nn комплексних розв'язків (їх називають nn-ми коренями з одиниці), і ці розв'язки мають вигляд wn,k=e2kπinw_{n, k} = e^{\frac{2 k \pi i}{n}} при k=0…n−1k = 0 \dots n-1. Крім того, ці комплексні числа мають деякі дуже цікаві властивості: наприклад, головний nn-й корінь wn=wn,1=e2πinw_n = w_{n, 1} = e^{\frac{2 \pi i}{n}} можна використати, щоб описати всі інші nn-ті корені: wn,k=(wn)kw_{n, k} = (w_n)^k.

Дискретне перетворення Фур'є (DFT) многочлена A(x)A(x) (або, що те саме, вектора коефіцієнтів (a0,a1,…,an−1)(a_0, a_1, \dots, a_{n-1})) означується як значення многочлена в точках x=wn,kx = w_{n, k}, тобто це вектор:

DFT(a0,a1,…,an−1)=(y0,y1,…,yn−1)=(A(wn,0),A(wn,1),…,A(wn,n−1))=(A(wn0),A(wn1),…,A(wnn−1))\begin{align} \text{DFT}(a_0, a_1, \dots, a_{n-1}) &= (y_0, y_1, \dots, y_{n-1}) \\ &= (A(w_{n, 0}), A(w_{n, 1}), \dots, A(w_{n, n-1})) \\ &= (A(w_n^0), A(w_n^1), \dots, A(w_n^{n-1})) \end{align}

Аналогічно означується обернене дискретне перетворення Фур'є: обернене DFT значень многочлена (y0,y1,…,yn−1)(y_0, y_1, \dots, y_{n-1}) — це коефіцієнти многочлена (a0,a1,…,an−1)(a_0, a_1, \dots, a_{n-1}).

InverseDFT(y0,y1,…,yn−1)=(a0,a1,…,an−1)\text{InverseDFT}(y_0, y_1, \dots, y_{n-1}) = (a_0, a_1, \dots, a_{n-1})

Отже, якщо пряме DFT обчислює значення многочлена в точках, що є nn-ми коренями, то обернене DFT може відновити коефіцієнти многочлена за цими значеннями.

Застосування DFT: швидке множення многочленів​

Нехай маємо два многочлени AA і BB. Обчислимо DFT для кожного з них: DFT(A)\text{DFT}(A) і DFT(B)\text{DFT}(B).

Що станеться, якщо ми перемножимо ці многочлени? Очевидно, що в кожній точці значення просто перемножуються, тобто

(A⋅B)(x)=A(x)⋅B(x).(A \cdot B)(x) = A(x) \cdot B(x).

Це означає, що якщо ми перемножимо вектори DFT(A)\text{DFT}(A) і DFT(B)\text{DFT}(B) — перемноживши кожен елемент одного вектора на відповідний елемент іншого, — то ми отримаємо не що інше, як DFT многочлена DFT(A⋅B)\text{DFT}(A \cdot B):

DFT(A⋅B)=DFT(A)⋅DFT(B)\text{DFT}(A \cdot B) = \text{DFT}(A) \cdot \text{DFT}(B)

Нарешті, застосувавши обернене DFT, отримаємо:

A⋅B=InverseDFT(DFT(A)⋅DFT(B))A \cdot B = \text{InverseDFT}(\text{DFT}(A) \cdot \text{DFT}(B))

Справа під добутком двох DFT ми розуміємо поелементний добуток елементів векторів. Його можна обчислити за час O(n)O(n). Якщо ми вміємо обчислювати DFT та обернене DFT за O(nlog⁡n)O(n \log n), то ми можемо обчислити добуток двох многочленів (а отже, і двох довгих чисел) з тією самою часовою складністю.

Слід зауважити, що обидва многочлени повинні мати однаковий степінь. Інакше два результуючі вектори DFT матимуть різну довжину. Цього можна досягти, додавши коефіцієнти зі значенням 00.

А також, оскільки результатом добутку двох многочленів є многочлен степеня 2(n−1)2 (n - 1), нам доводиться подвоїти степені кожного многочлена (знову ж таки, доповнивши 00). З вектора з nn значень ми не можемо відновити бажаний многочлен із 2n−12n - 1 коефіцієнтами.

Швидке перетворення Фур'є​

Швидке перетворення Фур'є — це метод, який дозволяє обчислити DFT за час O(nlog⁡n)O(n \log n). Основна ідея FFT — застосувати «розділяй і володарюй». Ми ділимо вектор коефіцієнтів многочлена на два вектори, рекурсивно обчислюємо DFT для кожного з них і комбінуємо результати, щоб обчислити DFT повного многочлена.

Отже, нехай маємо многочлен A(x)A(x) степеня n−1n - 1, де nn — степінь 22 і n>1n > 1:

A(x)=a0x0+a1x1+⋯+an−1xn−1A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}

Поділимо його на два менші многочлени: один містить лише коефіцієнти на парних позиціях, а інший — коефіцієнти на непарних позиціях:

A0(x)=a0x0+a2x1+⋯+an−2xn2−1A1(x)=a1x0+a3x1+⋯+an−1xn2−1\begin{align} A_0(x) &= a_0 x^0 + a_2 x^1 + \dots + a_{n-2} x^{\frac{n}{2}-1} \\ A_1(x) &= a_1 x^0 + a_3 x^1 + \dots + a_{n-1} x^{\frac{n}{2}-1} \end{align}

Легко бачити, що

A(x)=A0(x2)+xA1(x2).A(x) = A_0(x^2) + x A_1(x^2).

Многочлени A0A_0 і A1A_1 мають удвічі менше коефіцієнтів, ніж многочлен AA. Якщо ми вміємо обчислити DFT(A)\text{DFT}(A) за лінійний час, використовуючи DFT(A0)\text{DFT}(A_0) і DFT(A1)\text{DFT}(A_1), то для часової складності отримаємо рекурентне співвідношення TDFT(n)=2TDFT(n2)+O(n)T_{\text{DFT}}(n) = 2 T_{\text{DFT}}\left(\frac{n}{2}\right) + O(n), що за основною теоремою дає TDFT(n)=O(nlog⁡n)T_{\text{DFT}}(n) = O(n \log n).

Навчимося, як цього досягти.

Припустимо, що ми обчислили вектори (yk0)k=0n/2−1=DFT(A0)\left(y_k^0\right)_{k=0}^{n/2-1} = \text{DFT}(A_0) і (yk1)k=0n/2−1=DFT(A1)\left(y_k^1\right)_{k=0}^{n/2-1} = \text{DFT}(A_1). Знайдемо вираз для (yk)k=0n−1=DFT(A)\left(y_k\right)_{k=0}^{n-1} = \text{DFT}(A).

Для перших n2\frac{n}{2} значень ми можемо просто скористатися зазначеним раніше рівнянням A(x)=A0(x2)+xA1(x2)A(x) = A_0(x^2) + x A_1(x^2):

yk=yk0+wnkyk1,k=0…n2−1.y_k = y_k^0 + w_n^k y_k^1, \quad k = 0 \dots \frac{n}{2} - 1.

Однак для других n2\frac{n}{2} значень потрібно знайти дещо інший вираз:

yk+n/2=A(wnk+n/2)=A0(wn2k+n)+wnk+n/2A1(wn2k+n)=A0(wn2kwnn)+wnkwnn/2A1(wn2kwnn)=A0(wn2k)−wnkA1(wn2k)=yk0−wnkyk1\begin{align} y_{k+n/2} &= A\left(w_n^{k+n/2}\right) \\ &= A_0\left(w_n^{2k+n}\right) + w_n^{k + n/2} A_1\left(w_n^{2k+n}\right) \\ &= A_0\left(w_n^{2k} w_n^n\right) + w_n^k w_n^{n/2} A_1\left(w_n^{2k} w_n^n\right) \\ &= A_0\left(w_n^{2k}\right) - w_n^k A_1\left(w_n^{2k}\right) \\ &= y_k^0 - w_n^k y_k^1 \end{align}

Тут ми знову скористалися A(x)=A0(x2)+xA1(x2)A(x) = A_0(x^2) + x A_1(x^2) і двома тотожностями wnn=1w_n^n = 1 та wnn/2=−1w_n^{n/2} = -1.

Отже, ми отримуємо бажані формули для обчислення всього вектора (yk)(y_k):

yk=yk0+wnkyk1,k=0…n2−1,yk+n/2=yk0−wnkyk1,k=0…n2−1.\begin{align} y_k &= y_k^0 + w_n^k y_k^1, &\quad k = 0 \dots \frac{n}{2} - 1, \\ y_{k+n/2} &= y_k^0 - w_n^k y_k^1, &\quad k = 0 \dots \frac{n}{2} - 1. \end{align}

(Цей шаблон a+ba + b та a−ba - b іноді називають метеликом.)

Таким чином, ми навчилися обчислювати DFT за час O(nlog⁡n)O(n \log n).

Обернене перетворення FFT​

Нехай дано вектор (y0,y1,…yn−1)(y_0, y_1, \dots y_{n-1}) — значення многочлена AA степеня n−1n - 1 у точках x=wnkx = w_n^k. Ми хочемо відновити коефіцієнти (a0,a1,…,an−1)(a_0, a_1, \dots, a_{n-1}) многочлена. Ця відома задача називається інтерполяцією, і для її розв'язання існують загальні алгоритми. Але в цьому окремому випадку (оскільки ми знаємо значення в точках, що є коренями з одиниці) можна отримати значно простіший алгоритм (фактично той самий, що й пряме FFT).

Запишемо DFT, згідно з його означенням, у матричній формі:

(wn0wn0wn0wn0⋯wn0wn0wn1wn2wn3⋯wnn−1wn0wn2wn4wn6⋯wn2(n−1)wn0wn3wn6wn9⋯wn3(n−1)⋮⋮⋮⋮⋱⋮wn0wnn−1wn2(n−1)wn3(n−1)⋯wn(n−1)(n−1))(a0a1a2a3⋮an−1)=(y0y1y2y3⋮yn−1)\begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix}

Ця матриця називається матрицею Вандермонда.

Таким чином, ми можемо обчислити вектор (a0,a1,…,an−1)(a_0, a_1, \dots, a_{n-1}), помноживши вектор (y0,y1,…yn−1)(y_0, y_1, \dots y_{n-1}) зліва на обернену матрицю:

(a0a1a2a3⋮an−1)=(wn0wn0wn0wn0⋯wn0wn0wn1wn2wn3⋯wnn−1wn0wn2wn4wn6⋯wn2(n−1)wn0wn3wn6wn9⋯wn3(n−1)⋮⋮⋮⋮⋱⋮wn0wnn−1wn2(n−1)wn3(n−1)⋯wn(n−1)(n−1))−1(y0y1y2y3⋮yn−1)\begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^1 & w_n^2 & w_n^3 & \cdots & w_n^{n-1} \\ w_n^0 & w_n^2 & w_n^4 & w_n^6 & \cdots & w_n^{2(n-1)} \\ w_n^0 & w_n^3 & w_n^6 & w_n^9 & \cdots & w_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{n-1} & w_n^{2(n-1)} & w_n^{3(n-1)} & \cdots & w_n^{(n-1)(n-1)} \end{pmatrix}^{-1} \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix}

Швидка перевірка може підтвердити, що обернена матриця має такий вигляд:

1n(wn0wn0wn0wn0⋯wn0wn0wn−1wn−2wn−3⋯wn−(n−1)wn0wn−2wn−4wn−6⋯wn−2(n−1)wn0wn−3wn−6wn−9⋯wn−3(n−1)⋮⋮⋮⋮⋱⋮wn0wn−(n−1)wn−2(n−1)wn−3(n−1)⋯wn−(n−1)(n−1))\frac{1}{n} \begin{pmatrix} w_n^0 & w_n^0 & w_n^0 & w_n^0 & \cdots & w_n^0 \\ w_n^0 & w_n^{-1} & w_n^{-2} & w_n^{-3} & \cdots & w_n^{-(n-1)} \\ w_n^0 & w_n^{-2} & w_n^{-4} & w_n^{-6} & \cdots & w_n^{-2(n-1)} \\ w_n^0 & w_n^{-3} & w_n^{-6} & w_n^{-9} & \cdots & w_n^{-3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w_n^0 & w_n^{-(n-1)} & w_n^{-2(n-1)} & w_n^{-3(n-1)} & \cdots & w_n^{-(n-1)(n-1)} \end{pmatrix}

Таким чином, ми отримуємо формулу:

ak=1n∑j=0n−1yjwn−kja_k = \frac{1}{n} \sum_{j=0}^{n-1} y_j w_n^{-k j}

Порівнюючи це з формулою для yky_k

yk=∑j=0n−1ajwnkj,y_k = \sum_{j=0}^{n-1} a_j w_n^{k j},

ми помічаємо, що ці задачі майже однакові, тож коефіцієнти aka_k можна знайти тим самим алгоритмом «розділяй і володарюй», що й пряме FFT, тільки замість wnkw_n^k нам доведеться використати wn−kw_n^{-k}, а наприкінці потрібно поділити отримані коефіцієнти на nn.

Таким чином, обчислення оберненого DFT майже таке саме, як і обчислення прямого DFT, і його також можна виконати за час O(nlog⁡n)O(n \log n).

Реалізація​

Тут ми наведемо просту рекурсивну реалізацію FFT та оберненого FFT, обидва в одній функції, оскільки різниця між прямим та оберненим FFT настільки мінімальна. Для зберігання комплексних чисел ми використовуємо тип complex зі стандартної бібліотеки C++ STL.

using cd = complex<double>;
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;
    }
}

Функції передається вектор коефіцієнтів, і вона обчислює DFT або обернене DFT та знову зберігає результат у цьому ж векторі. Аргумент invert\text{invert} показує, чи слід обчислювати пряме, чи обернене DFT. Усередині функції ми спершу перевіряємо, чи довжина вектора дорівнює одиниці, і якщо так, то нам нічого не потрібно робити. Інакше ми ділимо вектор aa на два вектори a0a0 і a1a1 та рекурсивно обчислюємо DFT для обох. Потім ми ініціалізуємо значення wnwn та змінну ww, яка міститиме поточний степінь wnwn. Потім за наведеними вище формулами обчислюються значення результуючого DFT.

Якщо прапорець invert\text{invert} встановлено, то ми замінюємо wnwn на wn−1wn^{-1}, і кожне зі значень результату ділиться на 22 (оскільки це робитиметься на кожному рівні рекурсії, у підсумку кінцеві значення поділяться на nn).

Використовуючи цю функцію, ми можемо створити функцію для множення двох многочленів:

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(fa, false);
    fft(fb, false);
    for (int i = 0; i < n; i++)
        fa[i] *= fb[i];
    fft(fa, true);

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

Ця функція працює з многочленами з цілими коефіцієнтами, проте її можна налаштувати й для роботи з іншими типами. Оскільки під час роботи з комплексними числами виникає певна похибка, нам потрібно наприкінці заокруглити отримані коефіцієнти.

Нарешті, функція для множення двох довгих чисел практично не відрізняється від функції для множення многочленів. Єдине, що нам потрібно зробити після цього, — це нормалізувати число:

    int carry = 0;
    for (int i = 0; i < n; i++)
        result[i] += carry;
        carry = result[i] / 10;
        result[i] %= 10;
    }

Оскільки довжина добутку двох чисел ніколи не перевищує сумарної довжини обох чисел, розміру вектора достатньо, щоб виконати всі операції перенесення.

Покращена реалізація: обчислення на місці​

Щоб підвищити ефективність, ми перейдемо від рекурсивної реалізації до ітеративної. У наведеній вище рекурсивній реалізації ми явно розділяли вектор aa на два вектори — елементи на парних позиціях присвоювалися одному тимчасовому вектору, а елементи на непарних позиціях — іншому. Однак якщо ми переупорядкуємо елементи певним чином, нам не потрібно створювати ці тимчасові вектори (тобто всі обчислення можна виконати «на місці», прямо в самому векторі AA).

Зауважте, що на першому рівні рекурсії елементи, у яких молодший біт позиції дорівнював нулю, присвоювалися вектору a0a_0, а ті, у яких молодший біт позиції дорівнював одиниці, — вектору a1a_1. На другому рівні рекурсії відбувається те саме, але вже з другим молодшим бітом, і т. д. Тому якщо ми обернемо біти позиції кожного коефіцієнта й відсортуємо їх за цими оберненими значеннями, ми отримаємо бажаний порядок (його називають перестановкою з оберненням бітів, bit-reversal permutation).

Наприклад, бажаний порядок для n=8n = 8 має вигляд:

a={[(a0,a4),(a2,a6)],[(a1,a5),(a3,a7)]}a = \bigg\{ \Big[ (a_0, a_4), (a_2, a_6) \Big], \Big[ (a_1, a_5), (a_3, a_7) \Big] \bigg\}

Дійсно, на першому рівні рекурсії (обведеному фігурними дужками) вектор ділиться на дві частини [a0,a2,a4,a6][a_0, a_2, a_4, a_6] та [a1,a3,a5,a7][a_1, a_3, a_5, a_7]. Як ми бачимо, у перестановці з оберненням бітів це відповідає простому поділу вектора на дві половини: перші n2\frac{n}{2} елементів та останні n2\frac{n}{2} елементів. Потім для кожної половини відбувається рекурсивний виклик. Нехай отримане для кожної з них DFT повертається на місце самих елементів (тобто в першу половину та другу половину вектора aa відповідно).

a={[y00,y10,y20,y30],[y01,y11,y21,y31]}a = \bigg\{ \Big[y_0^0, y_1^0, y_2^0, y_3^0\Big], \Big[y_0^1, y_1^1, y_2^1, y_3^1 \Big] \bigg\}

Тепер ми хочемо об'єднати два DFT в одне для повного вектора. Порядок елементів ідеальний, і об'єднання ми також можемо виконати прямо в цьому векторі. Можемо взяти елементи y00y_0^0 та y01y_0^1 і виконати перетворення метелика. Місце двох отриманих значень те саме, що й місце двох початкових значень, тож ми отримуємо:

a={[y00+wn0y01,y10,y20,y30],[y00−wn0y01,y11,y21,y31]}a = \bigg\{ \Big[y_0^0 + w_n^0 y_0^1, y_1^0, y_2^0, y_3^0\Big], \Big[y_0^0 - w_n^0 y_0^1, y_1^1, y_2^1, y_3^1\Big] \bigg\}

Аналогічно можемо обчислити перетворення метелика для y10y_1^0 та y11y_1^1 і покласти результати на їхнє місце, і так далі. У результаті отримуємо:

a={[y00+wn0y01,y10+wn1y11,y20+wn2y21,y30+wn3y31],[y00−wn0y01,y10−wn1y11,y20−wn2y21,y30−wn3y31]}a = \bigg\{ \Big[y_0^0 + w_n^0 y_0^1, y_1^0 + w_n^1 y_1^1, y_2^0 + w_n^2 y_2^1, y_3^0 + w_n^3 y_3^1\Big], \Big[y_0^0 - w_n^0 y_0^1, y_1^0 - w_n^1 y_1^1, y_2^0 - w_n^2 y_2^1, y_3^0 - w_n^3 y_3^1\Big] \bigg\}

Таким чином, ми обчислили потрібне DFT з вектора aa.

Тут ми описали процес обчислення DFT лише на першому рівні рекурсії, але те саме, очевидно, працює й для всіх інших рівнів. Отже, застосувавши перестановку з оберненням бітів, ми можемо обчислити DFT на місці, без жодної додаткової пам'яті.

Це додатково дозволяє нам позбутися рекурсії. Ми просто починаємо з найнижчого рівня, тобто ділимо вектор на пари та застосовуємо до них перетворення метелика. У результаті отримуємо вектор aa з виконаною роботою останнього рівня. На наступному кроці ми ділимо вектор на вектори розміру 44 і знову застосовуємо перетворення метелика, що дає нам DFT для кожного блоку розміру 44. І так далі. Нарешті, на останньому кроці ми отримали результат DFT обох половин aa, і, застосувавши перетворення метелика, отримуємо DFT для повного вектора aa.

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

int reverse(int num, int lg_n) {
    int res = 0;
    for (int i = 0; i < lg_n; i++) {
        if (num & (1 << i))
            res |= 1 << (lg_n - 1 - i);
    }
    return res;
}

void fft(vector<cd> & a, bool invert) {
    int n = a.size();
    int lg_n = 0;
    while ((1 << lg_n) < n)
        lg_n++;

    for (int i = 0; i < n; i++) {
        if (i < reverse(i, lg_n))
            swap(a[i], a[reverse(i, lg_n)]);
    }

    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;
    }
}

Спершу ми застосовуємо перестановку з оберненням бітів, обмінюючи кожен елемент з елементом на оберненій позиції. Потім на log⁡n−1\log n - 1 станах алгоритму ми обчислюємо DFT для кожного блоку відповідного розміру len\text{len}. Для всіх цих блоків ми маємо той самий корінь з одиниці wlen\text{wlen}. Ми проходимо всі блоки й виконуємо для кожного з них перетворення метелика.

Ми можемо ще додатково оптимізувати обернення бітів. У попередній реалізації ми проходили всі біти індексу й створювали побітово обернений індекс. Однак біти можна обернути в інший спосіб.

Припустимо, що jj уже містить обернення ii. Тоді, щоб перейти до i+1i + 1, нам потрібно збільшити ii, а також збільшити jj, але в «оберненій» системі числення. Додавання одиниці у звичайній двійковій системі еквівалентне перетворенню всіх хвостових одиниць на нулі та перетворенню нуля безпосередньо перед ними на одиницю. Так само в «оберненій» системі числення ми перетворюємо всі провідні одиниці, а також наступний нуль.

Таким чином, ми отримуємо таку реалізацію:

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

void fft(vector<cd> & a, bool invert) {
    int n = 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)
            swap(a[i], a[j]);
    }

    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;
    }
}

Додатково ми можемо обчислити перестановку з оберненням бітів заздалегідь. Це особливо корисно, коли розмір nn однаковий для всіх викликів. Але навіть коли ми маємо лише три виклики (які необхідні для множення двох многочленів), ефект помітний. Також ми можемо заздалегідь обчислити всі корені з одиниці та їхні степені.

Теоретико-числове перетворення​

Тепер трохи змінимо мету. Ми все ще хочемо перемножити два многочлени за час O(nlog⁡n)O(n \log n), але цього разу хочемо обчислити коефіцієнти за модулем деякого простого числа pp. Звісно, для цього завдання ми можемо використати звичайне DFT і застосувати операцію взяття за модулем до результату. Однак це може призвести до похибок заокруглення, особливо коли йдеться про великі числа. Теоретико-числове перетворення (NTT) має ту перевагу, що воно працює лише з цілими числами, тому результат гарантовано буде правильним.

Дискретне перетворення Фур'є базується на комплексних числах і nn-х коренях з одиниці. Щоб ефективно його обчислити, ми активно використовуємо властивості коренів (наприклад, те, що існує один корінь, який породжує всі інші корені піднесенням до степеня).

Але ті самі властивості виконуються й для nn-х коренів з одиниці в модульній арифметиці. nn-й корінь з одиниці над простим полем — це таке число wnw_n, що задовольняє:

(wn)n=1(modp),(wn)k≠1(modp),1≤k<n.\begin{align} (w_n)^n &= 1 \pmod{p}, \\ (w_n)^k &\ne 1 \pmod{p}, \quad 1 \le k < n. \end{align}

Інші n−1n-1 коренів можна отримати як степені кореня wnw_n.

Щоб застосувати це в алгоритмі швидкого перетворення Фур'є, нам потрібно, щоб корінь існував для деякого nn, що є степенем 22, а також для всіх менших степенів. Можемо помітити таку цікаву властивість:

(wn2)m=wnn=1(modp),з m=n2(wn2)k=wn2k≠1(modp),1≤k<m.\begin{align} (w_n^2)^m = w_n^n &= 1 \pmod{p}, \quad \text{з } m = \frac{n}{2}\\ (w_n^2)^k = w_n^{2k} &\ne 1 \pmod{p}, \quad 1 \le k < m. \end{align}

Отже, якщо wnw_n — це nn-й корінь з одиниці, то wn2w_n^2 — це n2\frac{n}{2}-й корінь з одиниці. А отже, для всіх менших степенів двійки існують корені потрібного степеня, і їх можна обчислити за допомогою wnw_n.

Для обчислення оберненого DFT нам потрібен обернений елемент wn−1w_n^{-1} до wnw_n. Але для простого модуля обернений завжди існує.

Таким чином, усі властивості, які нам потрібні від комплексних коренів, доступні й у модульній арифметиці за умови, що ми маємо достатньо великий модуль pp, для якого існує nn-й корінь з одиниці.

Наприклад, ми можемо взяти такі значення: модуль p=7340033p = 7340033, w220=5w_{2^{20}} = 5. Якщо цього модуля недостатньо, нам потрібно знайти іншу пару. Ми можемо скористатися тим фактом, що для модулів вигляду p=c2k+1p = c 2^k + 1 (і pp — просте) завжди існує 2k2^k-й корінь з одиниці. Можна показати, що gcg^c — це такий 2k2^k-й корінь з одиниці, де gg — первісний корінь числа pp.

const int mod = 7340033;
const int root = 5;
const int root_1 = 4404020;
const int root_pw = 1 << 20;

void fft(vector<int> & a, bool invert) {
    int n = 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)
            swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        int wlen = invert ? root_1 : root;
        for (int i = len; i < root_pw; i <<= 1)
            wlen = (int)(1LL * wlen * wlen % mod);

        for (int i = 0; i < n; i += len) {
            int w = 1;
            for (int j = 0; j < len / 2; j++) {
                int u = a[i+j], v = (int)(1LL * a[i+j+len/2] * w % mod);
                a[i+j] = u + v < mod ? u + v : u + v - mod;
                a[i+j+len/2] = u - v >= 0 ? u - v : u - v + mod;
                w = (int)(1LL * w * wlen % mod);
            }
        }
    }

    if (invert) {
        int n_1 = inverse(n, mod);
        for (int & x : a)
            x = (int)(1LL * x * n_1 % mod);
    }
}

Тут функція inverse обчислює обернений елемент за модулем (див. Обернений елемент за модулем). Константи mod, root, root_pw визначають модуль і корінь, а root_1 — це обернений до root за модулем mod.

На практиці ця реалізація повільніша за реалізацію з комплексними числами (через величезну кількість операцій за модулем), але має деякі переваги, як-от менше використання пам'яті та відсутність похибок заокруглення.

Множення з довільним модулем​

Тут ми хочемо досягти тієї самої мети, що й у попередньому розділі. Перемножити два многочлени A(x)A(x) і B(x)B(x) та обчислити коефіцієнти за модулем деякого числа MM. Теоретико-числове перетворення працює лише для певних простих чисел. А як бути у випадку, коли модуль не має потрібного вигляду?

Одним із варіантів було б виконати кілька теоретико-числових перетворень з різними простими числами вигляду c2k+1c 2^k + 1, а потім застосувати китайську теорему про остачі, щоб обчислити остаточні коефіцієнти.

Інший варіант — розподілити кожен із многочленів A(x)A(x) і B(x)B(x) на два менші многочлени

A(x)=A1(x)+A2(x)⋅CB(x)=B1(x)+B2(x)⋅C\begin{align} A(x) &= A_1(x) + A_2(x) \cdot C \\ B(x) &= B_1(x) + B_2(x) \cdot C \end{align}

з C≈MC \approx \sqrt{M}.

Тоді добуток A(x)A(x) і B(x)B(x) можна подати як:

A(x)⋅B(x)=A1(x)⋅B1(x)+(A1(x)⋅B2(x)+A2(x)⋅B1(x))⋅C+(A2(x)⋅B2(x))⋅C2A(x) \cdot B(x) = A_1(x) \cdot B_1(x) + \left(A_1(x) \cdot B_2(x) + A_2(x) \cdot B_1(x)\right)\cdot C + \left(A_2(x) \cdot B_2(x)\right)\cdot C^2

Многочлени A1(x)A_1(x), A2(x)A_2(x), B1(x)B_1(x) і B2(x)B_2(x) містять лише коефіцієнти, менші за M\sqrt{M}, тому коефіцієнти всіх добутків, що з'являються, менші за M⋅nM \cdot n, що зазвичай достатньо мало, щоб обробити це типовими типами з рухомою комою.

Отже, цей підхід вимагає обчислення добутків многочленів із меншими коефіцієнтами (за допомогою звичайних FFT та оберненого FFT), а потім вихідний добуток можна відновити за допомогою модульного додавання та множення за час O(n)O(n).

Застосування​

DFT можна використовувати у величезній різноманітності інших задач, які на перший погляд не мають нічого спільного з множенням многочленів.

Усі можливі суми​

Нам дано два масиви a[]a[] і b[]b[]. Нам потрібно знайти всі можливі суми a[i]+b[j]a[i] + b[j] і для кожної суми порахувати, як часто вона трапляється.

Наприклад, для a=[1, 2, 3]a = [1,~ 2,~ 3] і b=[2, 4]b = [2,~ 4] отримуємо: суму 33 можна отримати 11 способом, суму 44 також 11 способом, 55 — 22 способами, 66 — 11, 77 — 11.

Будуємо для масивів aa і bb два многочлени AA і BB. Числа масиву виступатимуть як показники степеня в многочлені (a[i]⇒xa[i]a[i] \Rightarrow x^{a[i]}); а коефіцієнтом при цьому члені буде те, як часто число трапляється в масиві.

Тоді, перемноживши ці два многочлени за час O(nlog⁡n)O(n \log n), ми отримуємо многочлен CC, де показники степеня скажуть нам, які суми можна отримати, а коефіцієнти — як часто. Продемонструймо це на прикладі:

(1x1+1x2+1x3)(1x2+1x4)=1x3+1x4+2x5+1x6+1x7(1 x^1 + 1 x^2 + 1 x^3) (1 x^2 + 1 x^4) = 1 x^3 + 1 x^4 + 2 x^5 + 1 x^6 + 1 x^7

Усі можливі скалярні добутки​

Нам дано два масиви a[]a[] і b[]b[] довжини nn. Нам потрібно обчислити добутки aa з кожним циклічним зсувом bb.

Ми генеруємо два нові масиви розміру 2n2n: розвертаємо aa і дописуємо до нього nn нулів. А bb просто дописуємо сам до себе. Коли ми перемножуємо ці два масиви як многочлени й дивимося на коефіцієнти c[n−1], c[n], …, c[2n−2]c[n-1],~ c[n],~ \dots,~ c[2n-2] добутку cc, отримуємо:

c[k]=∑i+j=ka[i]b[j]c[k] = \sum_{i+j=k} a[i] b[j]

І оскільки всі елементи a[i]=0a[i] = 0 при i≥ni \ge n:

c[k]=∑i=0n−1a[i]b[k−i]c[k] = \sum_{i=0}^{n-1} a[i] b[k-i]

Легко бачити, що ця сума — це просто скалярний добуток вектора aa з (k−(n−1))(k - (n - 1))-м циклічним зсувом bb уліво. Таким чином, ці коефіцієнти є відповіддю на задачу, і ми все ж змогли отримати її за час O(nlog⁡n)O(n \log n). Зауважте, що c[2n−1]c[2n-1] також дає нам nn-й циклічний зсув, але це те саме, що й 00-й циклічний зсув, тож нам не потрібно враховувати його окремо у нашій відповіді.

Дві смуги​

Нам дано дві булеві смуги (циклічні масиви зі значень 00 і 11) aa і bb. Ми хочемо знайти всі способи прикласти першу смугу до другої так, щоб у жодній позиції 11 першої смуги не стояла поряд із 11 другої смуги.

Ця задача насправді не дуже відрізняється від попередньої. Прикладання двох смуг просто означає, що ми виконуємо циклічний зсув другого масиву, і ми можемо прикласти дві смуги, якщо скалярний добуток двох масивів дорівнює 00.

Зіставлення рядків​

Нам дано два рядки: текст TT і взірець PP, що складаються з малих літер. Нам потрібно обчислити всі входження взірця в текст.

Ми створюємо многочлен для кожного рядка (T[i]T[i] і P[I]P[I] — числа між 00 і 2525, що відповідають 2626 літерам алфавіту):

A(x)=a0x0+a1x1+⋯+an−1xn−1,n=∣T∣A(x) = a_0 x^0 + a_1 x^1 + \dots + a_{n-1} x^{n-1}, \quad n = |T|

з

ai=cos⁡(αi)+isin⁡(αi),αi=2πT[i]26.a_i = \cos(\alpha_i) + i \sin(\alpha_i), \quad \alpha_i = \frac{2 \pi T[i]}{26}.

І

B(x)=b0x0+b1x1+⋯+bm−1xm−1,m=∣P∣B(x) = b_0 x^0 + b_1 x^1 + \dots + b_{m-1} x^{m-1}, \quad m = |P|

з

bi=cos⁡(βi)−isin⁡(βi),βi=2πP[m−i−1]26.b_i = \cos(\beta_i) - i \sin(\beta_i), \quad \beta_i = \frac{2 \pi P[m-i-1]}{26}.

Зауважте, що вираз P[m−i−1]P[m-i-1] явно розвертає взірець.

(m−1+i)(m-1+i)-ті коефіцієнти добутку двох многочленів C(x)=A(x)⋅B(x)C(x) = A(x) \cdot B(x) скажуть нам, чи з'являється взірець у тексті в позиції ii.

cm−1+i=∑j=0m−1ai+j⋅bm−1−j=∑j=0m−1(cos⁡(αi+j)+isin⁡(αi+j))⋅(cos⁡(βj)−isin⁡(βj))c_{m-1+i} = \sum_{j = 0}^{m-1} a_{i+j} \cdot b_{m-1-j} = \sum_{j=0}^{m-1} \left(\cos(\alpha_{i+j}) + i \sin(\alpha_{i+j})\right) \cdot \left(\cos(\beta_j) - i \sin(\beta_j)\right)

з αi+j=2πT[i+j]26\alpha_{i+j} = \frac{2 \pi T[i+j]}{26} та βj=2πP[j]26\beta_j = \frac{2 \pi P[j]}{26}

Якщо є збіг, то T[i+j]=P[j]T[i+j] = P[j], а отже, αi+j=βj\alpha_{i+j} = \beta_j. Це дає (використовуючи основну тригонометричну тотожність):

cm−1+i=∑j=0m−1(cos⁡(αi+j)+isin⁡(αi+j))⋅(cos⁡(αi+j)−isin⁡(αi+j))=∑j=0m−1cos⁡(αi+j)2+sin⁡(αi+j)2=∑j=0m−11=m\begin{align} c_{m-1+i} &= \sum_{j = 0}^{m-1} \left(\cos(\alpha_{i+j}) + i \sin(\alpha_{i+j})\right) \cdot \left(\cos(\alpha_{i+j}) - i \sin(\alpha_{i+j})\right) \\ &= \sum_{j = 0}^{m-1} \cos(\alpha_{i+j})^2 + \sin(\alpha_{i+j})^2 = \sum_{j = 0}^{m-1} 1 = m \end{align}

Якщо збігу немає, то принаймні один символ відрізняється, через що один із добутків ai+1⋅bm−1−ja_{i+1} \cdot b_{m-1-j} не дорівнює 11, що призводить до коефіцієнта cm−1+i≠mc_{m-1+i} \ne m.

Зіставлення рядків із символами підстановки​

Це розширення попередньої задачі. Цього разу ми дозволяємо, щоб взірець містив символ підстановки \*\*, який може зіставитися з будь-якою можливою літерою. Наприклад, взірець a∗ca*c з'являється в тексті abccaaccabccaacc рівно у трьох позиціях: в індексі 00, індексі 44 та індексі 55.

Ми створюємо точнісінько ті самі многочлени, за винятком того, що встановлюємо bi=0b_i = 0, якщо P[m−i−1]=∗P[m-i-1] = *. Якщо xx — кількість символів підстановки в PP, то ми матимемо збіг PP у TT в індексі ii, якщо cm−1+i=m−xc_{m-1+i} = m - x.

Задачі для практики​

Відеоматеріали​