Как сделать преобразование фурье

Добавил пользователь Валентин П.
Обновлено: 04.10.2024

Основная задача, которая решается при помощи быстрого преобразования Фурье (Fast Fourier Transform, FFT) — это умножение многочленов за время $O(n \log n)$.

Тривиально (по определению) многочлены степеней $m$ и $n$ умножаются за время $O(nm)$. Долгое время считалось, что быстрее это сделать невозможно, А.Н. Колмогоров даже принимал попытки доказать нижнюю границу, пока в 1960 году А. Карацуба не придумал способ умножать многочлены степени $n$ за время $O(n^<\log_<2>3>)$.

Как известно, многочлен степени строго меньше $n$ однозначно определяется своими значениями в $n$ (вообще говоря комплексных) точках. Действительно, если есть два различных многочлена с одинаковыми значениями в $n$ точках, то их разность имеет $n$ комплексных корней, причём она является ненулевым многочленом степени строго меньше $n$, что противоречит основной теореме алгебры. С другой стороны, интерполяционный многочлен Лагранжа в явном виде предъявляет многочлен степени строго меньше $n$ по значениям в $n$ точках.

Таким образом, многочлены можно хранить не в виде вектора коэффициентов, а в виде набора значений в некоторых точках. Над многочленами в таком виде очень удобно производить арифметические операции, в том числе умножать их за время $O(n)$ (нужно просто перемножить значения в соответствующих точках). С другой стороны, непонятно, как считать значения в других точках, да и знать сами коэффициенты бывает полезно. А воостанавливать коэффициенты по значениям в некотором наборе точек сложно, тот же интерполяционный многочлен Лагранжа вычисляется за время $O(n^2)$. Да и многочлены нам обычно задаются в форме вектора коэффициентов, получить значеня в $n$ произвольных точках вряд ли можно быстрее, чем за $O(n^2)$.

Хитрость FFT в том, что точки, в которых считаются значения многочлена, выбираются отнюдь не произвольным образом.

Описание

Итак, FFT преобразует вектор $\langle a_, a_, \ldots, a_ \rangle$ в вектор $\langle b_, b_, \ldots, b_ \rangle$, где $b_ = \sum_^ a_ e^<2 \pi i \frac>$, $n=2^$, иначе говоря, преобразует вектор коэффициентов многочлена степени $n-1$ в набор его значений в точках $\omega_ = e^<2 \pi i \frac>$.

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

Что делать с многочленами других степеней?
Почему именно эти точки?

С первым вопросом всё просто: нужно дополнять коэффициенты нулями до ближайшей степени двойки.

Со вторым же вопросом дело чуть хитрее. $\omega_$ - это комплексные корни из $1$ $n$-й степени, то есть $\omega_^n = 1$. У них есть замечательные свойства: $\omega_ = \omega_^$, $\omega_ \omega_ = \omega_$.

Разделяй и властвуй

$A_$ и $A_$ - многочлены степени $\frac = 2^$, к ним можно применить FFT и получить набор значений в точках $e^> = e^> = \omega_$.

$A(\omega_) = A_(\omega_^) + \omega_ A_(\omega_^) = A_(\omega_) + \omega_ A_(\omega_)$. Таким образом, зная FFT для $A_$ и $A_$, можно вычислить FFT для $A$.

Осталось определить базу рекурсии — $n=1$. Для этого нужно посчитать значение константного многочлена $A=a_$ в точке $\omega_ = 1$. Это значение, очевидно, равно $a_$.

Таким образом, алгоритм работает за время $T(n) = 2 T(\frac) + O(n)$. По мастер-теореме $T(n) = O(n \log n)$.

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

Хорошо, мы научились по многочлену вычислять его значение в $n$ особых точках за $O(n \log n)$, потом мы можем перемножить значения и получить FFT от произведения многочленов. Но мы пока не умеем восстанавливать многочлен по его FFT.

Можно рассмотреть FFT как линейное преобразование:

Возведём матрицу преобразования в квадрат: $R_ = \sum_^ \omega_^ \omega_^ = \sum_^ \omega_ = \sum_^ \omega_^$

Если $\omega_ = 1$, то $R_ = n$. В противном случае можно посчитать сумму геометрической прогрессии и получить $0$.

Таким образом, обратное преобразование выглядит следующим образом:

Применить прямое преобразование
Разделить все элементы на $n$
Развернуть массив без первого элемента

Оптимизации

В общем-то, это всё; однако данная реализация работает не очень быстро и потребляет $O(n \log n)$ памяти.

Что мы глобально делаем:

Переставляем коэффициенты
Делаем рекурсивные запуски
Склеиваем результаты

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

Как переставляются индексы: в левую часть идут чётные индексы, в правую — нечётные, потом к каждой из половин рекурсивно примерняется то же правило. Если рассмотреть битовую запись числа, то мы сначала сортируем по младшему биту, потом “отрезаем” его и делаем дальше то же самое. Нетрудно догадаться, что мы на самом деле сортируем индексы по реверснутой битовой записи; строгое доказательство проводится мат.индукцией по номеру уровня.

Реверснутую битовую запись всех чисел от $0$ до $2^-1$ можно предподсчитать с помощью ДП, воспользовавшись следующей идеей: $rev(mask) = rev(mask \oplus 2^) \oplus 2^$.

Также можно заранее предподсчитать комплексные корни из $1$. Тут есть тонкий момент: синусы и косинусы весьма медленные, в то время как если считать $\omega_ = \omega \omega_$, то может набежать большая погрешность. Можно использовать промежуточные варианты, например, честно посчитать первые $2^>$ значений, а также значения по индексам, кратным $2^>$, а все остальные корни разложить в произведение двух уже посчитанных.

Наконец, можно уменьшить количество умножений комплексных чисел. Вспомним, что

Их можно считать одновременно, это уменьшит количество умножений вдвое.

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

Размер применяемого FFT должен быть строго больше, чем степень произведения многочленов.

В задачах часто многочлены имеют целочисленные коэффициенты, причём неотрицательные. Понятно, что при таких условиях коэффициенты у произведения будут тоже целыми и неотрицательными.

С таким способом округления следует быть осторожным, он верен только для неотрицательных чисел. Если коэффициенты многочлена могут быть отрицательными, следует округлять аккуратно:

Также нужно помнить, что double имеет точность около $15$ знаков, поэтому расчитывать на точное умножение многочленов можно только если коэффициенты произведения не превосходят $10^$, а лучше и ещё меньше.

Разные мелочи

При реализации я пользовался встроенным классом complex . Чтобы его использовать, нужно подключить заголовочный файл complex .

Разумеется, можно использовать не только complex , но и complex или complex double> . Первый вариант действительно может помочь уменьшить время работы и объём потребляемой памяти.

Удивительно, но написание своего класса Complex также может уменьшить время работы.

В общем случае нам нужно 3 вызова FFT, чтобы перемножить два многочлена. Существуют методы, позволяющие проводить 2 FFT одновременно. Но если у нас есть два набора из $n$ и $m$ многочленов, и мы хотим посчитать попарные произведения, то для этого достаточно $n + m + nm$ вызовов FFT.

Операции по модулю

Часто в задачах просят посчитать что-то по модулю некоторого числа. Это же может относиться и к произведению многочленов. Конечно, в таком случае у нас все коэффициенты целые неотрицательные. Однако если модуль порядка $10^$, то коэффициенты произведения могут получиться порядка $10^$, что нельзя не то что точно сохранить в double , а даже сохранить в long long , чтобы потом взять остаток по модулю. Что же делать?

“Хороший” модуль

Почему мы вообще использовали комплексные числа, если все коэффициенты исходных многочленов были действительными? Дело в том, что нам нужно было $n$ разных корней из $1$, чего действительные числа предоставить не могут. Если бы только был какой-то другой объект с такими свойствами…

Как известно из теории чисел, ненулевые остатки по модулю простого числа $P$ образуют циклическую группу порядка $P-1$ по умножению. Пусть $P-1$ делится на достаточно большую степень двойки, то есть $P-1 = nz$ для некоторого целого $z$. Тогда в этой группе тоже есть $n$ разных корней из $1$ степени $n$, и они тоже обладают всеми необходимыми нам свойствами.

Таким образом, мы можем выполнять все операции по модулю $P$, а комплексные корни из $1$ нужно заменить на корни из $1$ по модулю.

Таким образом, если нам повезло, то всё хорошо. Но нужно понимать, что в рамках соревнований “повезло” — это значит авторы задачи сделали так, чтобы нам повезло. Поэтому если вы видите в задаче “необычный модуль”, это может быть сильным намёком на то, что в задаче требуется FFT. С другой стороны, давать такую огромную подсказку — это нежелательно для авторов. Так или иначе, коротенький список часто встречающихся “необычных модулей”, подходящих для написания FFT по модулю: $7340033 = 7 \cdot 2^ + 1$, $998244353 = 119 \cdot 2^ + 1$.

Любой модуль

Тут есть два принципиально разных подхода:

Несколько хороших модулей + Китайская теорема об остатках

Название говорит само за себя. Можно незавсимо выполнить умножение по 3 разным хорошим модулям, а потом узнать искомые коэффициенты при помощи КТО.

Разбить на многочлены с меньшими коэффициентами

Выберем $Q \approx 1000$. Запишем

Коэффициенты многочленов $A_(x)B_(x)$ будут порядка $nQ^2 \le 10^$.

Применение

Понятно, что FFT используется для умножения многочленов, так что это будет скорее описание применений умножения многочленов.

Умножение длинных чисел

Длинные числа умножаются как многочлены (коэффициенты — это цифры), только потом нужно сделать переносы. Однако лучше разбивать цифры на группы по три, то есть рассматривать число в системе счисления по основанию $1000$, — таким образом степени многочленов уменьшаются в три раза.

Вычисление скалярных произведений

Есть массив $A$ длины $n$ и массив $B$ длины $m$ ($m \le n$). Нужно посчитать скалярные произведения $B$ со всеми подотрезками $A$ длины $m$.

Перевернём $B$ и умножим с $A$ как многочлены. Несложно понять, что после переворота скалярное произведение превращается в свёртку. Коэффициенты с $(m-1)$-го по $(n-1)$-й — это ответы.

Нечёткий поиск

Есть строчки $S$ и $P$ над алфавитом $\Sigma$. Для каждой подстроки $S$ длины $ | P | $ найти расстояние Хэмминга (количество позиций, в которых строки отличаются).

Будем считать не расстояние Хэмминга, а $ | P | $ - расстояние, то есть количество совпадающих позиций. Переберём символ из $\Sigma$. В каждой строчке заменим все вхождения данного символа на $1$, а всех остальных — на $0$. Мы свели задачу к предыдущей.

Сложность — $O( | \Sigma | n \log n)$.

Используя эту задачу, можно решать задачи вида “Найдите все вхождения с не более чем $k$ ошибками”.

Все суммы

Есть два множества чисел $A$ и $B$, все числа целые от $0$ до $n$. Выдать все числа из $\lbrace x = a + b \mid a \in A, b \in B \rbrace$.

Динамическое программирование

Иногда в ДП переход выглядит как свёртка, и его можно ускорить при помощи FFT.

Читайте также: