Быстрое преобразование Фурье

Материал из Algocode wiki
Перейти к: навигация, поиск

То что здесь написано не претендует на звание самого понятного и полного конспекта. На emaxx довольно хорошо написанно про FFT, поэтому можно почитать и там. Этот конспект нужен скорее как предисловие к длинной арифметике и чтобы показать простую оптимальную реализацию fft, которую реально запомнить. Если есть желающие перезаписать эту страницу, смело это делайте.

Умножение многочленов

Умножать многочлены — это круто, но долго. Самый тривиальный алгоритм делает это за $O(n^2)$. Есть конечно Карацуба, который справляется с этим чуть быстрее, но тоже не идеально. Но есть алгоритм быстрого преобразования Фурье, с помощью которого это достигается за $O(n \log n)$.

В классическом понимании для умножения многочленов надо перемножить их коэффициенты. Но можно подойти к этому с другой стороны. Посчитаем значения многочленов в $2n$ различных точках. Замечаем, что значение их произведения в каждой из этих точек будет равно произведению их значений в каждой из этих точек. А т.к. различные многочлены степени $n$ имеют различные значения в $n + 1$ точке, то если по значения многочлена в $n$ точках мы сможем восстановить многочлен, то так мы найдём произведение.

Прямое преобразование Фурье

Надо найти значения многочлене степени $n$ в $n$ точках. Для удобства увеличим $n$ так, чтобы $n$ стал точной степенью 2.

Пусть есть многочлен $A = a_0 + a_1 * x + a_2 * x^2 + \ldots + a_{n - 1} * x^{n - 1}$. Заметим, что $$A = (a_0 + a_2 * x^2 + x_4 * x^4 + \ldots + a_{n - 2} * x^{n - 2}) + x * (a_1 + a_3 * x^2 + \ldots + a_{n - 1} * x^{n - 2})$$

Мы хотим посчитать его значения в $x_0, x_1, \ldots x_{n - 1}$. Тогда заметим, что если бы множество из $x_0^2, x_1^2, \ldots x_{n - 1}^2$ состояло бы из $n / 2$ чисел, то мы бы посчитали бы за $O(n / 2 \log(n / 2))$ значения $(a_0 + a_2 * x^2 + x_4 * x^4 + \ldots + a_{n - 2} * x^{n - 2})$ в этих $n / 2$ точках, далее за $O(n / 2 \log(n / 2))$ посчитали бы значения $(a_1 + a_3 * x^2 + \ldots + a_{n - 1} * x^{n - 2})$ в этих же $n / 2$ точках и в конце за $O(n)$ построили бы нормальный ответ.

Такие $x_0, x_1, \ldots, x_{n - 1}$ существуют. Для этого нам потребуется первообразный $\sqrt[n]{1}$. Первообразным корнем $n$-й степени из 1 называется такое число, что оно во всех степенях меньших или равных $n$ принимает разные значения, а в степени $n$ оно равно 1. Вспомниим комплексные числа. Возьмём число $z = (\cos (\frac{2 * \pi}{n}),\: \sin (\frac{2 * \pi}{n}))$. Замечаем что $z$ соответствует вектору, имеющему длину 1 и угол $\frac{2 * \pi}{n}$. Тогда $z^k$ соответствует вектору, имеющему длину 1 и угол $k * \frac{2 * \pi}{n}$. Тогда $z^n$ равен $1$ и все числа $z^1, z^2, \ldots z^n$ различны.

Тогда возьмём $x_0 = z^0, x_1 = z^1, x_2 = z^2, \ldots, x_{n - 1} = z^{n - 1}$. Замечаем, что $x_i^2 = z^{2i}$. Тогда угол вектора, соответствующего $x_i$ делится на $2 * \frac{2 * \pi}{n}$. А таких углов всего $n / 2$.

Тогда быстрое преобразование устроено так:

Сначала многочлен делится на два многочлена размера $n / 2$.

$$A = (a_0 + a_2 * x^2 + x_4 * x^4 + \ldots + a_{n - 2} * x^{n - 2}) + x * (a_1 + a_3 * x^2 + \ldots + a_{n - 1} * x^{n - 2})$$

Для каждого из них считаются значения в разных чётных степенях $z = (\cos (\frac{2 * \pi}{n}),\: \sin (\frac{2 * \pi}{n}))$. После чего на основе этого считаются значения $A$ во всех степенях $z$.

Обратное преобразование Фурье

Нам надо по значениям многочлена степени $n - 1$ в $n$ различных точках, равных степеням $z = (\cos (\frac{2 * \pi}{n}),\: \sin (\frac{2 * \pi}{n}))$, восстановить многочлен.

Будем использовать прямое FFT к значения многочлена, но вместо $z$ возьмём $\overline{z} = (\cos (-\frac{2 * \pi}{n}),\: \sin (-\frac{2 * \pi}{n}))$. В конце каждый член надо будет поделить на $n$. Так мы получим ответ.

Доказательство правильности этого можете загуглить или прочитать на том же emaxx.

Реализация

Единственное полезное, что есть на этой странице.

 1 // Этот код ни разу не компилировался и создан для ознакомления
 2 void fft(
 3     vector<complex<double>> &s, // Массив коэффициентов, передаётся по ссылке
 4     vector<complex<double>> &res, // Массив результата
 5     int n, // Размер массива, обязательно степень 2
 6     complex<double> x, // Комплексный первообразный корень из 1
 7         // В случае прямого преобразования это complex<double>(cos(phi * 2 / n), sin(phi * 2 / n))
 8         // В случае обратного преобразования это complex<double>(cos(phi * 2 / n), -sin(phi * 2 / n))
 9     int bs = 0, // Начало, с которого надо рассматривать коэффициенты из s
10     int step = 1, // Шаг, с которым надо рассматривать коэффициенты из s
11     int rs = 0, // Начало, с которого надо записывать результат в массив res
12 ) {
13     if (n == 1) {
14         res[rs] = s[bs];
15         return;
16     }
17     fft(s, res, n / 2, x * x, bs, bstep * 2, rs);
18     // Теперь в res[rs:rs + n / 2] записан результат fft для каждого второго коэффициента
19     fft(s, res, n / 2, x * x, bs + bstep, bstep * 2, rs + n / 2);
20     // Теперь в res[rs + n / 2:rs + n] записан результат fft для оставшихся коэффициентов
21     complex<double> c = 1;
22     // Теперь надо сделать новый res из старого
23     for (int i = rs; i < rs + n / 2; i++) {
24         auto a = res[i];
25         auto b = res[i + n / 2];
26         res[i] = a + b * c;
27         res[i + n / 2] = a - b * c; // Это можно написать оптимальнее без переменной b, но так понятнее
28         c *= x;
29     }
30 }

Применения

FFT часто применяется для быстрого умножения чисел. Вот тут написано как написать оставшуюся часть длинной арифметики.



Автор конспекта: Филипп Грибов

По всем вопросам пишите в telegram @grphil