patnashev Aug 22 at 08:06

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

10 min

11K

В одной из предыдущих статей я рассказал о математических алгоритмах, позволяющих проверить простоту очень большого числа. Но в основе всех тех алгоритмов лежит одна базовая операция — перемножение двух больших чисел. Именно операции длинного умножения занимают 99,9% времени выполнения любого теста простоты. Как же умножение реализуется на практике? Говорят, что при помощи быстрого преобразования Фурье. Но беглое прочтение Википедии вызывает недоумение. Какое отношение преобразование Фурье имеет к умножению целых чисел? Давайте разбираться.

Преобразование Фурье широко используется нами в повседневной жизни. Его можно найти в таких технологиях как JPEG, MP3, Wi-Fi, связь 4G. Это один из базовых алгоритмов цифровой обработки сигналов. Любой сигнал, записанный в виде отдельных точек, будь то двумерная картинка или одномерный звук, можно разложить на составляющие при помощи дискретного преобразования Фурье.

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

В случае цифровой обработки сигналов, выявление синусоид является также одной из главных задач. Самый простой способ сделать это — взять синусоиду известной частоты и каждую её точку умножить на соответствующую точку сигнала. Если в сигнале точно такая же синусоида, то мы получим квадрат синусоиды (случай 1 на иллюстрации). Так как квадрат всегда положителен, то, просуммировав все точки результата, мы получим большое положительное число. Это явное указание на то, что в исходном сигнале такая синусоида есть. Если же частота синусоиды в сигнале отличается от заданной, то в произведении будут как положительные, так и отрицательные числа. Их сумма будет стремиться к нулю (случай 2). Ноль означает, что в сигнале нет синусоиды частоты , и нам надо попробовать другие частоты. Такой подход не даёт ответа, какая именно частота у исходного сигнала, но позволяет найти её путём перебора.

Однако, кроме частоты у синусоиды есть ещё такая характеристика, как фаза, или сдвиг относительно нулевой точки. В случае сдвига на половину периода, мы получаем зеркальную синусоиду (противофазу). Но произведение таких синусоид всё равно даёт квадрат, хоть и отрицательный (случай 3). Сумма результата даёт большое отрицательное число, что говорит о том, что в сигнале такая частота есть, но в противофазе. Однако в случае, если фаза сдвинута всего на четверть периода, синусоида превращается в косинусоиду, и сумма результата опять обращается в ноль (случай 4), как и в случае отсутствия искомой частоты.

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

Это оказывается настолько удобно, что умножение на синус и косинус одновременно стало стандартной операцией. И для этого даже есть удобный механизм — комплексные числа . При умножении на обычное число, каждое из пары умножается независимо, что нам и нужно. Традиционно, в качестве действительной части берётся косинус, а в качестве мнимой — синус. Такое комплексное число выглядит как $\cos{x}+i\sin{x}$ . Соответственно, можно пользоваться всеми формулами для работы с комплексными числами. Фаза и амплитуда вычисляются как аргумент и модуль комплексной суммы результата.

Более того, эту запись можно ещё упростить. Дело в том, что когда великий российский математик швейцарского происхождения Леонард Эйлер начал исследовать разложение функций в бесконечные ряды, он заметил, что сумма рядов $\cos{x}$ и $i\sin{x}$ совпадает с рядом $e^{ix}$ . Это равенство сейчас называется формула Эйлера:
$e^{ix} = \cos{x}+i\sin{x}$

Каков смысл функции $e^{ix}$ ? Долгое время это оставалось непонятным. Например, в XIX веке профессор Гарвардского университета Бенджамин Пирс после доказательства этой формулы на лекции заявил: "Это, наверное, правда, но она абсолютно парадоксальна; мы не можем понять её, и мы не знаем, что она значит, но мы доказали её, и поэтому мы знаем, что она должна быть достоверной." В настоящее время принята геометрическая интерпретация этой формулы:

$e^{ix}$ — это поворот вектора (1,0) на угол в комплексной плоскости. Например, если $x=\pi$ , что соответствует повороту на 180 градусов, мы получаем вектор (-1,0), что можно также выразить как очень известное тождество Эйлера:
$e^{i\pi} + 1 = 0$ .

Эта формула многими считается самой красивой формулой математики. В ней используются основные математические константы: $0, 1, \pi, e, i$ ; а также основные операции: =, +, *, степень. Если мы прилетим на другую планету с разумной жизнью и посмотрим, как инопланетяне записывают эту формулу, мы начнём понимать существенную часть их математики.

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

Таким образом, если записать в формальном виде сумму произведения дискретного сигнала на синусоиду и косинусоиду, у нас получится:

$\overline{a} = \sum_{n=0}^{N-1}{a_n e^{i2\pi f n/N}}$

Эта формула является мощным инструментом по обнаружению сигнала на частоте . Но что делать, если мы хотим просканировать диапазон частот? Самый простой способ — начать с частоты и идти с некоторым шагом, проверяя каждую (то же самое происходит, когда вы крутите ручку радиоприёмника). Но если мы можем выбирать, какие частоты использовать, то у нас есть весьма заманчивый вариант — частоты от до . Подставив их в предыдущую формулу, мы получим последовательность из результатов:

$\overline{a}_k = \sum_{n=0}^{N-1}{a_n e^{i2\pi k n/N}}, \quad k = 0,...,N-1.$

Если вы заглянете в Википедию в статью «Дискретное преобразование Фурье», то увидите практически идентичную формулу. Основное преимущество выбора таких частот — эффективность вычислений. Благодаря тому, что функция $e^{i2\pi k n/N}$ принимает всего разных значений, их можно выносить за скобки, группировать, использовать различные программистские хитрости, и получить алгоритм, который называется «Быстрое преобразование Фурье» (БПФ, FFT). Этот алгоритм вычисляет результатов из точек исходного сигнала за $N\log{N}$ умножений. В практическом плане это означает, что можно преобразовать сигнал из миллиона точек, получить фазы и амплитуды миллиона частот за 20 миллионов умножений (то есть, за долю секунды на современных вычислительных устройствах). Именно высокая скорость преобразования привела к повсеместному использованию этого алгоритма в нашей повседневной жизни.

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

Мы разобрались с тем, что такое быстрое преобразование Фурье, но что насчёт умножения больших чисел? Для начала, нужно определить задачу. Прежде всего, вспомним, что мы никогда не оперируем числами, а только их представлениями. Когда мы говорим о числе 3196, мы понимаем, что "3196" — это лишь представление этого числа в десятичной системе счисления. C7C₁₆ и 110001111100₂ являются представлениями этого же самого числа в 16-ричной и двоичной системах счисления. Эти представления ничуть не лучше и не хуже, они описывают то же самое число, и их можно использовать там, где это удобно. И у этих представлений есть общее свойство — они используют позиционную систему счисления, то есть цифры могут иметь разное значение в зависимости от их позиции. Мы неявно предполагаем, что каждая цифра умножается на соответствующую степень основания.
$3196 = 3\cdot10^3 + 1\cdot10^2 + 9\cdot10^1 + 6\cdot10^0$ .
Иными словами, цифры являются коэффициентами многочлена

а само число равно значению этого многочлена в точке, равной основанию: .

Таким образом, задача умножения двух чисел сводится к задаче умножения двух многочленов:
$P(x) = A(x) \times B(x)$
Результатом умножения будет значение . Самый простой способ нахождения многочлена преподаётся в начальной школе и называется «умножение столбиком». Каждая цифра одного числа умножается на каждую цифру другого числа, и результат встаёт в позицию, равную сумме позиций цифр. Для примера, вычислим произведение многочленов, соответствующих 3196 и 32:
$P(x) = A(x) \times B(x) = (3x^3 + 1x^2 + 9x + 6)(3x + 2) = 9x^4 + 9x^3 + 29x^2 + 36x + 12$
Результат умножения. Заметим, что мы можем применить к получившемуся многочлену операцию переноса (2 пишем, 1 переносим в следующий разряд, 36 + 1 = 37, 7 пишем, 3 переносим в следующий разряд и т.д.). В результате получится многочлен

Многочлены и являются совершенно разными (даже их степени разные), но имеют одно и то же значение в точке 10: . Операция переноса не является обязательной для получения верного результата, её можно делать или не делать в зависимости от ваших целей. И она выполняется очень быстро, поэтому не сильно влияет на общее время выполнения программы. Умножение столбиком же, наоборот, является очень медленной операцией. Время её выполнения пропорционально квадрату размера чисел.

Существует история, как Андрей Колмогоров на семинаре в 1960 году предположил, что перемножить два числа быстрее, чем за квадрат от их длины, невозможно. В аудитории сидел 23-летний Анатолий Карацуба, который по дороге домой придумал более быстрый алгоритм, известный сейчас как алгоритм Карацубы. Это сломало "квадратную стену" в сознании, и математики начали находить ещё более быстрые алгоритмы. В конце концов, в 2019 году был предложен алгоритм Харви — ван дер Хувена с теоретической сложностью $N\log{N}$ , которая считается оптимальной и неулучшаемой. Но этот алгоритм крайне сложен, и он становится быстрее всех остальных только при перемножении чисел, превышающих размер Вселенной. Поэтому на практике используют другие алгоритмы. Давайте рассмотрим один из них.

В основе этого алгоритма лежит тот факт, что значение многочлена-произведения в произвольной точке равно произведению значений исходных многочленов в этой точке. То есть, для любого
$P(r) = A(r) \cdot B(r)$
Таким образом, даже не зная многочлена , мы можем узнать его значение в любой точке. Но как нам узнать коэффициенты самого ? В этом нам поможет теорема Безу, которая утверждает, что значение многочлена в точке является на самом деле остатком от деления столбиком этого многочлена на многочлен :
$P(r) = P(x) \bmod (x - r)$
Иными словами, если у нас есть значений многочлена в разных точках $\{P(r_k)\}$ , то мы можем считать это множеством многочленов нулевой степени, являющихся остатками от деления на разные модули С помощью китайской теоремы об остатках, которая работает и для многочленов, мы можем сконструировать остаток от деления на произведение всех модулей:

$P_N(x) = P(x) \bmod \prod_{k=0}^{N-1}{(x - r_k)}.$

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

Для примера, возьмём всё те же многочлены и . Посчитаем их значения в точках 1, 2, 3, 4, 5, перемножим и получим значения многочлена в этих точках. Произведение всех модулей равно

Китайская теорема об остатках позволяет нам найти

$P_5(x) = P(x) \bmod (x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120) =\\= 9x^4 + 9x^3 + 29x^2 + 36x + 12.$

Заметим, что этот многочлен совпадает с тем, что получился при умножении столбиком.

Показать вычисления

Точки	1	2	3	4	5
A(x)	19	52	123	250	451
B(x)	5	8	11	14	17
P(x)	95	416	1353	3500	7667
deg 2		416(-1)(x - 3) + 1353(x - 2) = 937x - 1458 mod (x^2 - 5x + 6)		3500(-1)(x - 5) + 7667(x - 4) = 4167x - 13168 mod (x^2 - 9x + 20)
deg 4		((937x - 1458)(1/3x - 6/12) mod (x^2 - 5x + 6))(x^2 - 9x + 20) + ((4167x - 13168)(-1/3x + 22/12) mod (x^2 - 9x + 20)) (x^2 - 5x + 6) = (7286/12x - 13740/12)(x^2 - 9x + 20) + (-5666/12x + 43664/12)(x^2 - 5x + 6) = 135x^3 - 610x^2 + 1422x - 1068 mod (x^4 - 14x^3 + 71x^2 - 154x + 120)
deg 5	95(1/24)(x^4 - 14x^3 + 71x^2 - 154x + 120) + ((135x^3 - 610x^2 + 1422x - 1068)(-1/24x^3 + 13/24x^2 - 58/24x + 96/24) mod (x^4 - 14x^3 + 71x^2 - 154x + 120))(x - 1) = 95(1/24)(x^4 - 14x^3 + 71x^2 - 154x + 120) + (121/24x^3 + 1667/24x^2 - 4382/24x + 11112/24)(x - 1) = 9x^4 + 9x^3 + 29x^2 + 36x + 12 mod (x^5 - 15x^4 + 85x^3 - 225x^2 + 274x^1 - 120)

Вычислительную сложность этого алгоритма можно оценить как $N\log^2{N}$ . Логарифм оказывается в квадрате, потому что алгоритм является рекурсивным — в ходе работы китайской теоремы об остатках нам приходится перемножать многочлены меньшей степени. Но тут есть одна интересная возможность — мы вольны произвольно выбирать точки, в которых вычисляем значения многочленов. И если бы нам удалось подобрать точки так, чтобы все промежуточные многочлены имели вид , то это сильно ускорило бы алгоритм, потому что умножать на многочлен такого вида быстро и просто. Как этого добиться? Очень просто, достаточно вспомнить формулу из школьной алгебры . Возьмём в качестве произведения всех модулей самый простой многочлен и будем раскладывать его с помощью этой формулы до тех пор, пока не получим многочлены первой степени. На первом шаге он раскладывается на $x^{N/2} - 1$ и $x^{N/2} + 1$ . Первый из них раскладывается снова в $x^{N/4} - 1$ и $x^{N/4} + 1$ , а второй в $x^{N/4} - i$ и $x^{N/4} + i$ . Дальше мы получаем $x^{N/8} - 1$ , $x^{N/8} + 1$ , $x^{N/8} - i$ , $x^{N/8} + i$ , $x^{N/8} - \sqrt{i}$ , $x^{N/8} + \sqrt{i}$ , $x^{N/8} - i\sqrt{i}$ , $x^{N/8} + i\sqrt{i}$ . И так далее. Как следует из основной теоремы алгебры, процесс закончится в корнях многочлена . Корни этого многочлена также известны как корни из единицы, они равномерно расположены на единичной окружности в комплексной плоскости:

Тут как раз время вспомнить про функцию $e^{ix}$ . С её помощью мы можем легко записать формулу для корней из единицы:
$r_k = e^{i2\pi k/N}$
Таким образом, теперь можно формально определить, что значит вычислить значение многочлена в точках :

$A(r_k) = \sum_{n=0}^{N-1}{a_n r_k^n} = \sum_{n=0}^{N-1}{a_n e^{i2\pi k n/N }}, \quad k = 0,...,N-1.$

Сравните это с выведенной выше формулой дискретного преобразования Фурье. Да, мы случайно получили именно его. Причём, не похожее преобразование, а абсолютно идентичное. Если посмотреть на программную реализацию без контекста, то нельзя понять, какую именно задачу решал программист — занимался он цифровой обработкой сигнала или же перемножал большие числа. Код получается идентичный. Конечно, этот результат не случаен, если копнуть математику глубже, то можно найти связь между умножением и сигналами, вроде теоремы о свёртке. Но делать этого не стоит, и вот почему.

Дело в том, что дискретное преобразование Фурье совпадает с быстрым умножением только для многочлена . Но это не единственный используемый на практике многочлен. Все числа из списка самых больших простых чисел Прота были найдены с использованием многочлена . Корни этого многочлена также равномерно расположены на единичной окружности, но со сдвигом:

При использовании , умножение по-прежнему остаётся быстрым, алгоритм оказывается почти тем же, изменения минимальны. Но это уже не преобразование Фурье. В этом преобразовании нет нулевой частоты, точек и , которые не дают мнимой части. Все значения (при чётном ) являются комплексными числами. Наиболее близкими к этому алгоритму являются z-преобразование и чирп-алгоритм. Также, этот алгоритм называют «негациклическая свёртка» из-за того, что применение модуля эквивалентно вычитанию коэффициентов, больших , из младших коэффициентов. Но в целом, какого-то определённого собственого названия у алгоритма нет. Более того, его практические реализации часто продолжают называть быстрым преобразованием Фурье, хоть это и сомнительно с математической точки зрения.

Но даже этим разнообразие вариантов быстрого умножения не ограничивается. Программа Cyclo использует корни кругового многочлена для эффективного поиска простых чисел вида $b^{2^n} - b^{2^{n-1}} + 1$ . В данный момент, 7-е и 8-е по величине известные простые числа были найдены с использованием именно этой программы. Уже крайне трудно называть подобные алгоритмы быстрым преобразованием Фурье.

Что мы имеем в итоге? Существуют очень эффективные алгоритмы умножения больших чисел, которые либо совпадают, либо очень похожи на быстрое преобразование Фурье. Однако они имеют совершенно другую, алгебраическую природу. Понимание этого позволяет создавать новые алгоритмы и адаптировать старые для решения других, более интересных задач.

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

Tags:

Hubs: