DIGIT Аdd(

DIGIT C[ ], // результат

const DIGIT A[ ], // первое слагаемое

const DIGIT B[ ], // второе слагаемое

int n) // длина слагаемых

{              TWODIGIT T;

DIGIT d=0;

int i;

for(i=0; i<n; i++)

{              T = (TWODIGIT)A[i]+B[i]+d;

C[i] = LODIGIT(T);

d = HIDIGIT(T);

}

return d;

}

Сложность данного алгоритма O(n). Реализация этого алгоритма проста, так как он имитирует известную всем процедуру сложения в «столбик».

С умножением дело обстоит не так просто, как с операциями сложения и вычитания. Дело в том, что известный нам со школы алгоритм умножения “в столбик” не самый быстрый. Вот он.

Алгоритм 1.3. Умножение «в столбик» длинных чисел С = А*В

длины n цифр

1.       C := 0; (достаточно обнулить n младших цифр)

2.       для i = 0 … n-1 выполнить шаги 3-8;

3.       d :=0;

4.       для j = 0 … n-1 выполнить шаги 3-5;

5.       T :=C_i+j + A_i*B_i +d;

6.       C_i+j : = LODIGIT(T);

7.       d := HIDIGIT(T);

8.       C_i+n := d;

9.       конец.

Обратите внимание на то, что результат операции умножения в общем случае в 2 раза длиннее сомножителей, то есть имеет длину 2n цифр. Нетрудно видеть, что описанный алгоритм требует выполнения n операций умножения и еще некоторого количества операций сложения. И именно так умножали люди с давних времен, пока 39 лет назад московским математиком А.А. Карацубой не было совершено неожиданное открытие. Он придумал куда более эффективный метод.

Метод Карацубы

Пусть А и В — два n-значных длинных числа и n – четное. Эти числа можно представить в виде

А = b^n/2U₁ + U₀, B = b^n/2V₁ + V₀,

где U₀, U₁,V₀,V₁ — n-значные числа. Обычное умножение “в столбик” равносильно следующему:

АЧВ = (b^n/2U₁ + U₀)(b^n/2V₁ + V₀ ) = bⁿU₁V₁+ b^n/2 (U₁V₀ + U₀V₁) + U₀V₀,

те. требует выполнения 4 операций умножения -значных чисел.

Однако, воспользовавшись тождеством

(U₁ - U₀) (V₀ -V₁) = -U₁V₁-U₀V₀+ U₁V₀+U₀V₁,

получаем

АЧВ=(b^n/2U₁+U₀)(b^n/2V₁+V₀ )=(bⁿ+ b^n/2)U₁V₁+ b^n/2(U₁ - U₀)(V₀ -V₁)+(b^n/2+1)U₀V₀.

При этом умножение n-значных чисел потребует выполнения З операций умножения-значных чисел и еще некоторое количество операций сложения и вычитания. Итак, налицо выигрыш  по сравнению с умножением в столбик. Но поскольку и умножение -значных чисел тоже можно выполнять таким же образом, организован рекурсивный вызов процедуры умножения и т.д., то теоретически получим сложность около n^1.6 элементарных умножений, а не n², как раньше.

На практике не все столь радужно, как показалось. При малых значениях n скорее всего можно проиграть на накладных расходах (организация рекурсии и циклов, дополнительные операции вычитания, трудности с учетом знака результата и т.п.). Практические вычисления показали, что реальный выигрыш наступает при разрядности операндов где-то около 16384 битов, что не может интересовать, поскольку ГОСТ Р 34.10-94 регламентирует длину модуля не более 1024, да и в реальной практике пока модулярные вычисления с разрядностью выше 2048 битов используются редко. А по сему реализацию этого метода рассматривать здесь не будем.

Интересующимся же заметим, что метод Карацубы был в 1963 году обобщен А. Тоомом, который нашел связь этой задачи с интерполяционными полиномами. И еще более продвинулись в решении задачи умножения А. Шенхаге и Ф. Штрассен, которые помимо интерполяционных полиномов использовали быстрое преобразование Фурье и построили алгоритм умножения с числом элементарных операций умножения не более О(nЧlog nЧlog log n). Однако, как и метод Карацубы, эти алгоритмы слишком сложны в реализации и дают практический выигрыш только при очень большой разрядности сомножителей.

Умножение на цифру

Далее при делении понадобится функция умножения «длинного» n-значного числа на цифру. В принципе, можно поступить просто и обратиться к алгоритму умножения, предварительно обобщив его на операнды различной длины. Поступим по-другому

Алгоритм 1.4. Умножения числа А длины n цифр на цифру х:

С = А*х

1.       Т:=0;

2.       для i = 0.. .n-1 выполнить 3-4;

3.       Т:=Т/b+А_i Ч х;

4.       С_i := LODIGIT(Т);

5.       Сn := HIDIGIT(T);

б.      конец.

Кстати, при этом получается число С длины n+ 1. Даже, если старшая цифра результата равна нулю, требуем, чтобы для нее была выделена память.

При реализации алгоритма Рабина-Миллера был создан класс длинной арифметики long_ar, где в качестве функций были встроены операции сложения длинных чисел, вычитание длинных чисел, умножение длинных чисел на число, перемножение двух длинных чисел, поиск остатка от деления двух длинных чисел. Также в класс long_ar вошли такие функции, как обнуление длинных чисел, печать длинных чисел, создание длинных чисел при чтении из файла, определение количества битов в длинном числе, проверка делимости длинного числа на три и на два, сравнение с нулем, единицей и с 100000, проверка нечетности, побитовый сдвиг вправо с присвоением и др.

Программы состоит из двух подпрограмм: is_probably_prime() и is_sprp(int n, int b) – и начинается с вызова is_probably_prime().

Подпрограмма is_probably_prime()

Шаг 1: ввести проверяемое число num.

Шаг 2: если num – четное, то вывести сообщение, что число составное, т.к. единственное простое четное число – это 2, а программа рассчитана на проверку длинных чисел. Завершить программу.

Шаг 3: если число num делится на три – вывести сообщение, что число делится на три (неделимость на три – необходимое условие для теста Рабина-Миллера). Завершить программу.

Шаг 4: ввести i > 0 (количество прохождений теста Рабина-Миллера). Если тест даст результат, что num – простое, то от величины i будет зависеть погрешность полученного результата.

Шаг 5: сгенерировать случайное нечетное число base, 1 < base < num.

Шаг 6: запустить подпрограмму is_sprp со входными параметрами (num, base).

Шаг 7: если is_sprp возвращает значение false – вывести сообщение, что число составное и завершить программу. Иначе уменьшить i на единицу.

Шаг 8: если i больше нуля перейти к шагу 5.

Шаг 9: Вывести сообщение, что число простое и завершить программу.

Подпрограмма is_sprp(int n, int b)

Шаг 1: num = n.

Шаг 2: s = 1. Тип s – integer.

Шаг 3: e = num >> 1. e – вспомогательная переменная. Операция «>>» обозначает побитовый сдвиг вправо.

Шаг 4: пока e четно, выполнять побитовый сдвиг e с присвоением   (e >>= 1) и инкапсуляцию s (s++). Как только e станет нечетно – переход к шагу 5.

Шаг 5: factor = 1, congr = b (и factor и congr – длинные числа).

Шаг 6: если e меньше двух – переход к шагу 11.

Шаг 7: если e нечетно, factor = (factor * congr)%num.

Шаг 8: congr = (congr * congr)%num.

Шаг 9: e >>= 1 (побитовый сдвиг с присвоением).

Шаг 10: переход к шагу 6.

Шаг 11: congr = (congr * factor)%num.

Шаг 12: если congr == 1, вернуть значение true.

Шаг 13: r = 0.

Шаг 14: если выполняется условие (congr == num – 1), вернуть значение true.

Шаг 15: congr = (congr * congr)%num.

Шаг 16: увеличить r на единицу.

Шаг 17: если r < s перейти к шагу 14.

Шаг 18: вернуть значение false.

• Если n < 2·10¹² и является 2, 3, 5, 7 и 11-SPRP, то оно простое.
• Если n < 3·10¹⁴ и является 2, 3, 5, 7, 11 и 13-SPRP, то оно простое.
• Если n < 3,4·10¹⁴ и является 2, 3, 5, 7, 11, 13 и 17-SPRP, то оно простое.

Следует обратить внимание на то, что эти тесты не вероятностные, они доказывают простоту числа. Логичным завершением рассмотрения такого подхода является утверждение, предложенное Миллером:

Тест Миллера: если верна расширенная гипотеза Римана, то если n является a-SPRP для всех a: 1 < a < 2 (log n)², то n – простое.

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

Если рассматривать по отдельности, то все приведенные ранее проверки на простоту работают либо недостаточно надежно, либо недостаточно быстро. Вместе с тем, если объединить тест Рабина-Миллера с пробным делением, мы получаем алгоритм, работающий намного лучше и того и другого. Дело в том, что для большого n вычислительно проще провести пробное деление на небольшое простое число, а не тест Рабина-Миллера. При этом сразу же отбрасывается достаточно большая часть составных чисел. К примеру, при проверке делимости на 3, 5 и 7 отбрасывается ~50% всех составных нечетных чисел, при проверке делимости на все простые меньшие 256 отбрасывается уже ~80% всех составных нечетных чисел.

Объединенный алгоритм (для k-битного числа):

1. Пробное деление на все числа меньше некоторого граничного числа B. Это число определяется временем полного возведения в степень по модулю k-битного числа и полным временем проверки делимости k-битного числа на небольшое простое число и равно примерно их отношению.

2. Тест Рабина-Миллера, проводимый t раз, где t выбирается из соображений необходимой точности.

Можно подобрать такие значения B и t можно получить существенно лучшие результаты, чем при использовании просто теста Рабина-Миллера. К примеру, для k = 500 и t = 6 вероятность ошибки такого алгоритма ~2^-92, что существенно меньше, чем 2^-2t = 2^-12. Объединенный алгоритм (тест Рабина-Миллера + пробное деление) находит широкое применение в криптосистемах с открытым ключом для построения простых ключей длиной 512, 1024 и 2048 бит.

Метод пробных делений

Этот метод относится к группе детерминированных методов. Пусть n Є N. Как проверить, является ли n простым? Пока не существовало необходимости генерировать большие простые числа, можно было использовать методы проверки, которые достаточно легко реализуемы без применения вычислительной техники и не требуют больших усилий для проверки маленьких чисел. Первым из таких методов является, естественно, полный перебор всех возможных делителей. Чаще всего используют модификацию такого перебора, называемую пробным делением (trial division): чтобы проверить число на простоту, делим его на все простые меньше либо равные корню из этого числа. Если n—составное, то n =ab, где , причем . Поэтому для d = 2, 3, . . ., [] следует проверить, делится ли n на d? Если делитель числа n не будет найден, то n—простое. В противном случае будет найден минимальный простой делитель числа n, т.е. мы даже разложим n на два множителя. Сложность метода составляет C*n^1/2 арифметических операций с целыми числами.

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

Возможны модификации этого метода, которые работают быстрее. Например, мы можем проверить, делится ли n на 2 и на 3 , и если нет, то перебираем далее только числа d вида 1 + 6j и 5+ 6j, j =1, 2, … Сложность метода отличается от предыдущей лишь постоянной в C₁(·)

Решето Эратосфена

Для построения простых чисел, меньших какому-либо N, можно воспользоваться так называемым решетом Эратосфена. Если мы хотим составить таблицу всех простых чисел среди чисел 2, 3, . . . , N, то нужно сперва вычеркнуть все числа, делящиеся на 2, кроме 2. Затем взять число 3 и вычеркнуть все последующие числа, делящиеся на3 . Затем взять следующее не вычеркнутое число (т. е. 5) и вычеркнуть все последующие делящиеся на него числа, и так далее. В итоге останутся лишь простые числа. Для реализации метода нужен большой объем памяти ЭВМ, однако для составления таблиц простых чисел он является наилучшим. В книге  описаны эффективные алгоритмы, реализующие решето Эратосфена для построения таблиц простых чисел и вычисление факторных баз.  Этот метод также как и метод пробных делений является детерминированным.

Вероятностные тесты

Пусть n Є N, n нечетно, n > 1. Вероятностный тест на простоту проводится следующим образом. Выбирается случайное a Є N, 1 ≤ a <n, и для него проверяется выполнение некоторых условий. Если какое-то из условий не выполнено, то число n—составное, поскольку для простых чисел эти условия являются необходимыми. Если же все условия выполнены, то из этого еще не следует простота n. Однако можно будет считать, что «n —простое число с некоторой вероятностью». Кроме того, обычно доказывают оценку снизу для этой вероятности. Чем больше значений a мы протестируем, тем ближе эта вероятность к единице.

Малая теорема Ферма

Для проверки простоты чисел n порядка 10³⁰—10⁴⁰ метод пробных делений уже неприменим, потому что даже на самых мощных компьютерах вычисления займут много времени (несколько лет). В 17 веке французский математик Пьер Ферма выдвинул утверждение, которое лежит в основе практически всех возможных тестов на простоту. Малая теорема Ферма: Если n – простое и a – любое целое, то . В частности, если n не делит a, то .

На основании этой теоремы можно построить эффективный тест на простоту.

Тест Ферма: для n > 1 выбираем a > 1 и проверяем малую теорему Ферма. Если малая теорема Ферма не выполнена, то n – составное. Если же она выполнена, то мы еще не можем сделать вывод о простоте n, поскольку теорема дает лишь необходимое условие. Этот тест является эффективным для обнаружения составных чисел.

Тест Рабина-Миллера и сильно возможно простые числа.

Можно существенно улучшить тест Ферма, заметив, что если n – простое нечетное, то для 1 есть только два квадратных корня по модулю n: 1 и -1. Таким образом, квадратный корень из a^n-1, a^(n-1)/2 равен плюс или минус единице. Если (n – 1)/2 опять нечетно, то мы сможем снова извлечь корень и так далее. Первый вариант алгоритма, предлагает использовать только одно деление:

Тест Леманна (Lehmann): если для какого-либо целого числа a меньшего n не выполняется условие a^(n-1)/2 ±1 (mod n), то число n – составное. Если это условие выполняется, то число n – возможно простое, причем вероятность ошибки не превышает 50%. Тест Леманна не используется, так как он хуже аналогичного вероятностного теста Рабина-Миллера. Но этот тест можно естественным образом улучшить, если извлекать корень по модулю не один раз, а столько, сколько получится.

Самым быстрым оказался алгоритм построенный на методе Брента. Он смог работать  более быстро с большими числами, чем остальные алгоритмы. Классический метод Полларда перестал работать с числами порядка 10³¹. Метод ускорения с помощью китайской теоремы об остатках работает с большими числами, но медленнее чем метод ускорения методом Брента.

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

2.1. [Проверить, будет ли число простым.] Если n – простое число

2.2. Вывести n в качестве результата; на этом выполнение алгоритма завершается.

3.1. [Множитель найден?] Присвоить g←gcd(x’-x, n).

3.2. Если g=1, то перейти к шагу 4;

3.3. В противном случае вывести g.

3.4. Если g=n, алгоритм завершается (его выполнение прерывается, ибо нам известно, что n не является простым числом).

3.5. В противном случае присвоить n ← n/g, x ← x mod n, x’ ← x’ mod n и возвратится к шагу 2. (Заметим, что g может и не быть простым числом – это подлежит проверке. В каждом случае, когда g – не простое число, его простые множители не могут быть определены при помощи этого алгоритма.)

4.1. [Продвинуться.] Присвоить k ← k – 1.

4.2. Если k=0,

4.3. То присвоить x’ ← x, l ← 2l, k ← l. Присвоить x ← (x²+1) mod n и возвратится к шагу 3.

Более того, в случае, когда операция нахождения наибольшего общего делителя выполняется медленно, Поллард предложил ускорить процесс путем накапливания произведения по модулю n, скажем, десяти последовательных значений (х’ - х) перед тем, как искать наибольший общий делитель. Таким образом, приблизительно 90% операций нахождения наибольшего общего делителя заменяется одним умножением по моду­лю N ценой некоторого увеличения вероятности того, что решение при этом не будет найдено.

В тех редких случаях, когда для больших N результат не был найден, можно применить функцию f(x) = х² +с при некотором с ≠ 0 или 1,  Следует избегать также значения с = -2, поскольку рекуррентное уравнение x_m+1 = x_m² – 2 имеет решение в виде . Похоже, что другие значения параметра с не приводят к возникновению простых связей по модулю р и все они должны быть удовлетворительными при подходящих начальных значениях.

Ускорение с помощью китайской теоремы об остатках

1. (Инициализация) Вводим x₀, k=0, y = 0;

2. Цикл от 1 до 10;

2.1 k = k + 1;

2.2 x_k = P(x_k-1);

2.3. Если k – нечетное, идти в 2.1;

2.4. j = k/2;

2.5. m_i = |x_k – x_j|

2.6. y = y + m_i;

3 Найти НОД (n, y);

4. Если НОД = 1, идти в 2.1;

5.1 Если НОД < n,

5.2 То p = НОД, конец;

5.3 Вывод p;

6. (Если НОД делится не только на p, но и на n). Идти в 1.

1. (Инициализация) Вводим x₀, k=0;

2. k = k + 1;

3. x_k = P(x_k-1);

4. Если k – нечетное, идти в 2;

5.1 j = k/2.

5.2 Найти НОД (n, |x_k – x_j|);

6. Если НОД = 1, идти в 2;

7.1 Если НОД < n,

7.2 То p = НОД, конец;

7.3 Вывод p;

8. Если НОД делится не только на p, но и на n. Идти в 1.

Доказательство Пусть числа и определены согласно m1m2 × × × mk = Msms, причем  º 1(mod ms),  s = 1, 2, . . . , k и целые числа  попарно взаимнопростые, то есть ,

Пусть b₁, b₂, . . . , b_k тоже целые числа такие что ,

Тогда система сравнений:

x º b1(mod m1),

x º b2(mod m2),

. . . . . . . . . . .

x º bk(mod mk)

Имеет в интервале где M = m₁×m₂ × × × m_k одно общее решение

(7)

Если числа  и определены из условий m1m2 × × × mk = Msms, причем  º 1(mod ms),  s = 1, 2, . . . , k и целые числа  попарно взаимнопростые, то есть , то совокупность значений x которые удовлетворяют систему (1) определяется сравнением  где  +, х – в виде уравнения (7) удовлетворяет систему (6). Согласно другой теореме – если b1, b2, . . . , bk независимо друг от друга пробегают полную систему остатков по модулям . Число x₀ может принять значение одного из остатков по модулю M = m₁×m₂ × × × m_k., т.е. .

Докажем что для x₀ из интервала для разных наборов целых чисел b₁, b₂, . . . , b_k таких что , , общее решение (7) является единственным. Доказательство будет вестись от обратного допустим что существует еще один набор чисел  таких, что ,  Пускай для конкретности эти наборы отличаются числами что стоят на s-позиции т.е. . Но тогда мы приходим к тому что справедливы сравнения  и  так как по условию теоремы при всех  выходит . Поэтому из равенства чисел вытекает справедливость сравнения . При условии  это значит  а это противоречит начальному условию. Следовательно решение (6) единственно.

Второй метод ускорения строится на алгоритме нахождения периода последовательности методом Брента. Пусть F – конечное множество, x F и f: F → F отображение, которому соответствует последовательность x_n+1 = f(x_n), начинающая с x₀.Цели алгоритма Брента – сосчитать длину периода последовательности (x_n)_n≥0, использую при этом как можно меньше памяти . Сложность алгоритма имеет порядок квадратичного корня из F.

Отметим, что для h = 2ⁱ – 1 интервал соответствующих значений k равен [h+1, 2h +1]. Для того, чтобы теперь x_h = x_k, необходимо и достаточно, чтобы h ≥ μ (индекс вхождения в период) и k – h было кратно λ (длина периода) – условие, которое выполняется при достаточно большом h (так как интервал для k удваивается на каждом этапе). Пусть h – какой-либо индекс, для каждого существует k  [h+1,2h+1] c x_h = x_k. Если k – первый индекс интервала [h+1,2h+1], для которого x_h = x_k, то k – h равно длине периода λ. Наконец, если μ – наименьший индекс, такой, что x_μ = x_μ+λ, то μ – индекс вхождения в период. Это приводит к алгоритму, в котором переменная j  принимает значения h+1, т.е. степени 2. Хотя алгоритм находит только длину периода, дальнейшее применение метода Брента позволяет определить все 4 параметра (h,k,μ,λ), описанные выше: сложность измеряется при помощи параметра k, который равен числу итераций x ← f(x) и обозначается k_f,x0.