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

Наиболее наивный подход - просто делить на все числа подряд до корня из N. Если Вы хотели узнать это - дальше можно не читать ;-)

     Пусть мы хотим поделить N на все простые числа до корня из N. Для этого мы можем иметь или не иметь в своем распоряжении достаточно большую таблицу простых чисел. Если это не так, то, очевидно, мы можем делить N на числа d в заданных классах эквивалентности, например, 1 и 5 по модулю 6, или 1,7,11,13,17,19,23,29 по модулю 30. Тогда мы проделаем много бессмысленных делений (на составные числа), однако результат останется верным. Таким образом, можно составить следующий алгоритм:


     Предположим, что у нас уже есть таблица простых чисел: p[1] = 2, p[2] = 3, ... , p[k], где k > 3, массив t := [6,4,2,4,2,4,6,2], и индекс j, такой что если p[k] mod 30 равно 1,7,11,13,17,19,23 или 29, то j := 0,1,2,3,4,5,6 или 7 соответственно.
     Выберем некоторую верхнюю грань B, такую что B >= p[k], чтобы не тратить на этот алгоритм слишком много времени.

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

1. [Инициализация]
        Если N <= 5, вывести разложение 1 = 1, 2 = 2, 3 = 3, 4 = 2², 5 = 5 в зависимости от N и выйти. Иначе, i := -1, m := 0, L := [ N^1/2 ]
   
   
2. [Cледующее простое]
        Пусть m := m + 1. Если m > k, то i : = j - 1 и идти к шагу 5, иначе d : = p[m].
   
   
3. [Простое деление]
        Пусть r : = N mod d. Если r = 0, то вывести d как нетривиальный делитель N и выйти ( или N := N / d, L : = [ N^1/2 ] и повторить шаг 3, если мы хотим продолжить нахождение делителей N).
   
   
4. [Простое?]
        Если d >= L, то в случае N > 1 вывести сообщение, что оставшаяся часть N - простое и выйти. Иначе, если i < 0 идти на шаг 2.
   
   
5. [Следующий делитель]
        Пусть i := i + 1 mod 8, d := d + t[i]. Если d > B, то вывести сообщение о том, что оставшиеся простые делители N больше B, иначе идти на шаг 3.

     Заметим, что i = -1 пока мы используем нашу таблицу простых и i >= 0 если нет. Этот тест не следует использовать для полной факторизации, кроме тех случаев, когда N - очень маленькое ( например, N < 10⁸), так как для этого есть лучшие методы. Однако, с другой стороны, он весьма полезен для удаления малых множителей.


Заметки к практической реализации.
     Предположим, что мы используем таблицу простых до 500'000 ( 41'538 простых чисел ). Простое деление при таком ограничении никогда не занимает больше нескольких секунд на современных компьютерах. Причем хранить можно только разности между простыми (или их половины), а не сами числа, так как p[k] - p[k-1] можно поместить в 1 байт вместо четырех, если p[k] <= 1'872'851'947 ( а половину этой разности - если p[k] <= 1'999'066'711'391).

     Кроме того, не нужно делать больше делений после того, как кончится таблица простых, так как есть лучшие методы для удаления малых делителей.

     И, наконец, заметим, что нет необходимости в вычислении L := [ N^1/2 ] во время инициализации, так как проверка d >= L в шаге 4 может быть заменена на q <= L, где q - евклидово частное при делении N на d : N = q * d + r, обычно вычисляемое одновременно с остатком в шаге 3.

Методы Монте-Карло

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

  Метод Монте-Карло 1

Время работы этого метода порядка корень квадратный из минимального простого числа, делящего m. То есть в худшем случае, когда m есть произведение двух простых чисел примерно одного порядка, число m раскладывается этим методом на множители за время корень четвертой степени из m. Алгоритм очень прост, его запись намного короче, чем объяснение, почему он работает. Слова "Монте-Карло" присутствуют в его названии потому, что работа алгоритма зависит от случайного выбора начального числа.

Рассмотрим отображение f кольца Z_m в Z_m:

f: Z_m --> Z_m
f(x) = x^2 + 1 (mod m)

Выберем случайным образом элемент b_0 кольца Z_m.
Рассмотрим бесконечную последовательность элементов кольца Z_m:

Последовательность представляет собой орбиту элемента b_0 при отображении f. Поскольку все элементы последовательности принадлежат конечному множеству, последовательность циклическая -- точнее, она содержит начальный апериодический отрезок и далее бесконечно повторяющийся период.

Пусть p -- делитель числа m. Рассмотрим элементы последовательности (2) по модулю p (т.е. образы элементов b_i при каноническом эпиморфизме Z_m --> Z_p):

Так как в Z_p меньше элементов, чем в Z_m, то с большой вероятностью период последовательности (3) меньше, чем период последовательности (2).
Следовательно, найдется пара индексов i, j такая, что

Это означает, что

Отсюда вытекает, что (b_i - b_j) делится на p, но не делится на m. Следовательно, НОД(b_i - b_j, m) нетривиален, и нам удалось разложить m на множители.

Итак, алгоритм Полларда 1 сводится к поиску цикла в бесконечной рекурсивной последовательности, состоящей из элементов конечного множества. При этом вместо того, чтобы сравнивать между собой два элемента, мы вычисляем наибольший общий делитель их разности и числа m. Алгоритм завершается, когда наибольший общий делитель нетривиален.

Можно предложить 2 способа решения задачи поиска цикла в последовательности. Первый способ наиболее простой. Второй чуть-чуть сложнее, но зато более быстрый.

Способ 1

---

Выполняется следующая последовательность сравнений:
 b0 <---> b1
 b1 <---> b3
 b2 <---> b5
 b3 <---> b7
 b4 <---> b9
 . . .
 b_i <---> b_{2*i+1}
 . . .
 
 

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

Выпишем алгоритм нахождения делителя.

алгоритм факторизация1(вход: целое число m, выход: целое число d): успех | дано: целое число m | надо: получить нетривиальный делитель d числа m | возвращаемое значение: true, если удалось разложить, | false в противном случае начало | maxSteps := 1000000 // Максимальное число шагов | step := 0 | | b0 := случайное число в интервале 0..m | b1 := mod(b0 * b0 + 1, m) | d := gcd(b1 - b0, m) | | цикл пока step < maxSteps && d == 1 выполнять // Пока НОД тривиален | | b0 = mod(b0 * b0 + 1, m) // Применяем отображение f | | b1 = mod(b1 * b1 + 1, m); // один раз к b0 и дважды | | b1 = mod(b1 * b1 + 1, m) // к b1 | | d := gcd(b1 - b0, m) | | step := step + 1 | конец_цикла | | вернуть (d != 1) // Успех := d нетривиален конец_алгоритма

На каждом шаге цикла мы трижды вычисляем значение отображения f.

Небольшая модификация алгоритма позволяет делать это только один раз.

Способ 2

---

Выполняется следующая бесконечная последовательность сравнений

b0 <---> b1

b1 <---> b2

b1 <---> b3  

b2 <---> b4

b2 <---> b5

b2 <---> b6

b2 <---> b7

b4 <---> b8

b4 <---> b9

b4 <---> b10

. . .  

b4 <---> b15

b8 <---> b16

b8 <---> b17

. . .

b8 <---> b31

. . .

Вся последовательность сравнений разбивается на серии. В очередной серии мы сравниваем элемент b_s, где s -- степень двойки, последовательно с элементами b_{2s}, b_{2s+1}, b_{2s+2}, ..., b_{4s-1}. Серия содержит 2s сравнений.

  Выпишем алгоритм.

алгоритм факторизация2(вход: целое число m, выход: целое число d): успех | дано: целое число m | надо: получить нетривиальный делитель d числа m | возвращаемое значение: true, если удалось разложить, | false в противном случае начало | maxSteps := 19 // Максимальная длина серии 2^19 | step := 0 | | b0 := случайное число в интервале 0..m | b1 := mod(b0 * b0 + 1, m) | a := b1 // Первый элемент серии | seriesLength := 1 // Длина серии | | цикл пока step < maxSteps && d == 1 выполнять // пока НОД тривиален | | Инвариант: b0 -- элемент последовательности с индексом, | | равным нулю или степени двойки | | a -- элемент, индекс которого равен удвоенному индексу | | элемента b0 (или 1, если индекс b0 равен 0) | | seriesLength == удвоенному индексу элемента a | | d := gcd(b1 - b0, m) | | len := 0 | | | | цикл пока d == 1 и len < seriesLength выполнять | | | b1 = mod(b1 * b1 + 1, m); | | | d := gcd(b1 - b0, m) | | | len := len + 1 | | конец_цикла | | | | b0 := a | | a := b1 | | seriesLength := seriesLength * 2 | конец_цикла | | вернуть (d != 1) // Успех := d нетривиален конец_алгоритма    

 Метод Монте-Карло 2

---

Пусть m --- целое число, которое мы раскладываем на множители. Оно представимо в виде произведения степеней простых чисел

m = p1^e1 p2^e2 ... pk^ek

Предположим, что p1 - 1 представимо в виде произведения степеней простых чисел, причем каждая из этих степеней не очень велика. Более точно, существует N такое, что

p1 - 1 = q1^a1 q1^a2 ... qr^ar q1^a1 < N, q2^a2 < N, ..., qr^ar < N.

Рассмотрим всевозможные максимальные степени простых чисел, не превосходящие N. Например, пусть N = 20, тогда рассматриваются степени простых 16, 9, 5, 7, 11, 13, 17, 19. Обозначим эти степени простых через t1, t2, ..., ts.

Выберем произвольное целое число b. Рассмотрим последовательность

b0 = b, b1 = b0^t1 (mod m), b2 = b1^t2 (mod m), ... bs = b_{s-1}^ts (mod m)

Каждый раз, вычислив bi, вычисляем одновременно

НОД(bi - 1, m).

Утверждается, что с большой вероятностью на каком-то шаге этот НОД будет нетривиальным делителем N. Действительно, покажем, что

p | НОД(bs - 1, m).

Действительно,

bs = b^(t1 t2 ... ts),

и, поскольку по предположению, p1 - 1 | t1 t2 ... t3, то есть t1 t2 ... ts = (p1 - 1)g, то

bs = b^(t1 t2 ... ts) = b^((p1 - 1) g) = (b^(p1 - 1))^g = 1 (mod p1)

по малой теореме Ферма. Значит, bs - 1 делится на p1, число m также делится на p1, следовательно, НОД(bs - 1, m) делится на p1.

Проиллюстрируем алгоритм на простом примере. Возьмем N = 20. Выпишем все степени простых, не превосходящие 20:

t1 = 16, t2 = 9, t3 = 5, t4 = 7,

t5 = 11, t6 = 13, t7 = 17, t8 = 19.

Попытаемся разложить на множители число m = 41779 = 41 * 1019. Выберем b = 2. Последовательно вычисляем

2 ^ 16 (mod 41779) = 23757, gcd(23757 - 1, 41779) = 1, 23757 ^ 9 (mod 41779) = 7970, gcd(7970 - 1, 41779) = 1, 7970 ^ 5 (mod 41779) = 33580, gcd(33580 - 1, 41779) = 41.

Мы получили нетривиальный делитель на третьем шаге, поскольку 41 - 1 = 8 * 5 делит (t1 * t2 * t3) = 16 * 9 * 5.

Мощность алгоритма зависит от числа N --- чем больше оно, тем большие числа можно разложить с помощью этого алгоритма. Работа алгоритма разбивается на 2 шага. Сначала мы генерируем все максимальные степени простых чисел, не превосходящие N. Этот шаг выполняется только один раз и не зависит от входного числа m, поэтому сгенерированные степени можно, к примеру, записать в файл и в дальнейшем использовать многократно. Затем мы выбираем случайным образом число b и вычисляем указанную выше последовательность степеней b. Для каждой степени bi вычисляется НОД(bi - 1, m).

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

(b_i - 1) (b_{i+1} - 1) (b_{i+2} - 1) ... (b_{i+99} - 1) = n (mod m)

и затем вычисляется НОД(n, m).