Быстрое умножение многочленов при помощи преобразования Фурье — это просто

Если Вы знакомы с комплексными числами, то можете пропустить этот пункт, в противном случае, вот краткое определение: x=a+ib, где i 2 =-1 Здесь a называется вещественной (Real) частью, а b — мнимой (Imaginary). В этих числах, как нетрудно заметить, можно извлекать корень из отрицательных (да и вообще любых) чисел — это очень удобно при работе с многочленам — как следует из основной теоремы алгебры, у каждого многочлена степени n имеется ровно n комплексных корней (с учётом кратности). Также их очень удобно представлять в виде точек на плоскости:

Еще одним замечательным свойством комплексных чисел является то, что их можно представить в виде x=(cosα+isinα)r, где α — полярный угол «числа» (называется аргументом), а r — расстояние от нуля до него (модуль). А при умножении двух чисел: a=(cosα+i⋅sinα)ra b=(cosβ+i⋅sinβ)rb ab=(cosα+i⋅sinα)(cosβ+i⋅sinβ)rarb ab=(cosα⋅cosβ-sinα⋅sinβ+i(sinα⋅cosβ+cosβ⋅sinα))rarb ab=(cos(α+β)+i⋅sin(α+β))rarb Их модули перемножаются, а аргументы складываются.

Теперь давайте поймём, как выглядят комплексные корни n-ой степени из 1. Пусть x n =1, тогда его модуль, очевидно, равен единице, а n⋅argx=2πk, где k — целое. Это обозначает, что после n умножений числа на самого себя (т.е. возведения в n-ю степень) его аргумент станет «кратен» 2π (360 градусам). Вспомним формулу числа, если известен аргумент и модуль, получаем: α=2π⋅x/n, где 0x ωi=cosα+i⋅sinα Т.е. если порисовать, то мы получим просто точки на окружности через равные промежутки:

Прошу заметить три вещи, которыми мы будем активно пользоваться (без них ничего не получится): ωa⋅ωb=ω(a+b)modn ω0+ω1+ω2+. +ωn-1=0 ω0 n/2 =ω2 n/2 =ω4 n/2 =. =1 (при чётном n) Из-за этих свойств именно в этих точках мы и будем считать значение многочлена. Разумеется, результаты необязательно будут вещественными, поэтому в программе потребуется работать с комплексными числами.

Как работает

Напишем что-нибудь

Для того, чтобы сделать процедуру нерекурсивной, придётся подумать. Легче всего, как мне кажется, провести аналогию с MergeSort (сортировка слиянием) и нарисовать картинку, на которой показаны все рекурсивные вызовы:

1. Пробежимся циклом от 0 до n-1
2. Будем хранить и динамически пересчитывать номер старшего единичного бита числа. Он меняется, только когда текущее число — степень двойки: увеличивается на 1.
3. Когда мы знаем старший бит числа, перевернуть всё число не составляет труда: «отрезаем» старший бит (XOR), переворачиваем остаток (уже посчитанное значение) и добавляем «отрезанную» единицу

Аргументами процедуры для слияния двух блоков будут два vector'а (естесственно, по ссылке, исходный и результат), номер стартового элемента первого блока (второй идёт сразу после) и длина блоков. Можно было бы конечно сделать и iterator'ами — для большей STL'ности, но мы ведь всё равно будем переносить эту процедуру внутрь основной для краткости.

И основная процедура преобразования:

На данный момент весь код преобразования занимает порядка 55 строк. Не сотню, но это достаточно много — можно короче. Дляначала избавимся от кучи лишних переменных и операций в fft_merge:

Теперь можно переместить цикл из fft_merge в основную процедуру (также можно убрать p2, поскольку p2=p1+len — у меня это также дало небольшой выигрыш по времени. Что любопытно, если убрать p1=pdest, то у меня лично выигрыш по времени убивается):

В результате теперь код выглядит так:

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

Пусть мы применяем обратное преобразование к многочлену P(x) с коэффициентами vi (исходный многочлен имел коэффициенты ai): vi=a0+ωia1+ω2ia2+ω3ia+. Посмотрим на результат преобразования: bi=v0+ωiv1+ω2iv2+ω3iv3+. Подставим значения vj (помним, что ωaωb=ωa+bmodn:

Теперь давайте докажем один замечательный факт: при x≠0, ω0+ωx+ω2x+. +ω(n-1)x=0. Доказывается аналогично тому, что сумма корней — ноль: обозначим за φ сумму, домножим обе части на ωx и посмотрим, что получилось. Теперь применим этот факт к вычислению значения bi. Заметим, что все строки, кроме одной, в которой содержится an-i, обнулятся.

bi=an-i⋅n

Что и требовалось доказать.

Применение

Легко заметить, что время работы этой программы — O(n⋅log2n) и самые трудоёмкие операции — преобразования Фурье. Также можно заметить, что если нам требуется вычислить более сложное выражение с двумя многочленами, то по-прежнему можно выполнять лишь три приобразования — сложение и вычитание также будут работать за линейное время. К сожалению, с делением не всё так просто, поскольку многочлен может случайно принять значение 0 в какой-нибудь из точек. UPD2: не забудьте, что степень произведения двух многочленов степени n будет равна 2n, поэтому при вводе следует добавить «лишние» нулевые старшие коэффициенты. Если представить число в десятичной (или более) системе счисления, как многочлен с коэффициентами — цифрами, то умножение длинных чисел также можно выполнять очень быстро. И, напоследок, задача, которую я разберу в следующем посте: у вас есть две строки одинаковой длины порядка 10 5 из букв A, T, G, C. Требуется найти такой циклический сдвиг одной из строк, чтобы совпало максимальное количество символов. Очевидно наивное решение за O(n 2 ), но есть решение при помощи БПФ. Удачи!