:: алгоритмы и методы :: :: олимпиадные задачи :: :: связь :: :: о сайте ::
Путь: Математика » Теория чисел » Факторизация » Методы Монте-Карло
  Методы Монте-Карло

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


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

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

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

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

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

b_0, b_1 = f(b_0), b_2 = f(b_1), b_3 = f(b_2), ... (2)

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

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

c_0 == b_0 (mod p), c_1 == b_1 (mod p), c_2 == b_2 (mod p), ... (3)

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

c_i = c_j, b_i != b_j.

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

b_i == b_j (mod p), b_i != b_j (mod m).

Отсюда вытекает, что (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).