Алгоритмы быстрого умножения чисел: от столбика до Шенхаге-Штрассена03.03.2023 14:01

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

И уж конечно, никогда при написании a * b мы не задумываемся о том, как реализовано умножение чисел a и b в нашем языке. Какие вообще есть алгоритмы умножения? Это какая-то нетривиальная задача?

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

Его величество столбик
• Про О-нотацию
Разделяй и властвуй: алгоритм Карацубы
Многочлены vs преобразование Фурье: алгоритм Шенхаге-Штрассена
• Про быстрое преобразование Фурье
Модульная арифметика и второй алгоритм Шенхаге-Штрассена
• Про модульную арифметику чисел
• Про модульную арифметику многочленов

Зачем быстро умножать числа?

Программы постоянно перемножают числа. Перемножения 32-битных, 64-битных, а иногда и более длинных чисел встроены напрямую в арифметико-логические устройства микропроцессоров; генерация оптимальных цепей для кремния — отдельная инженерная наука. Но не всегда встроенных возможностей хватает.

Например, для криптографии. Криптографические алгоритмы (вроде повсеместно используемого RSA) оперируют числами длиной в тысячи бит. Если мы сводим операции над 4096-битными числами к операциям над 64-битными словами, разница в количестве операций между алгоритмами за N^2 и $N\log N$ уже составляет десять раз!

Деление тоже сводится к умножению; в некоторых процессорах даже нет инструкции для целочисленного деления.

Но сразу признаюсь: не все существующие алгоритмы быстрого умножения достаточно практичны для широкого применения — по крайней мере на сегодняшний день. Помимо практического здесь силён академический интерес. Как вообще математически устроена операция умножения? Насколько быстро её можно делать? Можем ли мы найти оптимальный алгоритм и чему научимся в процессе поиска?

Его величество столбик

Когда-то очень давно мне попалось на глаза видео с кричащим названием «How to Multiply», описывающее так называемый японский метод умножения. Что меня удивило — так это то, что метод этот подаётся как более простой и интуитивно понятный. Оказывается, если умножение в столбик делать не в столбик, а занимая половину листа бумаги, получается проще и понятнее!

Да-да, принципиально это тот же самый алгоритм умножения столбиком. Можем посмотреть на запись в столбик и сравнить её с картинкой:

92 — это две больших группы красных точек на нижне-правой диагонали. 23 — две маленьких группы на верхне-левой диагонали. Точно так же, как и в столбике, мы вынуждены перемножить попарно все разряды, потом сложить, и потом сделать перенос. Умножение в столбик — не что иное, как компактный на бумаге и относительно удобный для человеческого сознания способ проделать алгоритм, объединяющий в себе практически все придуманные до XX века способы умножения, за редким исключением вроде умножения египетских дробей.

Асимптотическая сложность умножения в столбик не зависит от того, в какой системе счисления мы производим умножение, и составляет $O(M\cdot N)$ операций, где и — количество разрядов в множителях. Каждый из разрядов одного множителя нужно умножить на каждый из разрядов второго множителя, после чего получившиеся маленьких чисел (в каждом не больше двух разрядов) сложить.

Почему асимптотика не зависит от системы счисления?

В разных системах счисления у чисел разное количество разрядов, это верно. Количество десятичных разрядов в числе x равно $\lceil\log_{10}x\rceil$ в двоичной — $\lceil\log_{2}x\rceil$ (скобки-уголки означают округление вверх). Но мы знаем, что логарифмы по разному основанию связаны друг с другом константным множителем:

$\log_{10}x=\log_{10}2\cdot\log_{2}x$

А константные множители не имеют значения при анализе асимптотики.

Обычно при анализе сложности под и подразумевается количество двоичных разрядов, а числа предполагаются примерно одинаковой длины; тогда формула сложности упрощается до O(N^2) .

Что означает запись O (N²)?

Многих программистов О-нотации учит улица — они сталкиваются с ней сразу в оценке сложности каких-то алгоритмов. Но О-формализм происходит из математического анализа, и у него есть строго определённый смысл, причём не очень хорошо согласующийся с остальными привычными для математики обозначениями.

Если у нас есть функция f(N) — в наших примерах это, как правило, будет количество операций, нужных для обработки алгоритмом входных данных в размере бит — то запись f(N) = O(N^2) означает, что существует некоторая константа C такая, что начиная с достаточно больших

$f(N) \le C\cdot N^2$

То есть функция растёт не быстрее, чем N^2.

При этом сама запись O(N^2) означает класс функций, для которых выполнено это условие. В привычкой теоретико-множественной нотации это было бы правильно записать как $f(N)\in O(N^2)$ — функция входит во множество функций, растущих не быстрее чего-то там.

При этом ничего не сказано о том, может ли функция расти медленнее, чем N^2 . Например, в случае алгоритма умножения в столбик мы могли бы без зазрения совести записать f(N) = O(N^3) — и формально всё ещё были бы правы. Действительно, операций нужно меньше, чем N^3 . Полезная ли это информация? Не очень.

Соответственно, одна и та же функция может быть «равна» разным $O(\ldots)$ ; при этом между друг другом эти $O(\ldots)$ не становятся равны:

$f(N) = O(N^2),\ f(N) = O(N^3)\ \ \not{\!\!\!\Longrightarrow}\ \ O(N^2) = O(N^3)$

Так что обращаться с $O(\ldots)$ в выражениях нужно осторожно.

Помимо $O(\ldots)$ — наиболее, пожалуй, распространённого среди программистов — есть и другие классы функций с аналогичной записью. f(N) = o(N^2) означает, что функция растёт медленее, чем N^2 ; да не просто медленнее, а настолько, что их частное становится с ростом всё ближе и ближе к нулю:

$\frac{f(N)}{N^2} \rightarrow 0 \mbox{ при } N\rightarrow\infty$

Например, в матанализе постоянно используют запись o(1) , чтобы обозначить незначительную величину, стремящееся к нулю. Единица в скобках при этом не имеет какого-то существенного значения — там могла бы быть любая константа; можете написать o(69) и формально это будет верно. Иногда бывает важно указать, с какой именно скоростью величина стремится к нулю. Тогда пишут $o(1/\sqrt{N})$ , o(1/N) и так далее.

Другая встречающаяся нотация $f(N) = \Theta(N^2)$ означает, что функция растёт точно со скоростью N^2 , то есть существуют константы C_1 и C_2 такие, что начиная с достаточно больших

$C_1\cdot N^2 \le f(N) \le C_2\cdot N^2$

Помимо простого использования со знаком равенства, как в примерах выше, классы функций можно использовать в арифметических операциях. Например, запись

означает, что функция равна N^2 плюс нечто маленькое, стремящееся к нулю при росте . В отличие от f(N) = O(N^2) или $f(N) = \Theta(N^2)$ , эта запись не допускает произвольных множителей при N^2 — ровно N^2 и всё тут.

Работает (при некоторых ограничениях) и привычная арифметика в равенствах. Так, если f(N) = O(N^3) , обе части равенства можно поделить, например, на константу:

$\frac{f(N)}{2} = \frac{O(N^3)}{2} = O(N^3)$

$O(\ldots)$ съедает любые константы, поэтому справа ничего не изменилось. А можно поделить на :

$\frac{f(N)}{N} = \frac{O(N^3)}{N} = O(N^2)$

Не-константу нельзя просто выкинуть, поэтому она ушла внутрь $O(\ldots)$ и сократилась с тем, что там было.

Разделяй и властвуй: алгоритм Карацубы

Первый шаг на пути ускорения умножения совершил в 1960-м году советский математик Анатолий Карацуба. Он заметил, что если длинные числа поделить на две части:

${\color{Red}{235739}}\,{\color{DarkBlue}{098113}}\ \times\ {\color{DarkOrange}{187129}}\,{\color{Magenta}{102983}} \\[4mm] \left({\color{Red}{a_1}}\cdot 10^6 + {\color{DarkBlue}{a_0}}\right)\ \times\ \left({\color{DarkOrange}{b_1}}\cdot 10^6 + {\color{Magenta}{b_0}}\right) \\[4mm] {\color{Red}{a_1}}\cdot{\color{DarkOrange}{b_1}}\cdot 10^{12} + ({\color{Red}{a_1}}\cdot {\color{Magenta}{b_0}} + {\color{DarkBlue}{a_0}}\cdot {\color{DarkOrange}{b_1}})\cdot 10^6 + {\color{DarkBlue}{a_0}}\cdot {\color{Magenta}{b_0}}$

То можно обойтись тремя умножениями этих более коротких частей друг на друга, а не четырьмя, как можно было бы подумать. Вместо прямого подсчёта ${\color{Red}{a_1}}\cdot {\color{Magenta}{b_0}} + {\color{DarkBlue}{a_0}}\cdot{\color{DarkOrange}{b_1}}$ , требующего двух умножений, достаточно посчитать $\left({\color{Red}{a_1}} + {\color{DarkBlue}{a_0}}\right) \times \left({\color{DarkOrange}{b_1}} + {\color{Magenta}{b_0}}\right)$ и вычесть из результата числа ${\color{Red}{a_1}}\cdot{\color{DarkOrange}{b_1}}$ и ${\color{DarkBlue}{a_0}}\cdot {\color{Magenta}{b_0}}$ , нужные нам в любом случае для получения младших и старших разрядов.

После этих расчётов остаётся лишь пробежаться по разрядом справа налево и провести суммирование:

$\begin{array}{c@{\quad=\quad}rrrr} {\color{DarkBlue}{a_0}}\cdot {\color{Magenta}{b_0}} & &&\!\!\!\!\!\!10103&\!\!\!\!\!\!9710\overset{\displaystyle\leftarrow}{79} \\ {\color{Red}{a_1}}\cdot {\color{Magenta}{b_0}} + {\color{DarkBlue}{a_0}}\cdot {\color{DarkOrange}{b_1}} & &\!\!\!\!\!\!42636&\!\!\!\!\!\!897014& \\ {\color{Red}{a_1}}\cdot{\color{DarkOrange}{b_1}} & 44113&\!\!\!\!\!\!603331&& \\ \hline a \cdot b & 44113&\!\!\!\!\!\!645967&\!\!\!\!\!\!907117&\!\!\!\!\!\!971079 \end{array}$

По сравнению с умножениями это быстрая операция, не заслуживающая особого внимания — как и два лишних сложения в скобках.

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

Поскольку мы каждый раз делим задачу на три задачи с вдвое меньшими (по количеству бит) числами, для перемножения двух чисел длины 2^k нам потребуется рекурсия глубины и суммарно 3^k умножений чисел наименьшей длины; отсюда получаем оценку сложности в $O(N^{\log_2 3}) \approx O(N^{1.58}).$

Занятный факт — аналогичные фокусы с экономией за счёт одновременных умножений используются ещё в ряде мест. Так, два комплексных числа и также можно перемножить за три вещественных умножения вместо четырёх, вычислив , и $(a+b)\cdot(c+d)$ . А алгоритм Пана для быстрого перемножения матриц основывается на том, что можно одновременно вычислить произведения двух пар матриц и за меньшее число умножений, чем по отдельности [1, 2]. Не говоря уж о том, что само по себе перемножение матриц быстрее, чем за — это быстрое одновременное умножение матрицы на векторов.

Числа vs многочлены: алгоритмы Тоома-Кука

При виде магического сокращения вычислительной сложности при разбиении множителей на две части как-то сам собой возникает вопрос:, а можно разбить множители на бóльшее число частей и получить бóльшую экономию?

Ответ — да; и подход, позволяющий это сделать, включает в себя алгоритм Карацубы как частный случай. Но для его формулировки нам придётся проделать некий фокус.

Давайте попробуем разбить те же самые числа на три части:

${\color{Red}{2357}}\,{\color{DarkBlue}{3909}}\,{\color{DarkGreen}{8113}}\ \times\ {\color{DarkOrange}{1871}}\,{\color{Magenta}{2910}}\,{\color{Golden}{2983}} \\[5mm] \left({\color{Red}{a_2}}\cdot 10^8 + {\color{DarkBlue}{a_1}}\cdot 10^4 + {\color{DarkGreen}{a_0}}\right)\ \times\ \left({\color{DarkOrange}{b_2}}\cdot 10^8 + {\color{Magenta}{b_1}}\cdot 10^4 + {\color{Golden}{b_0}}\right)$

Если мы внимательно посмотрим на коэффициенты при степенях десятки, которые возникают при перемножении этих скобок:

$\begin{array}{r} {\color{Teal}{a_0\cdot b_0}} \\ {\color{Violet}{a_1\cdot b_0 + a_0\cdot b_1}} \quad\quad \\ {\color{Orchid}{a_2\cdot b_0 + a_1\cdot b_1 + a_0\cdot b_2}} \quad\quad\quad\quad \\ {\color{Gray}{a_2\cdot b_1 + a_1\cdot b_2}}\quad\quad\quad\quad\quad\quad \\ {\color{DarkRed}{a_2\cdot b_2}}\quad\quad\quad\quad\quad\quad\quad\quad \\ \hline a\cdot b \end{array}$

То (при наличии опыта в алгебре) заметим, что где-то мы такое уже видели. При перемножении многочленов!

Если взять части наших чисел a и b и объявить их коэффициентами многочленов f(x) и g(x) при соответствующих степенях, числа в таблице выше будут не чем иным, как коэффициентами многочлена $f(x) \cdot g(x)$ :

$\begin{array}{l} f(x) \cdot g(x) = \\[2mm] \quad\quad = (a_0 + a_1x + \ldots)\ \times\ (b_0 + b_1x + \ldots) = \quad\quad\quad\quad \\[2mm] \qquad\qquad\begin{array}{r@{\ }l} = & ({\color{Teal}{a_0\cdot b_0}}) \\ + & ({\color{Violet}{a_1\cdot b_0 + a_0\cdot b_1}})\cdot x \\ + & ({\color{Orchid}{a_2\cdot b_0 + a_1\cdot b_1 + a_0\cdot b_2}})\cdot x^2 \\ + & ({\color{Gray}{a_3\cdot b_0 + a_2\cdot b_1 + a_1\cdot b_2 + a_0\cdot b_3}})\cdot x^3 \\ + & \ldots \end{array} \end{array}$

Сами же числа получаются из многочленов путём вычисления значения в некоторой точке. Какой именно — зависит от системы счисления и размера частей:

$a = {\color{Red}{a_2}}\cdot 10^8 + {\color{DarkBlue}{a_1}}\cdot 10^4 + {\color{DarkGreen}{a_0}} = f(10^4)$

Хорошо, свели одну задачу к другой. В чём преимущество многочленов? Помимо того, что это некоторые формальные конструкции с параметрами, это также функции; чтобы вычислить $f(x) \cdot g(x)$ в произвольной точке, не нужно перемножать многочлены — достаточно вычислить оба значения в этой точке и перемножить их. Если бы я писал продакшн-код, в котором по какой-то причине нужно перемножать многочлены, я бы и вовсе сделал перемножение ленивым.

Но нам не нужно вычислять значение в точке. Нам нужны коэффициенты. А коэффициенты многочлена можно восстановить по значениям в точках, решив систему линейных уравнений:

$\left\{\begin{matrix} c0 + c_1\cdot x_1 + \ldots = f(x_1)\cdot g(x_1) \\ c0 + c_1\cdot x_2 + \ldots = f(x_2)\cdot g(x_2) \\ \vdots \end{matrix}\right.$

Здесь x_i — это точки, в которых нам известны значения многочлена, в правой части — сами эти значения, а c_i — неизвестные нам коэффициенты многочлена. Количество неизвестных соответствует степени многочлена + 1; чтобы решение было единственным, нужно, соответственно, знать значения в таком же количестве точек. В нашем примере это пять точек, потому что у многочлена-произведения степень 4:

$f(x) \cdot g(x) = {\color{DarkRed}{a_2\cdot b_2}}\,\cdot x^4 + \ldots$

Если переписать эту систему уравнений в матричном виде, получим

$\begin{bmatrix} 1 & x_1 & x_1^2 & x_1^3 & \cdots \\[1mm] 1 & x_2 & x_2^2 & x_2^3 & \cdots \\[1mm] 1 & x_3 & x_3^2 & x_3^3 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix} \times \begin{bmatrix} c_0 \\[1mm] c_1 \\[1mm] c_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} f(x_1)\cdot g(x_1) \\[1mm] f(x_2)\cdot g(x_2) \\[1mm] f(x_3)\cdot g(x_3) \\ \vdots \end{bmatrix}$

Соответственно, для создания алгоритма быстрого умножения чисел нам достаточно:

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

В ходе алгоритма нам нужно будет два раза умножить на прямую матрицу (по разу для каждого из чисел, которые мы хотим перемножить), перемножить результаты алгоритмом умножения чисел меньшего размера, и один раз умножить получившийся вектор на обратную матрицу.

Давайте посмотрим на матрицу для точек 0, 1, -1, 2, -2 :

from sympy import Rational
from sympy.matrices import Matrix

# Точки, в которых будем вычислять значения многочленов:
x = [0, 1, -1, 2, -2]

# Матрица Вандермонда, построенная по этим точкам:
v = Matrix([
    [
        Rational(x_i ** j )
        for j, _ in enumerate(x)
    ] for x_i in x
])

$\left[\begin{matrix}1 & 0 & 0 & {\color{Gray}{0}} & {\color{Gray}{0}}\\1 & 1 & 1 & {\color{Gray}{1}} & {\color{Gray}{1}}\\1 & -1 & 1 & {\color{Gray}{-1}} & {\color{Gray}{1}}\\1 & 2 & 4 & {\color{Gray}{8}} & {\color{Gray}{16}}\\1 & -2 & 4 & {\color{Gray}{-8}} & {\color{Gray}{16}}\end{matrix}\right]$

(Серым я отметил числа, не имеющие значения, поскольку соответствующие коэффициенты в многочленах-множителях заведомо равны нулю.)

И на её обратную:

$\frac{1}{24}\times\left[\begin{matrix}24 & 0 & 0 & 0 & 0\\0 & 16 & -16 & -2 & 2\\-30 & 16 & 16 & -1 & -1\\0 & -4 & 4 & 2 & -2\\6 & -4 & -4 & 1 & 1\end{matrix}\right]$

За исключением некоторой проблемы с делением на 3, подавляющая часть вычислений здесь представляется битовыми сдвигами, сложениями и вычитаниями «коротких» чисел, а это операции «простые» — имеющие линейную сложность.

Как проанализировать итоговую сложность получившегося алгоритма? Давайте обозначим количество операций через C(N) , где — количество бит в наших множителях, и выведем рекуррентную формулу:

$C(N) = 5\cdot C\left(\frac{N}{3}\right) + \mathrm{const}\cdot \frac{N}{3}$

Что здесь написано: умножение двух чисел длиной бит мы умеем сводить к 5 умножениям чисел длиной N/3 бит и ещё какому-то количеству простых операций над числами длиной N/3 бит. Применяем основную теорему о рекуррентных соотношениях и получаем итоговую сложность

$C(N) = O(N^{\log_3 5}) \approx O(N^{1.46})$

Таким образом мы ускорили умножение с $O(N^{1.58})$ до $O(N^{1.46})$ .

Алгоритмы, получаемые таким образом, называются алгоритмами Тоома-Кука. Они активно используются на практике; в одной только библиотеке GMP поддержано 13 разных вариантов разбиения чисел.

Закономерный следующий вопрос: что дальше? Если число бьётся на частей, сложность получается

$C(N) = O\left(N^{{\color{DarkBlue}{\log_M (2M+1)}}}\right)$

Показатель степени при большом можно упростить:

${\color{DarkBlue}{\log_M (2M+1)}} \approx \log_M (2M) = \frac{\log_2{2} + \log_2 M}{\log_2 M} = 1 + \frac{1}{\log_2 M} \rightarrow 1$

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

Во-первых, это непрактично. Логарифм возрастает крайне медленно; асимптотика, соответственно, тоже будет падать медленно, а константа при этом будет быстро расти, потому что числа в матрицах будут всё больше и разнообразнее.

Во-вторых, это неинтересно. Да-да! Дальше мы будем обсуждать алгоритмы, сложность которых лучше $O(N^{1+\varepsilon})$ , каким маленьким бы ни было $\varepsilon$ .

Занятный факт — аналогичной методу Тоома-Кука техникой строятся быстрые алгоритмы перемножения матриц. Однако процедура перемножения матриц по своей математической природе оказалась намного сложнее перемножения многочленов; поэтому до сих пор нет уверенности, что можно перемножать матрицы за $O(N^{2+\varepsilon})$ для сколь угодно малого $\varepsilon$ . Лучшая асимптотика на сегодняшний день составляет примерно $O(N^{2.37})$ .

Многочлены vs преобразование Фурье: алгоритм Шенхаге-Штрассена

Раз рекуррентное деление на несколько частей нам не интересно, давайте попробуем сделать радикальный шаг. Что, если делить числа не на фиксированное количество частей, а на части фиксированной длины?

Например, на десятичные разряды. В таком случае нам нужно будет 2N-1 точек для -разрядного числа; в остальном подход как будто бы тот же самый, что в предыдущем разделе. Что ж, давайте возьмём x_i = i и проведём численный эксперимент, чтобы убедиться, что всё работает!

import math
import numpy as np

# Числа, которые мы изначально собирались перемножить:
a = 235739098113
b = 187129102983
m = len(str(a))

# Коэффициенты соответствующих многочленов в порядке 
# возрастания степени. Поскольку результирующий многочлен 
# будет степени 2*m-1, добиваем нулями до нужной ширины:
a_coefs = [int(a_i) for a_i in str(a)[::-1]] + [0] * (m - 1)
b_coefs = [int(b_i) for b_i in str(b)[::-1]] + [0] * (m - 1)
n = len(a_coefs)

# Точки, в которых будем вычислять значения многочленов:
x = np.arange(n)

# Матрица Вандермонда, построенная по этим точкам:
v = np.vander(x, increasing=True)

# Вычисление значения многочлена в точке — не что иное, как
# умножение этой матрицы на вектор коэффициентов:
a_y = v @ a_coefs
b_y = v @ b_coefs

# Поточечно перемножаем значения, чтобы получить 
# значения многочлена-произведения:
c_y = a_y * b_y

# Восстанавливаем коэффициенты многочлена-произведения:
c_coefs = np.linalg.solve(v, c_y)

# Для сверки считаем коэффициенты "в лоб":
actual_c_coefs = [
    sum(a_coefs[j] * b_coefs[i-j] for j in range(i+1)) 
    for i in range(n)
]

# Считаем длину вектора-разности и делим на длину настоящего,
# чтобы получить относительную погрешность:
print(np.linalg.norm(c_coefs - actual_c_coefs) /
    np.linalg.norm(actual_c_coefs))

Запускаем и…

2.777342817120168e+17

Хо-хо, да это не просто мимо, это фантастически мимо! Но почему?

Оказывается [3], у матриц Вандермонда есть неприятное свойство — решать системы с ними в подавляющим большинстве случаев очень плохая идея, потому что погрешность результата растёт экспоненциально с размером матрицы. Можно улучшить ситуацию, взяв вместо x_i = i комплексные точки на единичной окружности:

x = np.exp(1j * np.random.uniform(0, math.pi, size=n))

...

0.00015636299542432733

Но лучшим выбором оказывается взять точки, равномерно распределённые на единичной окружности и являющиеся корнями из единицы степени .

x = np.exp(-1j * np.linspace(0, 2*math.pi, n, endpoint=False))

Что ещё за корни из единицы?

Запускаем новый вариант кода и получаем

1.5964929527133826e-15

Отлично! Мы дошли до предела машинной точности.

Но самое изумительное — при таком выборе иксов матрица представляет собой не что иное, как матрицу дискретного преобразования Фурье! А значит, мы можем умножать векторы на эту матрицу и решать систему уравнений с ней алгоритмом быстрого преобразования Фурье за $O(N \log N)$ операций сложения и умножения — вместо наивного алгоритма за O(N^2) .

Что такое дискретное преобразование Фурье?

Преобразование Фурье досталось нам из высшей математики и волновой физики. В физике оно известно как способ разложить сигнал (например, аудио) какой-то произвольной формы

на сумму элементарных сигналов — гармоник с длиной волны, кратной ширине отрезка:

Одна гармоника соответствует постоянному звуку некой фиксированной частоты и громкости; сумма гармоник может дать любую мелодию. Сколь угодно сложный сигнал можно разложить на гармоники, если взять их достаточно много. Совокупность гармоник (частоты + амплитуды) называется спектром сигнала.

Преобразование Фурье — построение в первую очередь теоретическое; Жозеф Фурье открыл его задолго до появления компьютеров. На практике мы не можем оперировать сигналами как функциями; у нас есть только дискретизации, замеры значения сигнала с каким-то шагом по времени.

Согласно теореме Котельникова-Найквиста-Шеннона при такой ограниченной информации мы всё ещё можем восстановить гармоники (а с ними — и весь сигнал), если предположим, что сигнал был достаточно простым (в нём не было каких-то сверхвысоких частот). Этим и занимается дискретное преобразование Фурье — фактически это матрица, которую нужно умножить на вектор со значениями дискретизированного сигнала, чтобы получить вектор с амплитудами, соответствующими разным частотам гармоник.

Так, дискретизация на картинке выше недостаточно частая, чтобы восстановить самую высокую частоту:

А более частая, удовлетворяющая условию теоремы — достаточна:

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

В целом анализ Фурье — очень интересный раздел математики, связывающий воедино множество ранее разобщённых концепций. Его подробное описание, конечно, уходит совсем далеко за рамки статьи; можно начать со статей на Хабре «Простыми словами о преобразовании Фурье» и «Преобразование Фурье в действии: точное определение частоты сигнала и выделение нот».

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

Дискретное преобразование Фурье применяется повсеместно для обработки сигналов (например, шумоподавления), часто реализуется аппаратно и встречается с совершенно неожиданных устройствах (например, МРТ-сканерах [4]). Это создаёт потребность в создании настолько быстрого алгоритма, насколько возможно.

Умножение на матрицу «в лоб» стоит O(N^2) — это считается довольно медленным в практических алгоритмах. Быстрое преобразование Фурье — название семейства алгоритмов, достигающих асимптотики $O(N\log N)$ за счёт сведения задачи к нескольким задачам меньшего размера с линейной стоимостью одного шага рекурсии. Самый простой вариант — алгоритм для N=2^k , сводящий задачу к двум задачам размера $2^{k-1}$ ; рекурсия глубины приводит к общей сложности $O(k\cdot 2^k)$ .

Про него в интернете, конечно, написано много материалов разной степени читабельности. Но давайте я попробую очень быстро на пальцах показать, за счёт чего вычисление можно ускорить до $O(N\log N)$ .

Дискретное преобразование Фурье и обратное дискретное преобразование Фурье — восстановление сигнала по спектру — очень похожи математически, и алгоритм почти идентичный; объяснить его будет проще на примере обратного преобразования.

Обратное преобразование заключается в том, что нужно взять гармоник, у которых мы знаем частоту и амплитуду, посчитать их значения в точках отрезка $[0, 2\pi]$ и сложить, получив вектор из значений сигнала-суммы гармоник. На рисунке пример для N=8 :

Поделив отрезок пополам, можно заметить, что гармоники делятся на два типа:

В верхней части — кривые, половинки которых одинаковы в первой и второй половине отрезка; в нижней — выглядящие вертикально отражёнными, то есть одна половина получается из другой умножением на . А значит, можно сэкономить, вычисляя значения всех гармоник только на половине отрезка!

При этом задача для произвольной гармоники на половине отрезка соответствует задаче для гармоники с вдвое меньшей частотой и N/2 точек:

Получается, что для решения исходной задачи восстановления сигнала в точках по данным гармоникам нам достаточно:

просуммировать гармоник первого типа (взяв вдвое меньшую частоту) в точках;
просуммировать гармоник второго типа (также взяв вдвое меньшую частоту) в точках;