О длинной арифметике, БПФ, и пользе алгоритмов

Изменено: . Разместил: FadeDemon | Меню
Количество просмотров4`283 / FD1804-6-0

Итак, началось всё с чужого курсового по программированию, темой которого являлась как раз таки длинная арифметика на целых числах, а именно, нужно было всего-то реализовать 4 базовые арифметические операции: сложение, вычитание, умножение, деление. Ну, деление с остатком, понятное дело. С реализацией проблем в общем-то не возникло, для умножения был использован обычный алгоритм умножения в столбик, с квадратичной асимптотикой, для деления классика — сдвиги и вычитания, тоже, вроде-бы работающие за квадрат. Курсовую быстро сдали, и было я забыл уже про неё, но тут посчастливилось познакомиться мне с такой штукой, как алгоритм Диффи — Хеллмана, тот самый, который используется для безопасного обмена ключами. Не знаю, что меня вдруг потянуло в криптографию. Как известно, для практической реализации этого алгоритма, нужно иметь возможность выполнять операции с длинными (порядка 300 знаков) числами, а именно, уметь очень быстро выполнять возведение в степень по модулю. И мне, естественно, было известно про готовые библиотеки, например GMP… Но первым что я вспомнил, оказалось именно та самая курсовая, и ни что иное. Вооружившись напильником, я принялся за работу. И надо сказать, работа была эта очень не простая. Но обо всём по порядку.


Первые опыты


На первый взгляд, имея в арсенале все арифметические операции, написать функцию быстрого возведения в степень было делом плёвым. Однако сразу я столкнулся с проблемой: имеющийся класс работал с числами в десятичной системе счисления (для удобства конвертирования в строку и обратно), на один байт приходилось 10 значений, что сильно затрудняло получения значения определённого бита нашего числа в двоичной записи. А значения битов для классического алгоритма быстрого возведения в степень, как известно, необходимы. Поэтому, первым что пришлось сделать, систему счисления поменять на кратную 2, и я выбрал шестнадцатеричную, то-есть, систему с основанием 16.


Шестнадцатеричная система счисления на самом деле дала помимо возможности проверки битов, еще и тучу вкусностей, например взятие числа по модулю 16 (остатка от деления на 16) можно было заменить операцией побитового AND, и появлялась возможность быстро делить и умножать наше число на числа вида 2n с помощью битовых сдвигов, а поскольку две эти операции очень активно использовались при выполнении арифметических действий, надо-ли говорить о многократном приросте производительности… Но главное, теперь не было проблем написать быстрое возведение в степень. Вот, кстати, что у нас получилось:


Код: C++
//беззнаковое помножалко в столбик
void __fastcall CHeavyNumber::op_mul(const CHeavyNumber &x){
	//тут вот все чуть сложнее
	long oll = max(nlen, x.nlen) * 2 + 1;
	this->set_length(oll);
	long i, j, k, l, ai, tmp;

	//тут придется использовать временный буфер
	unsigned char *bui = (char *) malloc(oll);
	memset(bui, 0, oll);

	//типо такое умножение в стлобик
	for(i = 0; i < x.nlen; i++){
		for(j = 0; j < nlen; j++){
			tmp = x.nbuf[i] * this->nbuf[j];
			k = 0;

			while(tmp){

				bui[i + j + k] += tmp & 15;
				l = 0;

				while(bui[ai = i + j + k + l] >= 16){

					bui[ai + 1] += bui[ai] >> 4;
					bui[ai] &= 15;
					l++;
				}

				tmp >>= 4;
				k++;
			}
		}
	}

	memcpy(nbuf, bui, oll);
	free(bui);

	this->trim_right();
}

void CHeavyNumber::qpowm(const CHeavyNumber &n, const CHeavyNumber &p){
	if(n.sign || p.sign) return;

	CHeavyNumber x(*this);
	*this = 1;
	long ma = n.nlen << 2;
	long i;

	for(i = 0; i < ma; i++){
		//как можем заметить мы, в 16ричной системе проверка битов вопросов не вызыввает
		if(n.nbuf[i >> 2] & (1 << (i & 3))){
			*this *= x;
			*this %= p;
		}
		x *= x;
		x %= p;
	}
}

Ну, тут стоит пояснить что переменная класса nbuf — указатель на массив из цифр числа, записанных в порядке от младших к старшим. В остальном всё тривиально, все операции здесь должным образом перегружены. Казалось-бы, наша проблема на этом решена, но когда я попытался с помощью этой функции возвести число порядка 10100 в степень около 10300 по такому-же модулю… Я, в общем, понял, что на самом деле всё здесь только начинается. Возведение выполнялось почти 3 минуты, я даже не стал ждать конца…


Начинаем ускоряться


Конечно первой реакцией на такой расклад был ступор, хотя, в общем-то, это было довольно предсказуемо. Вся проблема заключалась в алгоритмах, взять хотя-бы умножение в столбик, ведь это по сути самый медленный алгоритм из имеющихся. Но как его ускорить? Беглый поиск по проблеме открыл для меня такую интересную штуку, как алгоритм умножения Карацубы, работающий за O(Nlog23). Но мной был избран более интересный путь. Краем уха я слышал про метод умножения многочленов (а значит и длинных чисел) с помощью великого и ужасного быстрого преобразования Фурье, который работал за ещё лучшую асимптотику — O(N ⋅ log(N)). В общем, решено было остановиться на нём.


Для начала попробуем коротко пояснить, что вообще это такое, преобразование Фурье. Если пытаться объяснять простым языком, то это такая операция, которая из некоторого сигнала, или функции, нигде не стремящейся в бесконечность, извлекает спектр этого сигнала, или выражаясь по другому, раскладывает его на частоты, именуемые ещё гармониками. Именно этот процесс можно пронаблюдать в большинстве медиаплееров, которые при проигрывании звука рисуют диаграмму его частот, “эквалайзер”, как некоторые его называют. Диаграмма эта как раз таки строиться с помощью быстрого преобразования Фурье, или сокращённо, БПФ.


Сразу нужно заметить, что как медиаплеер, так и мы, работать будем с дискретными сигналами, то-есть такими, которые можно задать с помощью ограниченного количества измеренных его значений, идущих через определённый промежуток времени. В случае медиаплеера, этими значениями будут единичные импульсы в PCM модуляции, идущие со скоростью частоты дискретизации на устройство вывода звука, и превращающийся за счет этого в, собственно, звук, ну а в нашем случае этими значениями будут одинокие циферки нашего длинного числа, или по другому, коэффициенты нашего многочлена. Ведь думаю не нужно долго разъяснять, что длинное число из цифр a1 … an с основанием системы счисления x можно представить в виде многочлена примерно таким вот образом:


\[ P_n(x) = a_1x^0 + a_2x^1 + a_3x^2 + \mathellipsis + a_nx^{n-1} \]
Многочлен степени n

Чтобы это можно было понять, возьмём для примера число 435 в десятичной системе счисления. Цифры у нас здесь таковы: a1 = 5, a2 = 3, a3 = 4, идут они всегда от младших разрядов к старшим. По-скольку система счисления десятичная, то x = 10. Подставим теперь всё это в формулу выше, и получим:


\[ 5 \cdot 10^0 + 3 \cdot 10^1 + 4 \cdot 10^2 = 5 + 3 \cdot 10 + 4 \cdot 100 = 5 + 30 + 400 = 435 \]


Что, в общем-то и требовалось получить. Таким образом, в нашем случае, число будет представляться как такой-же многочлен, только с другим основанием — x = 16. Теперь нам нужно разобраться, что мы собираемся делать с нашими дискретными цифрами, чтобы получить из них спектр частот. А главное — зачем. Но по порядку обо всём.


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


Чтобы понять на пальцах, в чём же суть ДПФ, нужно вспомнить одну важную вещь из комплексного анализа, а именно, что комплексных корней степени n из единицы, может существовать ровно n штук. Например, квадратные корни из 1, это, очевидно, ±1, а корни 4-й степени из 1, это ±1 и ±i. В школе, все, наверное, извлекали корни третьей степени из единицы. Ну так вот, на самом-же деле, дискретным преобразованием Фурье коэффициентов некоторого многочлена степени n являются n комплексных значений этого многочлена во всех точках комплексной плоскости, совпадающих с корнями степени n из единицы. Важно здесь то, что коэффициентами этими может быть, в принципе, любой дискретизированный сигнал, в том числе и цифры нашего длинного числа. Многочлен мы здесь как-бы додумываем.


Двигаясь дальше, нужно вспомнить ещё две вещи — формулу n-го комплексного корня из единицы, ну и тот факт, что возводя в некоторую степень один из этих корней (так называемый, главный корень) можно быстро получить все остальные. Обозначив за wn,k k-тый корень степени n из единицы, за wn — главный корень, формулы для корней будут выглядеть следующим образом:


\[ w_n = e^{i{2\pi \over n}} \;\;\;\; w_{n,k} = e^{i{2\pi k \over n}} \;\;\;\; w_{n,k} = \left (w_n\right )^k \] \[ n,k \in \mathbb{Z} \;\;\;\; w_{n,k} \in \mathbb{C}\]
Свойства комплексных корней из единицы

Помня эти простые правила и обозначения, легко записать, что собственно представляет собой ДПФ. Допустим, у нас есть многочлен Pn(x) степени n, который мы записывали выше, и у него, соответственно, есть коэффициенты a1 … an, тогда ДПФ вектора его коэффициентов будет представлять собой вектор, состоящий из комплексных значений, и определённый следующим образом:


\[ {\rm DFT} \left [a_1, a_2, a_3, \; \mathellipsis \; a_n \right ] = \left [P(w_n^0), P(w_n^1), \; \mathellipsis \; P(w_n^{n-1}) \right ] \]
ДПФ для вектора коэффициентов многочлена P

Сразу надо пояснить физический смысл получаемых значений. Как мы уже знаем, преобразование Фурье раскладывает сигнал на частотные составляющие, а каждая составляющая описывается тремя характеристиками: частотой, амплитудой (интенсивностью), и фазой (смещением). Полученные комплексные значения с этими характеристиками связаны следующим простым образом: если комплексное значение имеет индекс k, то тогда его модуль делённый на n — значение амплитуды составляющей сигнала с частотой k/T (T — промежуток времени, за который были получены входные данные), а его аргумент — фаза этой составляющей. Зная это, уже достаточно легко построить, например, график спектра частот звукового сигнала, который строит нам Winamp. Впрочем именно эти знания нам по некоторой причине не пригодятся. Мы лучше возьмём, и запишем наконец-то окончательную формулу вычисления ДПФ. Пусть на входе у нас будет тот-же вектор a0 … an−1, а на выходе мы собираемся получить вектор b0 … bn−1, тогда для значения bk будет справедливо равенство:


\[ b_j = \sum\limits_{k=0}^{n-1} a_k \left (e^{i{2\pi \over n}}\right )^{kj} \]
Прямое ДПФ

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


\[ a_j = {1 \over n}\sum\limits_{k=0}^{n-1} b_k \left (e^{i{2\pi \over n}}\right )^{-kj} \]
Обратное ДПФ

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


Умножение многочленов через ДПФ


Итак, в чём-же суть нашего быстрого умножения? А суть, дорогие товарищи, очень проста. Оказывается, что результат перемножения двух многочленов, это есть ни что иное, как их линейная свёртка, а вычисление линейной свёртки двух многочленов эквивалентно обычному поэлементному перемножению их частотной области, то-есть, поэлементному перемножению комплексных значений их Фурье-образов. В нашем случае, остаётся только в полученном результате сделать обыкновенный перенос разрядов по модулю системы счисления, и всё. Результат готов. Красота, не правда-ли?


Остаётся рассмотреть алгоритм более подробно, отметив все тонкости. Во первых, многочлены нужно дополнить нулями до длины бо́льшего из них, умноженной на 2, поскольку только в этом случае, в результирующем многочлене с такой-же длиной 100% уместиться сколь угодно большой результат перемножения. Далее, от обоих многочленов нужно вычислить ДПФ, и поэлементно их умножить, то-есть результирующий элемент с индексом k будет равен произведению элементов с индексами k из первого и второго векторов. Затем, от того что получилось, нужно вычислить обратное ДПФ, взять его действительные значения, и сделать перенос разрядов. Всё. Формально, имея многочлены Pn(x) и Qn(x), это можно записать примерно так:


\[ P_n (x) \times Q_n (x) = {\rm DFT}^{-1}({\rm DFT}(P_{2n}) \cdot {\rm DFT}(Q_{2n})) \]
Умножение многочленов с помощью ДПФ

Здесь, как можно догадаться, DFT−1 означает обратное ДПФ, × имеет смысл полного перемножения, а ⋅ (точка) — поэлементного (скалярного). Описанный алгоритм, кстати, носит название “алгоритм Шёнхаге — Штрассена”, считается одним из самых быстрых алгоритмов перемножения длинных чисел, и имеет асимптотику O(N ⋅ log(N)). Правда пытливый читатель сразу заметит, что когда мы рассматривали ДПФ, то там логарифмическим временем и не пахло, формулы ДПФ и обратного ДПФ имели квадратичный порядок времени. Но, если честно, то ДПФ в том виде, в котором мы его рассмотрели — это не совсем то, что нам нужно. И именно поэтому мы идём дальше…


Быстрое преобразование Фурье, и его оптимизация


Задача ускорения ДПФ, на самом деле, сводиться к применению стандартной стратегии “разделяй и властвуй”, которая за счёт рекурсивного разбиения имеющейся задачи на более мелкие подзадачи позволяет существенно снижать асимптотическую оценку для алгоритма, обычно сводя её к логарифмической — O(N ⋅ log(N)) для различных сортировок, и нашего БПФ, O(log(N)) для, к примеру бинарного поиска. Давайте теперь разберёмся, как это будет происходить в нашем случае.


Снова возьмём наш многочлен Pn(x), с коэффициентами a0 … an−1, и сделаем из этих коэффициентов два других многочлена, An/2(x) и Bn/2(x). Сразу надо сказать, что с этого места будем считать, что n у нас всегда есть степень двойки, то-есть, n = 2m, где m — некоторое целое положительное число. Многочлены наши будут выглядеть следующим образом:


\[ \begin{aligned} & A_{n \over 2}(x) = a_0 x^0 + a_2 x^1 + a_4 x^2 + \mathellipsis + a_{n - 2} x^{{n \over 2}-1} \\ & B_{n \over 2}(x) = a_1 x^0 + a_3 x^1 + a_5 x^2 + \mathellipsis + a_{n - 1} x^{{n \over 2}-1} \\ \end{aligned} \]


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


\[ P_n(x) = A_{n/2}(x) + x B_{n/2}(x) \]


Теперь, предположим, что мы уже вычислили DFT(A) и DFT(B), и попробуем с помощью них вычислить интересующую нас DFT(P). Обозначим коэффициенты многочлена An/2(x) как c0 … cn/2−1, коэффициенты Bn/2(x) как d0 … dn/2−1, а значения интересующего нас DFT(P) как y0 … yn−1, и тогда, поглядев на формулы выше, мы уже можем, без больших усилий, получить формулы, так называемого преобразования бабочки, с помощью которых уже не составит никакого труда, рекурсивно определить БПФ. Формулы эти выглядят следующим образом:


\[ \begin{aligned} & y_j = c_j + w_n^j d_j, & j \in [0; {n \over 2} - 1] \\ & y_{j+{n \over 2}} = c_j - w_n^j d_j, & j \in [0; {n \over 2} - 1] \\ \end{aligned} \]
Преобразование бабочки

Вывод второй формулы мы не будем приводить, надо только заметить, что в нём используются такие тождества, как:


\[ \left (w_n\right )^n = 1 \;\;\;\; \left (w_n\right )^{n \over 2} = -1 \]
Тождества корней из единицы

Теперь, не составит никакого труда определить окончательный алгоритм быстрого преобразования Фурье, с использованием рекурсии, если учесть тот факт, что ДПФ от вектора, содержащего только 1 элемент, этим-же вектором и является, то-есть для его вычисления делать ничего не надо. Попробуем это записать:


\[ \left \lbrace \begin{aligned} & \mathfrak{F}\left [a_0\mathellipsis a_{n-1}\right ] = \underbrace{\mathfrak{F}\left [a_{j}\right ] + \left (w_n\right )^{j \over 2} \cdot \mathfrak{F}\left [a_{j + 1}\right ] \cup \mathfrak{F}\left [a_{j}\right ] - \left (w_n\right )^{j \over 2} \cdot \mathfrak{F}\left [a_{j + 1}\right ]} \limits_{j = 0, 2, 4 \; \mathellipsis \; n - 2} \\ & \mathfrak{F}\left [a\right ] = \left [a\right ] \end{aligned} \right . \]
Быстрое преобразование Фурье в рекурсивной форме

Вот такая вот устрашающая получилась у нас формула. Здесь под значком объединения ∪ надо подразумевать простую конкатенацию двух векторов, то-есть их склейку, ну а под готичной F — сокращение для FFT. Однако даже сейчас рано переходить к написанию работающего кода. Всё дело в том, что БПФ в рекурсивном виде является далеко не самым оптимальным решением, так как вынуждает нас на каждом шаге копировать элементы входных данных в момент их разделения. Поэтому, от рекурсии в явной форме мы откажемся, и рассмотрим ещё одну вещь: поразрядно обратную перестановку.


Чтобы понять, для чего эта перестановка нужна, давайте проследим, на какие именно векторы последовательно разбивает рекурсия исходный вектор данных. Допустим, если взять вектор A из 8-ми элементов a0 … a7, то эти разбиения будут для него выглядеть следующим образом:


\[ \overline{A} = \left [ \left [ \left [ a_0, a_4 \right ]_2, \left [ a_2, a_6 \right ]_2 \right ]_1, \left [ \left [ a_1, a_5 \right ]_2, \left [ a_3, a_7 \right ]_2 \right ]_1 \right ]_0 \]


Индексами после скобок здесь обозначены уровни рекурсии, на которых это разбиение возникает. И здесь хорошо видно, что при такой записи, два соседних вектора последнего уровня рекурсии, вернее ДПФ от них, являются входными данными для предыдущего уровня рекурсии, а два соседних 4-х элементных вектора, или ДПФ от них, входные данные для еще более высокого уровня рекурсии, и так далее. Получается, что если бы мы могли переставить элементы вот так, как здесь нарисовано, то мы бы смогли, во первых, избавиться от рекурсии в явном виде, вычисляя БПФ итерационно, а во вторых, обходиться без дополнительной памяти, поскольку результаты всегда бы вставали туда, где они должны быть в итоге. То-есть, после такой перестановки, для вектора из 8-ми элементов, вычисляться сначала БПФ для подвекторов [a0,a1], [a2,a3], [a4,a5], [a6,a7], причём результаты останутся там же, где были и входные данные, затем вычисляться БПФ для [a0 … a3] и [a4 … a7], и в самом конце преобразуется уже весь вектор [a0 … a7]. А вот здесь, как раз и приходит на помощь поразрядно обратная перестановка: оказывается, что если переставить элементы местами так, чтобы каждый элемент с индексом k обменялся местом с элементом с индексом bit_reverse(k), то мы, собственно, требуемую перестановку элементов и получим. Операция bit_reverse, или операция реверса битов определяется довольно просто: нужно число записать в двоичной системе, и слева дополнить его нулями так, чтобы полученное количество цифр могло вместить любое число вплоть до n − 1, где n это, собственно, длина входных данных, (здесь, кстати, легко заметить, что числа вида 2n − 1, по совместительству являющиеся числами Мерсенна, будут состоять в двоичной записи из одних единиц), и биты в получившейся записи переставить в обратном порядке. Допустим, если число k у нас равно 5 или 01012, n = 16, то bit_reverse(k) примет значение 10102, или 10.


Переходя от формул к коду


Ну вот, наконец, только теперь, можно попробовать что-нибудь реализовать. Вообще, когда реализовывал БПФ я, то мой код претерпел на себе OVER 9000 переписываний и допиливаний, поэтому выложу сразу окончательную версию, переписанную в бо́льшей половине своей на ассемблер с SSE. Сразу надо заметить, что там, где не используется ассемблер — используются обычные массивы, и рукотворные комплексные числа, представляющие собой обыкновенные структуры с перегруженными операторами. В ассемблерной части главную роль играет эпичная инструкция ADDSUBPD, значительно ускоряющая процедуру умножения комплексных чисел, в паре с инструкциями UNPCKHPD и UNPCKLPD. Вот, наконец, как всё это в итоге стало выглядеть:


Код: C++
extern long *bit_rev_table[bit_rev_table_size];

//наконец-то. встречайте, быстрое преобразование Фурье...
//а теперь еще и с поддержкой SSE2! эгегее!
bool fftcx(TFFComplex *dest, long n, bool inv){

	if(1 == n) return true;

	long l2;

	//быстрый логарифм по основанию 2
	ffm_flog2(n, &l2);

	if(!l2) return false;

	TFFComplex *wroots = (TFFComplex *) malloc_a16(sizeof(TFFComplex) * n);

	long i, j, jl, jl2;
	long *ttb;
	ttb = bit_rev_table[l2];

	__asm {
		movd mm0, esi
		movd mm1, edi
		movd mm2, ebx

		mov esi, dest
		mov edi, ttb

		mov ecx, n
		xor edx, edx

		reswaper:

			cmp edx, [edi]
			jnl noswap
			jz noswap

				mov eax, edx
				sal eax, 4

				mov ebx, [edi]
				sal ebx, 4

				movaps xmm0, [esi + eax]
				movaps xmm1, [esi + ebx]

				movaps [esi + eax], xmm1
				movaps [esi + ebx], xmm0

			noswap:

			add edi, 4

		inc edx
		dec ecx
		jnz reswaper

		movd ebx, mm2
		movd edi, mm1
		movd esi, mm0

		emms
	}

	// альтернатива без ассемблера
	// for(i = 0; i < n; ++i) dest[i] = src[i];
	// for(i = 0; i < n; ++i) if(i < bit_rev_table[l2][i]) swap(dest[i], dest[bit_rev_table[l2][i]]);

	double ang;

	TFFComplex *wlen = (TFFComplex *) malloc_a16(sizeof(TFFComplex));

	for(jl = 2; jl <= n; jl <<= 1){
		ang = 2 * M_PI / (double) jl * (inv ? -1 : 1);
		jl2 = jl >> 1;

		wlen->re = cos(ang);
		wlen->im = sin(ang);

		wroots[0].re = 1;
		wroots[0].im = 0;

		__asm {
			push esi
			mov esi, wroots
			mov ecx, jl2
			dec ecx
			test ecx, ecx
			jz okeento

			mov eax, wlen
			movaps xmm6, [eax]

			multiplento:

				movaps xmm1, [esi]
				movaps xmm2, xmm6

				// xmm1,xmm2 = input
				movaps xmm3, xmm1
				mulpd xmm1, xmm2
				shufpd xmm3, xmm3, 1
				mulpd xmm3, xmm2

				movaps xmm2, xmm1
				unpckhpd xmm1, xmm3
				unpcklpd xmm2, xmm3

				// xmm2 = result
				addsubpd xmm2, xmm1

				movaps [esi + 0x10], xmm2

                add esi, 0x10

			dec ecx
			jnz multiplento

			okeento:

			pop esi
		}

		/*for(i = 1; i < jl2; ++i){
			wroots[i] = wroots[i - 1] * *wlen;
		} */

		for(i = 0; i < n; i += jl){
			__asm {
				push edi
				push esi

				mov esi, wroots

				mov ecx, i
				sal ecx, 4
				add ecx, dest

				mov edx, jl2
				sal edx, 4
				mov eax, edx
				add edx, ecx

				// преобразование бабочки, однако
				loop_butterfly:

					// xmm0 = U
					movaps xmm0, [ecx]
					movaps xmm1, [ecx + eax]
					movaps xmm2, [esi]

					// умножение комплексных чисел
					movaps xmm3, xmm1
					mulpd xmm1, xmm2
					shufpd xmm3, xmm3, 1
					mulpd xmm3, xmm2

					movaps xmm2, xmm1
					unpckhpd xmm1, xmm3
					unpcklpd xmm2, xmm3

					// xmm2 = V
					addsubpd xmm2, xmm1
					movaps xmm1, xmm0

					addpd xmm0, xmm2
					subpd xmm1, xmm2

					movaps [ecx], xmm0
					movaps [ecx + eax], xmm1

					add ecx, 0x10
					add esi, 0x10

					cmp ecx, edx

					jc loop_butterfly

				pop esi
				pop edi
			}

			/*for(j = 0; j < jl2; ++j){
				u = dest[i + j];
				v = dest[i + j + jl2] * wroots[j];

				dest[i + j] = u + v;
				dest[i + j + jl2] = u - v;
			}*/
		}
	}

	if(inv){

		__asm {
			push esi
			mov esi, dest
			mov ecx, n
			cvtsi2sd xmm1, ecx
			shufpd xmm1, xmm1, 0

			test ecx, ecx
			jz nodivo

			redivo:

				movaps xmm0, [esi]
				divpd xmm0, xmm1
				movaps [esi], xmm0
                add esi, 0x10

			dec ecx
			jnz redivo

			nodivo:

			pop esi
		}

		// for(i = 0; i < n; ++i) dest[i] /= n;
	}

	free_a16(wlen);
	free_a16(wroots);
	return true;
};

Конечно, здесь требуются кое-какие разъяснения. Ну, во первых, в комментариях остались кусочки альтернативной реализации без ассемблера. Как мы можем заметить, ассемблер раздувает объём кода раз наверно в 10… Во вторых, переменная bit_rev_table — это таблица предподсчитанных значений реверса битов для всех возможных случаев. Если окунуться в детали, то это указатель на массив указателей, которые показывают на массивы значений bit_reverse для всех возможных значений для заданной длинны числа в двоичной записи, или для заданной степени двойки, кстати именно степень двойки является индексом в массиве указателей. Далее, функция ffm_flog2 здесь быстро вычисляет целочисленный логарифм по основанию 2 с округлением вниз, при помощи поистине волшебных инструкций BSR и BSF, вот её код:


Код: C++
//считает степень 2, или ноль, если это не степень 2
void __fastcall ffm_flog2(long a, long *l2){
	__asm {
		xor ecx, ecx
		bsr edx, a
		bsf eax, a
		cmp eax, edx
		cmove ecx, edx
		mov edx, l2
		mov [edx], ecx
	}
}

Кстати здесь используется ещё интересная инструкция CMOVE, которая копирует значение источника в приёмник, только если установлен флаг ZF. Дальше мы рассмотрим функцию malloc_a16, которая делает ни что иное, как выделяет блок памяти, выровненный по границе в 16 байт, при использовании SSE это довольно распространенная и часто геморройная задача… Особенно если учесть, что мне так и не удалось заставить компилятор C++ в моей Embarcadero RAD Studio XE2 выравнивать переменные по заданной границе, через __declspec(align(16))… Ну да впрочем не важно. Мне пришла в голову идея дописать некоторую надстройку над имеющимся менеджером памяти, вот в таком виде:


Код: C++
//----------------------------------------------------------------------------//

void * __fastcall malloc_a16(size_t count){
	void *temp = malloc(count + 0x20);
	if(NULL == temp) return NULL;
	void *alt = (void *)((unsigned long) temp + 0x20 & ~0xf);
	*(void **)((unsigned long)alt - 0x10) = temp;
	*(size_t *)((unsigned long)alt - 0xc) = count;
	return alt;
}

void __fastcall free_a16(void *ptr){
	free(*(void **)((unsigned long)ptr - 0x10));
}

void * __fastcall realloc_a16(void *ptr, size_t count){
	size_t blosi = *(size_t *)((unsigned long)ptr - 0xc);
	if(blosi < count) return ptr;

	void *temp = malloc_a16(count);
	if(NULL == temp) return NULL;

	memcpy(temp, ptr, blosi);
	memset((void *)((unsigned long) temp + blosi), 0, count - blosi);

	free_a16(ptr);

	return temp;
}

//----------------------------------------------------------------------------//

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


Ещё один грязный хак, который здесь используется, связан с тем, что когда проект компилировался с максимальным уровнем оптимизации, то в больших кусках кода, вместо того чтобы к переменным обращаться с помощью регистра EBP, компилятор без лишних рассуждений адресовался напрямую через ESP, что в итоге, при наличии в коде команд PUSH/POP приводило к полной белиберде. Выход из этой ситуации я нашел в виде использования вместо стека регистров MMX, засовывая туда значения регистров с помощью команды MOVD, и вытаскивая таким-же образом обратно. Именно поэтому, в коде выше кое-где красуются команды EMMS.


Про TFFComplex рассказывать мы не будем, поскольку это обычная структура из двух double, с перегруженными соответствующим образом операторами. Лучше сразу перейдём к коду, который был призван использовать вышеиспечённое БПФ применительно к умножению наших длинных чисел. Вот, собственно, и оно:


Код: C++
//долгожданное умножение с помощью БПФ
void CHeavyNumber::op_mul_fft(const CHeavyNumber &x){
	TFFComplex *a, *b;
	long max_sz = max(nlen, x.nlen);

	//степеньдвойкизируем
	__asm {
		bsr ecx, max_sz
		xor eax, eax
		inc eax
		shl eax, cl

		notgood:

			cmp eax, max_sz
			jge okay

				shl eax, 1
				jmp notgood

		okay:

		shl eax, 1
		mov max_sz, eax
	}

	a = (TFFComplex *) malloc_a16(max_sz * sizeof(TFFComplex));
	b = (TFFComplex *) malloc_a16(max_sz * sizeof(TFFComplex));

	memset(a, 0, sizeof(TFFComplex) * max_sz);
	memset(b, 0, sizeof(TFFComplex) * max_sz);

	long i;

	for(i = 0; i < this->nlen; ++i) a[i].re = this->nbuf[i];
	for(i = 0; i < x.nlen; ++i) b[i].re = x.nbuf[i];

	fftcx(a, max_sz);
	fftcx(b, max_sz);

	//for(i = 0; i < max_sz; ++i) a[i] *= b[i];
	__asm {
		push esi
		push edi

		mov esi, a
		mov edi, b

		mov ecx, max_sz

		multiplento:

			movaps xmm1, [esi]
			movaps xmm2, [edi]

			//xmm1,xmm2 = input
			movaps xmm3, xmm1
			mulpd xmm1, xmm2
			shufpd xmm3, xmm3, 1
			mulpd xmm3, xmm2

			movaps xmm2, xmm1
			unpckhpd xmm1, xmm3
			unpcklpd xmm2, xmm3

			//xmm2 = result
			addsubpd xmm2, xmm1

			movaps [esi], xmm2

			add esi, 0x10
			add edi, 0x10

		dec ecx
		jnz multiplento

		okeento:

		pop edi
		pop esi
	}
	//получили линейную свертку

	fftcx(a, max_sz, true);

	this->set_length(max_sz);

	long mam1 = max_sz - 1;
	long t = 0;

	//вот это дерьмо очень было бы неплохо очеловечить
	for(i = 0; i < mam1; ++i){
		t = (long) floor(a[i].re + 0.5);
		a[i + 1].re += t >> 4;
		//a[i].re = t & 15;
		this->nbuf[i] = t & 15;
	}

	nbuf[mam1] = floor(a[mam1].re + 0.5);
	//for(i = 0; i < max_sz; ++i) nbuf[i] = (long) floor(a[i].re + 0.5);

	//normalize();
	this->trim_right();

	free_a16(a);
	free_a16(b);
}

Здесь, в общем-то довольно ярко прослеживаются все действия, необходимые для перемножения многочленов. Первая ассемблерная вставка занимается тем, что добивает размер до ближайшей сверху степени двойки, и увеличивает его на 2. Далее, собственно, выделяется память, и копируются значения цифр, в шестнадцатеричной, кстати, системе счисления, которую использует наша библиотека. Затем вычисляются БПФ, вторая ассемблерная вставка их с помощью SSE перемножает, вычисляется обратное БПФ, и производиться поразрядный перенос, не совсем идеально, кстати, оптимизированный, хотя тут это уже мало значимо.


Результатом всех этих усилий стало то, что два числа, умножаемые классическим алгоритмом умножения в столбик за 35 миллисекунд, умножались кодом выше за 4 миллисекунды. На более длинных числах разница была больше. На наше же возведение в степень код оказал тоже довольно ощутимое влияние: время работы упало с нескольких минут до 8-10 секунд. Смешно конечно, но прогресс определённо есть. Все дальнейшие усилия были уже направлены на оптимизацию деления, а точнее взятия остатка от деления. Впрочем в этой области сильно далеко уйти не удалось…


В поисках быстрого деления


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


Код: C++
//хитрожопое деление со всеми премудростями
void CHeavyNumber::op_divide(const CHeavyNumber &x, CHeavyNumber &r){
	long comp = this->op_cmp(x);

	if(comp < 0){
		r = *this;
		this->init_from_long(0);
		return;
	}

	CHeavyNumber z(-0), t(-0), d(x), tmp(-0);

	long led = 0;

	z = *this;
	z.sign = false;
	d.sign = 0;

	led = z.nlen - x.nlen;
	d = x << led;
	tmp = 1;
	tmp <<= led;

	//обязательно сраывнение без учота знака
	while(z.op_cmp(x) >= 0){
		while(z >= d){
			t += tmp;
			z -= d;
		}

		tmp >>= 1;
		d >>= 1;
	}

	r = z;
	bool si = this->sign;
	this->assign(t);
	this->sign = si;
}

Конечно, для быстрого взятия остатка здесь было довольно много всего лишнего, например вычисление частного, которым можно было пренебречь. В результате обработки напильником, получилась в итоге вот такая функция:


Код: C++
//ориентировано на максимально быстрое получение остатка
void __fastcall CHeavyNumber::op_divmod(const CHeavyNumber &x){
	long comp = op_cmp(x);

	if(comp < 0) return;

	CHeavyNumber z(-0), d(x);

	long led = 0;

	z = *this;
	z.sign = false;
	d.sign = false;

	led = z.nlen - x.nlen;
	d = x << led;

	while(z.op_cmp(x) >= 0){
		while(z.op_sub(d, true));

		d >>= 1;
	}

	*this = z;
}

Здесь сто́ит по подробнее остановиться на нюансах. Для начала нужно прояснить тот факт, что операции сдвига в общем-то делают то-же самое что и обычные сдвиги, только сдвигают число сразу на 4 бита, то-есть во внутреннем представлении числа это получается как раз 1 байт. Приведём заодно и код этих сдвигов, уже, кстати, местами оптимизированный разбавленный ассемблером:


Код: C++
//сдвиги по модулю (умножалово на 16 или 1/16)
void __fastcall CHeavyNumber::op_shift_left(long x){
	//умножение на 16
	if(is_zero()) return; //сколько ноль не умножай...
	if(x <= 0) return;
	this->set_length(nlen + x);
	char *si = nbuf + nlen - 1;
	char *su = si - x;
	do{
		*(si--) = *su;
		*(su--) = 0;
	}while(su >= nbuf);
}

void __fastcall CHeavyNumber::op_shift_right(long x){
	//деление на 16
	if(nlen <= x){
		this->init_from_long(0);
		return;
	}

	if(x <= 0) return;

	/*char *si = nbuf;
	char *su = nbuf + x;
	char *sa = nbuf + nlen - 1;
	do{
		*(si++) = *(su++);
	}while(si <= sa);*/

	long loto = nlen;
	unsigned char *lizo = nbuf;

	// немного убыстренная версия
	__asm {
		cld

		movd mm0, edi
		movd mm1, esi

		mov esi, lizo
		mov edi, esi

		mov ecx, loto
		sub ecx, x
		add esi, x
		mov edx, ecx

		sar ecx, 2

		rep movsd

		mov ecx, edx
		and ecx, 3

		rep movsb

		movd esi, mm1
		movd edi, mm0

		emms
	}

	//ненужное откусываем
	this->set_length(nlen - x);
}

Про операцию же op_sub, которая здесь представляет собой обычное вычитание, придётся рассказать чуть более подробно, поскольку именно она являлась в данном алгоритме чрезвычайно узким местом, так как выполнялась просто огромное количество раз, примерно log2x раз, где x — делимое. Сначала внесём ясность в op_cmp, эта функция, как следует из названия, просто сравнивала два числа, и возвращала −1 если первое число меньше, 1 — если больше, и 0 в случае равенства. Код очень простой:


Код: C++
//обычная сравнивалка, работающая по обычным правилам
long __fastcall CHeavyNumber::op_cmp(const CHeavyNumber &x) const{
	if(nlen != x.nlen){
	   return nlen < x.nlen ? -1 : 1;
	}

	long i;
	for(i = nlen - 1; i >= 0; i--){
		if(nbuf[i] != x.nbuf[i]){
			return nbuf[i] < x.nbuf[i] ? -1 : 1;
		}
	}
	//ну равны они
	return 0;
}

Кстати как я потом понял, здесь даже оптимизировать нечего, поскольку вероятность выполнения n итераций приведённого здесь цикла убывает по закону 161−n, так как именно этот закон описывает вероятность совпадения n + 1 подряд идущих с конца символов двух случайных чисел в шестнадцатеричной системе счисления. Совсем другие просторы для оптимизации открываются в функции op_sub, то-есть вычитания. На самом деле, совокупность таких факторов, как система счисления с основанием 16, и размещение одной цифры на один байт позволило для вычитания (и, кстати, сложения) привлечь инструкции расширенного MMX, входящие в набор SSE2, которые позволяли одновременно вычитать сразу по 16 байт за один присест, с помощью команды PSUBB, а после некоторых ухищрений, примерно с такой-же скоростью, выполнять и поразрядный перенос. Попробуем разобрать, как этого удалось достичь:


Код: C++
//вычиталко
bool __fastcall CHeavyNumber::op_sub(const CHeavyNumber &x, bool is_unneg){
	//здесь надо быть хитрожопыми, и всегда вычитать из большего
	//меньшее. тогда будет кошерно, результат корректировать с помощью знака
	//рациональнее всего будет использовать тут, короче, операции сравнения
	//а оные работают не сложнее, чем ассемблерная repz cmpsb, только в обратную сторону
	//другой вариант - в случае получения "отрицательного кода" - вычитать число из 999999999999...
	long comp = this->op_cmp(x);

	//не вычитаем, если вдруг число второе меньше
	if(comp < 0 && is_unneg) return false;

	long oll = max(this->nlen, x.nlen);
	long sdl = min(oll, x.nlen);
	this->set_length(oll);
	long i;

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

	//тут все сразу понятно
	if(0 == comp){
		// x - x = 0
		this->init_from_long(0);
		return true;
	}

	/*if(comp < 0){
		for(i = 0; i < oll; i++) nbuf[i] = (x.nbuf[i] - nbuf[i]);
	}else{
		for(i = 0; i < oll; i++) nbuf[i] = (nbuf[i] - x.nbuf[i]);
	}*/


	char *bu1 = nbuf, *bu2 = x.nbuf;

	long sdss = 5;

	//без ассемблера никуда
	__asm {
		//push edi
		//push esi
		cld

		movd mm0, edi
		movd mm1, esi

		mov esi, bu1
		mov edi, bu2

		cmp comp, 0

		jge hessy
			xchg esi, edi
		hessy:

		mov ecx, sdl
		mov eax, ecx

			mov ecx, eax
			sar ecx, 4

		test ecx, ecx
		jz shory1

			//убираем обработаные биты
			and eax, 0xf

			//грузим нужные регистры
			mov edx, hvn_a_reg_0
			movaps xmm7, [edx]

			mov edx, hvn_a_reg_1
			movaps xmm6, [edx]

			mov edx, hvn_a_reg_2
			movaps xmm5, [edx]


		//блоки по 16 байт
		shorty0_iter:

			movaps xmm0, [esi]
			psubb xmm0, [edi]

			movaps xmm1, xmm0

			//get signs mask
			pcmpgtb xmm1, xmm7
			pxor xmm1, xmm7

			movaps xmm2, xmm1
			movaps xmm3, xmm1

			//load ones
			pand xmm2, xmm5

			//load 0x0f's
			pand xmm3, xmm6

			pslldq xmm2, 1

			paddb xmm0, xmm3
			psubb xmm0, xmm2

			movaps [esi], xmm0

			add edi, 0x10
			add esi, 0x10

		dec ecx
		jnz shorty0_iter

		shory1:

			mov ecx, eax
			test ecx, ecx
			jz exiorty

		//побайтово
		shory1_iter:

			mov al, [esi]
			sub al, [edi]
			mov [esi], al

			inc esi
			inc edi

		dec ecx
		jnz shory1_iter

		exiorty:

		mov esi, bu1
		mov ecx, oll

		mov edx, ecx
		sar ecx, 2

		test ecx, ecx
		jz carry_finish

		carry_iter:

			mov eax, [esi]
			and eax, 0x80808080
			test eax, eax
			jz skipolo2

				mov al, [esi + 0]
				test al, al
				jns sign_okay0
					add byte ptr [esi + 0], 0x10
					dec byte ptr [esi + 1]
				sign_okay0:

				mov al, [esi + 1]
				test al, al
				jns sign_okay1
					add byte ptr [esi + 1], 0x10
					dec byte ptr [esi + 2]
				sign_okay1:

				mov al, [esi + 2]
				test al, al
				jns sign_okay2
					add byte ptr [esi + 2], 0x10
					dec byte ptr [esi + 3]
				sign_okay2:

				mov al, [esi + 3]
				test al, al
				jns sign_okay3
					add byte ptr [esi + 3], 0x10
					dec byte ptr [esi + 4]
				sign_okay3:

			skipolo2:

			add esi, 4

		dec ecx
		jnz carry_iter

		/*carry_iter:

			//cmovs ecx, eax
			mov eax, [esi]
			and eax, ebx
			test eax, eax
			jnz fail4

				add esi, 4
				sub ecx, 4
				cmovs ecx, edx

			fail4:

			//lodsb
			mov al, [esi]
			test al, al
			jns sign_okay

				add byte ptr [esi], 0x10
				dec byte ptr [esi + 1]

			sign_okay:

			inc esi

		dec ecx
		jnz carry_iter   */

		carry_finish:

		mov ecx, edx
		and ecx, 3
		test ecx, ecx
		jz carry_false

		carry_byte:

			mov al, [esi]
			test al, al
			jns sign_okay

				add byte ptr [esi], 0x10
				dec byte ptr [esi + 1]

			sign_okay:

			inc esi

		dec ecx
		jnz carry_byte

		carry_false:

		movd edi, mm0
		movd esi, mm1

		emms
		//pop esi
		//pop edi
	}

	//последний цифр в этом случае никогда не будет отрицательным,
	//и поэтому баланс во вселенной не нарушиться
	/*for(i = 0; i < oll; i++){
		//это какбы значит что он меньше нуля
		while(this->nbuf[i] & 0x80){
			this->nbuf[i + 1]--;
			this->nbuf[i] += 16;
		}
	} */

	//меняем знак, если вдруг компейр был отрицательным
	if(comp < 0){
		this->sign = !this->sign;
	}

	//тримаем ненужные нули справа
	this->trim_right();

	return true;
}


#pragma pack(push, 1)
typedef struct {
	char bytez[16];
} TSSERegister;
#pragma pack(pop)

//структуры для реализации ултра-быстрого вычитания
static const TSSERegister 
SSEC_80x16x0 = {0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF,0xFF};

static const TSSERegister 
SSEC_0Fx16x1 = {0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x10,0x00};

static const TSSERegister 
SSEC_01x16x2 = {0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x01,0x00};

static const TSSERegister 
SSEC_SFx16x3 = {0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80,0x80};

void hvn_init(){
	hvn_a_regs = (TSSERegister *) malloc_a16(4 * 16);

	memcpy(&hvn_a_regs[0], &SSEC_80x16x0, 16);
	memcpy(&hvn_a_regs[1], &SSEC_0Fx16x1, 16);
	memcpy(&hvn_a_regs[2], &SSEC_01x16x2, 16);
	memcpy(&hvn_a_regs[3], &SSEC_SFx16x3, 16);

	hvn_a_reg_0 = hvn_a_regs + 0;
	hvn_a_reg_1 = hvn_a_regs + 1;
	hvn_a_reg_2 = hvn_a_regs + 2;
	hvn_a_reg_3 = hvn_a_regs + 3;
}

Код здесь, конечно, выглядит ещё более убийственно, но попробуем разобраться. В самом начале выбирается бо́льшее из чисел, чтобы из него начать вычитать меньшее. Затем, собственно, начинает производиться вычитание блоков по 16 байт с помощью PSUBB, но не простое вычитание. После проведения вычитания, надо полагать, что некоторые цифры, то-есть наши байты окажутся меньше нуля, для них, собственно, нужно выполнить поразрядный перенос, в обратную сторону. Для их поиска используется команда PCMPGTB, которая в нашем случае сравнивает все байты с байтом 0xFF, и в случае если число больше — в байт результата тот-же 0xFF и записывает. Затем байты инвертируются с помощью PXOR и 0xFF. В регистре XMM1 получается карта байтов, которые требуют переноса и коррекции. Что-же нужно сделать, для проведения коррекции? К переполненным байтам прибавить 0x10, а у левых их соседей отнять единичку. Именно это дальше и происходит. Из регистра XMM1 с помощью команды PAND создаются необходимые маски, маска с единичками сдвигается на 1 байт влево командой PSLLDQ, и, собственно “десятки” прибавляются, единицы вычитаются. Только после этого результат вычитания сохраняется.


Сразу надо сказать, что здесь есть несколько подвохов. Во первых, число может быть в общем-то и не кратно 16 — на этот случай за этим циклом следует побайтовый цикл, который занимается вычитанием как раз таких не кратных 16 хвостов. Но кроме того, нужно заметить, что во первых мы не имеем возможности сделать вычитание единички из соседнего блока (именно поэтому при проведении коррекции не обрабатывается крайний правый байт — в рабочих константах он обнулён), а во вторых, после вычитания единичек, левые соседи тоже могут переполниться, и чтобы наша коррекция стала полной, мы, по идеи, операцию коррекции должны повторить 15 раз. Но нам это делать лень. Вся хитрость использования этой коррекции заключается в том, чтобы значительно разгрузить алгоритм полной коррекции, который можно обнаружить в конце функции. Дело всё опять в вероятности. Ведь не трудно догадаться, что вероятность переполнения левого соседа составляет примерно 16−1 и в итоге получается, что не скорректированных символов остаётся примерно 1 на 16 байт. За счёт этого, появилась отличная возможность алгоритм финальной коррекции сделать с пропусками, чтобы он обрабатывал сразу 4 байта, и только в случае необходимости коррекции начинал побайтовую обработку. Эффективность такого подхода себя оправдала — после добавления такой коррекции к простому вычитанию через PSUBB, возведение в степень стало работать значительно шустрее, почти на секунду.


Подводим скромные итоги


Итак… В целом, за счёт всех представленных оптимизаций, удалось уменьшить время возведения в степень с 10 секунд, до одной секунды, а в более простом случае возведения числа с степень порядка 10100 по модулю 10300, время работы составляло рекордные 350 миллисекунд. Результат, конечно, неплохой. Но есть одно но. Дело в том, что возведение в степень по модулю является неотъемлемой частью такого алгоритма, как тест Миллера — Рабина, вероятностного теста чисел на простоту, без которого трудно обойтись работая с криптографией. Всё дело в том, что в этом тесте возведение в степень выполняется десятки, и даже сотни раз, причём число порядка 10300 здесь уже выступает показателем степени, что замедляет выполнение одного возведения до 1-2 секунд. А возведений десятки. Вот и получается, что в итоге, даже с использованием всех вышеописанных оптимизаций, ждать генерации случайного простого числа порядка 10300 всё равно приходиться несколько минут… Нужно согласиться, что для применения в реальных программах это неприемлемо. Ну согласитесь, кто будет ждать подключения программы к серверу 3 минуты, каждый раз?…


Решить эту проблему, мне честно говоря, не удалось. Возможно, здесь бы мог помочь алгоритм модульного умножения Монтгомери, или подобные, но терпения это проверять у меня не хватило. Для быстрого возведения в степень пришлось использовать GMP, в котором то-же возведение в степень работает за 20 миллисекунд. Каким образом это там достигается, разобраться тоже не особо получилось. Довольно печальный конец такой казалось-бы захватывающей истории… Впрочем, справедливости ради, наш оптимизированный алгоритм всё-же остался, но использовался там, где время было не так критично — в алгоритме Диффи — Хеллмана. Тест Миллера — Рабина же использовал GMP. Кстати, вот как он выглядел:


Код: C++
#include "gmp\mini-gmp.h"

//делявки на 2 и умножавки на 2
void __fastcall CHeavyNumber::div2n(long n){
	this->op_shift_right(n >> 2);
	long srn = n & 3, mxs = this->nlen - 1, i;
	for(i = 0; i < mxs; i++) this->nbuf[i] = this->nbuf[i] >> srn | (this->nbuf[i + 1] << (4 - srn) & 0xf);
	this->nbuf[mxs] >>= srn;
	this->trim_right();
}


void __fastcall CHeavyNumber::mul2n(long n){
	this->op_shift_left(n >> 2);
	long srn = n & 3, mxs = this->nlen, i;
	this->set_length(mxs + 1);
	for(i = mxs; i > 0; --i) this->nbuf[i] = this->nbuf[i] << srn & 0xf | (this->nbuf[i - 1] >> (4 - srn));
	this->nbuf[0] <<= srn;
	this->nbuf[0] &= 0xf;
	this->trim_right();
}

long __fastcall CHeavyNumber::is_even() const {
	//return (this->nbuf[0] & 1) == 0;
	long ris = 0;
	long i, j;
	for(i = 0; i < this->nlen; i++){
		if(0 == (0xf & this->nbuf[i])){
			ris += 4;
		}else{
			for(j = 0; j < 4; j++){
				if(0 == (this->nbuf[i] & (1 << j))){
					ris++;
				}else{
					return ris;
				}
			}
		}
	}

	return ris;
}

////////////////////////////////////////////////////////////////////////////////

void CHeavyNumber::qpowm_gmp(const CHeavyNumber &n, const CHeavyNumber &p){
	//mpz_t
	mpz_t dest, src, nxp, mod;
	mpz_init(dest);
	mpz_init(src);
	mpz_init(nxp);
	mpz_init(mod);

	long lep_src = this->nlen & 1 ? this->nlen + 1 >> 1 : this->nlen >> 1;
	long lep_nxp = n.nlen & 1 ? n.nlen + 1 >> 1 : n.nlen >> 1;
	long lep_mod = p.nlen & 1 ? p.nlen + 1 >> 1 : p.nlen >> 1;

	unsigned char *bu_src = new unsigned char[lep_src];
	unsigned char *bu_nxp = new unsigned char[lep_nxp];
	unsigned char *bu_mod = new unsigned char[lep_mod];

	this->pack(bu_src);
	n.pack(bu_nxp);
	p.pack(bu_mod);

	mpz_import(src, lep_src, -1, 1, 0, 0, bu_src);
	mpz_import(nxp, lep_nxp, -1, 1, 0, 0, bu_nxp);
	mpz_import(mod, lep_mod, -1, 1, 0, 0, bu_mod);

	delete bu_src;
	delete bu_nxp;
	delete bu_mod;

	mpz_powm(dest, src, nxp, mod);

	long lep_dest = mpz_sizeinbase(dest, 0x10);
	lep_dest = lep_dest & lep_dest ? lep_dest + 1 >> 1 : lep_dest >> 1;

	this->set_length(lep_dest << 1);
	unsigned char *bu_dest = new unsigned char[lep_dest];
	mpz_export(bu_dest, NULL, -1, 1, 0, 0, dest);
	this->unpack(bu_dest, lep_dest);
	delete bu_dest;

	mpz_clear(dest);
	mpz_clear(src);
	mpz_clear(nxp);
	mpz_clear(mod);
}

bool CHeavyNumber::miller_rabin(long r) const {
	if(this->is_zero()) return false;
	if(this->is_even()) return false;

	//small divisors test
	long k, si[10] = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31};

	CHeavyNumber tzeta(0);

	for(k = 0; k < 10; k++){
		tzeta = *this % si[k];
		if(tzeta.is_zero()) return false;
	}

	CHeavyNumber m1(*this);
	m1--;

	long s = m1.is_even();
	CHeavyNumber t(m1);
	t.div2n(s);
	CHeavyNumber a(0), x(0), one(1);

	bool is_fail;

	long i, j;
	for(i = 0; i < r; ++i){
		a.random_long(2, m1.nlen - 1);
		x = a;
		x.qpowm_gmp(t, *this);
		if(x == one || x == m1) continue;

        is_fail = true;

		for(j = 0; j < s; j++){
			x *= x;
			x %= *this;

			if(x == one) return false;
			if(x == m1){
				is_fail = false;
				break;
			}
		}

		if(is_fail) return false;
	}

	return true;
}

Впрочем, и для умножения на БПФ место здесь нашлось, и для быстрого взятия по модулю, но всё-равно, грустно как-то…


Мораль сей басни такова


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

Ключевые слова Программизм, математика, длинная арифметика, преобразование Фурье, многочлен, ДПФ, оптимизация, БПФ, SSE, BSF, PCMPGTB, тест Миллера — Рабина, алгоритмы
QR-код
Интересное
Итак, началось всё с чужого курсового по программированию, темой которого являлась как раз таки длинная арифметика на целых числах, а именно, нужно было всего-то реализовать 4 базовые арифметические операции: сложение, вычитание, умножение, деление.

Читать »»
Случайные фото