Введение

Стандартный способ умножения двух чисел (или полиномов) требует времени O(n2). Ниже будет описано, как уменьшить время умножения до O(nlogn) при помощи быстрого преобразования Фурье (Fourier Jean Baptiste Joseph).

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

Представление полиномов

Для начала определим ключевое для дальнейших рассуждений понятие полинома. Полином A(x) от переменной x над полем F имеет вид:

A(x) = a0x0 + a1x1 + … + an1xn1

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

A(x) = a0x0 + a1x1 + … + anxn,

а

B(x) = b0x0 + b1x1 + … + bmxm,

тогда

C(x) = A(x) + B(x) = (a0+b0)x0 + (a1+b1)x1 + … + (as+bs)xs,

где s = max(n, m).

Результатом умножения двух полиномов будет их произведение — C(x) = A(x)B(x). Умножение выполняется так: каждое слагаемое многочлена A(x) умножается на каждое слагаемое многочлена B(x), после чего выполняется приведение подобных слагаемых с одинаковыми степенями.

Утверждение

Степень многочлена-произведения равна сумме степеней сомножителей:

degree(C) = degree(A) + degree(B)

Задание многочлена вектором коэффициентов

Пусть дан полином

A(x) = a0x0 + a1x1 + … + anxn

тогда вектором коэффициентов будет являться a = (a0, a1, …, an).

Задание многочлена набором значений

Фиксируем n различных точек x0, x1, …, xn1. Многочлен A(x) степени ниже n однозначно определяется своими значениями в этих точках, то есть набором из n пар аргумент-значение.

{(x0, y0), (x1, y1), …, (xn1, yn1)},

где yk = A(xk), для k = 0, 1, …, n1.

Представление многочлена с помощью значений в заданных точках удобно для многих операций: в частности, для умножения достаточно умножить значения многочленов в каждой из точек (значения многочленов должны быть заданы в одних и тех же точках). Однако, надо иметь в виду, что при умножении степени многочленов складываются, и произведение двух многочленов степени меньше n может иметь степень больше n. При этом она заведомо меньше 2n, так что для восстановления произведения достаточно 2n точек. Таким образом, умножая два многочлена степени меньше n, полезно с самого начала иметь значения многочленов не в n, а в 2n точках. Затем эти значения можно перемножить за время O(n) и получить представление произведения в виде набора пар аргумент-значение.

Теорема (однозначность интерполяции)

Для любого множества {(x0, y0), (x1,y1), …, (xn1,yn1)} пар аргумент-значение (все xi различны) существует единственный многочлен A(x) степени ниже n, для которого yk = A(xk), для k = 0, 1, …, n1.

Доказательство. Доказать данную теорему можно, представив уравнение yk = A(xk) в матричном виде.

Затем нужно проверить, что определитель данной матрицы, содержащей все корни полинома, ненулевой при условии, что все xi различны. Из этого следует, что данная матрица невырождена. А это значит, что матрица обратима. Иначе говоря, мы можем по столбцу значений найти столбец коэффициентов.

Примечание

Можно не доказывать данную теорему, а просто сослаться на истинность интерполяционной формулы Лагранжа.

Умножение полиномов

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

При этом можно использовать любой набор из n точек, но выбрав их удобным образом, можно сократить время преобразования в ту и другую сторону до O(nlogn). Как мы увидим позднее, удобно взять в качестве точек комплексные корни из единицы, в этом случае оба перехода сведутся к так называемому дискретному преобразованию Фурье и обратному к нему преобразованию, которые выполняются за O(nlogn) операций.

Представим схему умножения полиномов, заданных вектором коэффициентов.

Здесь символы i2n обозначают комлексные корни из единицы степени 2n, возведенные в степень i. При умножении двух многочленов степени меньше n получается многочлен степени меньше 2n, поэтому для начала мы дополняем многочлены-сомножители нулевыми коэффициентами старших степеней. После этого мы имеем дело с многочленами степени меньше 2n, и потому используем комплексные корни степени 2n из единицы, обозначаемые i2n при i = 0, 1, …, 2n1.

Алгоритм быстрого умножения двух полиномов можно записать так:

1. Удвоение количество коэффициентов. Дополняем полиномы A(x) и B(x) нулевыми старшими коэффициентами так, чтобы количество элементов в каждом полиноме было 2n элементов.
2. Вычисление значений. При помощи быстрого преобразования Фурье вычислить значения многочленов A(x) и B(x) в точках, являющихся корнями степени 2n из единицы.
3. Поточеченое умножение. Поточечно умножить полученные значения многочленов A(x) и B(x) друг на друга. В результате получаются значения многочлена C(x) = A(x)B(x) в корнях степени 2n из единицы.
4. Интерполяция. Получить коэффициенты многочлена C(x) при помощи обратного преобразования Фурье.

Комплексные корни из единицы

Комплексным корнем степени n из единицы называют такое комплексное число , что n = 1. Очевидно, что имеется ровно n комплексных корней для данного уравнения. Эти корни равномерно распределены на окружности единичного радиуса с центром в нуле, как показано ниже на иллюстрации.

Решения уравнения n = 1 имеют вид exp(2ikn), для k = 0, 1, …, n1. Значение = exp(2in) называют главным значеним корня степени n из единицы.

Введем несколько вспомогательных утверждений.

Лемма 1 (о сокращении)

Для любых целых k 0, n 0, d > 0 верно

dkdn = kn

Доказательство.

dkdn = exp(2i (dn))dk = exp(2i n)k = kn

Следствие

Для любого четного m 0

m2m = 12 = 1

Доказательство. Вытекает из леммы о сокращении при d = m, k = 1, n = 2

Лемма 2 (о делении пополам)

Если n > 0 четно, то, возведя в квадрат все n комплексных корней степени n из единицы, мы получим все n 2 комплексных корней степени n 2 из единицы (каждый — по два раза).

Доказательство.

kn = knn2n = k+n2n

2kn = kn2

Лемма 3 (о сложении)

Для любого целого n 0 и неотрицательного целого k, не кратного n, выполнено равенство

(kn)0+(kn)1+…+(kn)n-1 = 0

Доказательство.

Используем для доказательства формулу для геометрической прогрессии.

(kn)0+(kn)1+…+ (kn)n-1 = ((kn)n-1)/(kn-1) = ((nn)k-1)/(kn-1) = (1k-1)/(kn-1) =0

Знаменатель не обращается в нуль, так как по условию k не кратно n.

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

Пусть n является степенью двойки и пусть задан вектор коэффициентов a = (a0, a1, …, an1). Тогда

yk = A(kn) = a0+ a1(kn)1+…+ an1(kn)n1

и вектор y = (y0, y1, …, yn1) называется дискретным преобразованием Фурье.

y DFTn(a)

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

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

Выделим в многочлене A(x) отдельно члены четных и нечетных степеней.

A(x) = A[0](x2) + xA[1](x2)

A[0](x) = a0x0 + a2x1 + … + an2xn21

A[1](x) = a1x0 + a3x1 + … + an1xn21

Тем самым задача вычисления A(x) в точках 0n, 1n, …, n1n сводится к вычислению значений многочленов A[0] и A[1] степени меньше n2 в точках (0n)2, (1n)2, …, (n1n)2 и последующему комбинированию результатов.

Примечание

Лемма о делении пополам гарантирует, что список (0n)2, (1n)2, …, (n1n)2 содержит всего n 2 различных значений, а именно, комплексных корней степени n / 2 из единицы (каждый входит дважды).

Таким образом, нам нужно вычислить значения многочленов A[0] и A[1] (степени меньше n 2) в n 2 комплексных корнях степени n 2 из единицы. Эти подзадачи имеют тот же вид, что исходная, но вдвое меньший размер. Рассуждая дальше подобным образом, мы придем к рекурсивному алгоритму вычисления преобразования Фурье вектора a = (a0, a1, …, an1), (где n — степень двойки).

Реализации быстрого преобразования Фурье

Рекурсивная реализация

Рассмотрим рекурсивную реализацию быстрого преобразования Фурье.

Алгоритм

array Recursive_FFT(array a) { 1 n = length(a) 2 if (n == 1) return a 3 n = exp(2i / n) 4 = 1 5 y[0] = Recursive_FFT(a0, a2, …, an2) 6 y[1] = Recursive_FFT(a1, a3, …, an1) 7 for (k = 0; k < n / 2; k++) { 8 yk = y[0]k + y[1]k 9 yk+n/2 = y[0]k y[1]k 10 = n 11 } 12 return y }

Процедура Recursive_FFT работает следующим образом.

В строке 1 переменной n присваивается количество элементов в векторе a.

Строка 2 образует «базис рекурсии». Очевидно, что для вектора единичной длины дискретным преобразованием Фурье является сам вектор, так как степень корня нулевая.

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

В строке номер 4 инициализируется единицей переменная, в которой будет храниться текущее значение корня.

В пятой и шестой строках вычисляется дискретное преобразование Фурье для векторов, состоящих из членов только при четных или при нечетных степенях.

В цикле, состоящем из строк 7, 8, 9, 10, собираются вместе результаты рекурсивных вычислений.

Восьмая строка верна на основании рассуждений из параграфа «Быстрое преобразование Фурье».

Девятая же строка следует из того, что восьмая строка верна и из леммы о делении пополам.

Десятая и седьмая строки гарантируют нам то, что будут перебраны все комплексные корни степени n из единицы. Обоснование приведено в параграфе «Комплексные корни из единицы»

Трудоемкость данного алгоритма O(nlogn). Высота дерева рекурсивных вызовов равна logn, а трудоемкость операций, выполняемых на каждом уровне O(n).

Ниже приведено дерево входных векторов рекурсивных вызовов для n = 8.

Итеративная реализация

Рекурсию можно заменить вычислением снизу вверх, расположив элементы исходного вектора в том порядке, в котором они появляются в листьях дерева. Этот порядок можно назвать «перевернутым двоичным». (В нашем случае числа 0, 4, 2, 6, 1, 5, 3, 7 соответствуют обратной двоичной записи номеров позиций, на которых они стоят: 000, 100, 010, 110, 001, 101, 011, 111). Эта особенность связана исключительно со способом построения дерева.

Интерполяция по значениям в корнях из единицы

Операция интерполяции многочлена по заданным значениям является обратной к дискретному преобразованию Фурье (y = DFTn(a)), то есть

a = DFT1n(y)

Исходя из этого можно доказать, что для вычисления обратного преобразования Фурье можно применить тот же алгоритм быстрого преобразования Фурье, поменяв в нем местами нужные элементы. Следовательно, DFT1n также можно вычислить за время O(nlogn).

Таким образом, мы научились переходить от коэффициентов многочлена к набору его значений и обратно за время O(nlogn).

Заключение

Несмотря на то, что данный метод умножения полиномов и длинных чисел не требует больших познаний в области математики, первое упоминание об употреблении преобразования Фурье для умножения длинных чисел было лишь в 1968 г. Этод метод был предложен Штрассеном (Strassen Volker). Позже вместе с Шёнхаге (Schnhage Arnold) он уточнил метод и опубликовал модифицированные алгоритмы в Computing 7 (1971), 281-292. Похожие идеи, но для произвольных целых чисел, были опубликованы Поллардом (Pollard John Michael) в Math. Comp. 25 (1971), 365-374.

Дискретное преобразование Фурье использует вычисления с комплексными числами, что может привести к потере точности из-за ошибок округления. Это особенно нежелательно, если исходными данными и результатами являются целые числа. Для таких случаев можно использовать вариант преобразования Фурье в кольце вычетов (в котором вычисления проводятся точно). См. [2].

В то же время, погрешность результата относительно невелика. Подробно ознакомиться с материалом по этой теме можно в [1].

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

Источники

1. Кнут Д. Искусство программирования, том 2. Получисленные алгоритмы, 3-е изд. — М.: Издательский дом «Вильямс», 2001.
2. Кормен Т., Лейзерсон Ч., Ривест Р. Алгоритмы: построение и анализ. — М.:МЦНМО, 2001.
3. Математика. Большой энциклопедический словарь. — М.: Большая Российская Энциклопедия, 2000.
4. Fidonet, конференция SU.HARDW.

Александр Горюнов

Кирилл Сомик / 2008-02-28 17:53:59

Спасибо за статью! Лаконично и предельно ясно изложен важный и сложный метод ВМ. В этой связи имеются предложения по сотрудничеству с автором - Александром Горюновым. Александр! Прошу Вас выйти на связь.

Вряд ли автор статьи заглядывает на эту страницу. А его нынешний электронный адрес мы, к сожалению, не знаем.

Антон Москвичёв / 2009-05-29 10:34:28

Есть специальные методы быстрого умножения полиномов. Многие из них описаны в книге Нуссбаумера "Быстрое преобразование Фурье и алгоритмы вычисления сверток". То, что описано здесь, подходит толькодля полиномов с небольшими коффициентами и не слишком больших длин. Если у полинома коэффициенты большие (от нескольких десятков значащих цифр) или его длина большая (от нескольких тысяч элементов), то вышеописанный метод становится очень медленным. Специальные полиномиальные преобразования же вообще не используют умножений (точнее, они используются, но формально или , как говорят, тривиально), только сложения-вычитания и деление пополам. Если в связи с этим имеются предложения по сотрудничеству - обращайтесь.