Тест простоты Рабина

Тест простоты Рабина 
 © Борисенко, мех-мат МГУ
== Напомним необходимые нам результаты из элементарной теории чисел и алгебры ==
=== Малая теорема Ферма. ===
Пусть p - простое число. Тогда для всякого целого числа b, отличного от нуля, справедливо сравнение bp-1 = 1 (mod p).
Малая теорема Ферма является непосредственным следствием теоремы Лагранжа (порядок любого элемента группы делит порядок группы) и того факта, что кольцо Zp в случае простого p является полем, т.е. все его ненулевые элементы принадлежат группе обратимых элементов. Порядок группы обратимых элементов кольца Zp равен p-1.
Простейший тест проверки простоты числа m состоит в проверке малой теоремы Ферма. Выберем произвольное целое число b (например, b = 2), и возведем его в степень m - 1 по модулю m. Если мы получим не единицу, то по малой теореме Ферма число m составное. Беда состоит в том, что если
 bm-1 = 1 (mod m),
то ничего нельзя сказать об m. Древние греки ошибочно полагали, что все числа m, удовлетворяющие обращению малой теоремы Ферма для основания 2, простые: если
 2m-1 = 1 (mod m),
то m - простое число. Минимальный контрпример к этому утверждению был найден только в XVII веке:
 2340 = 1 (mod 341),
но число 341 - не простое, 341 = 11 * 31.
(Действительно, 2340 = (2^10)34 = 102434, но 1024 = 3 * 341 + 1 = 1 (mod 341), поэтому 102434 = 1 (mod 341).)
То, что 341 не удовлетворяет малой теореме Ферма, может быть показано с помощью других оснований:
 3340 = 56 (mod 341)
Тем не менее существуют числа, которые не являются простыми, но которые ведут себя как простые в малой теореме Ферма. Такие числа называются кармайкловыми.
=== Определение. ===
Число m называется кармайкловым, если оно не простое и для всякого b, взаимно простого с m, выполняется утверждение малой теоремы Ферма:
 bm-1 = 1 (mod m).
Минимальные кармайкловы числа -- это 561, 1105, 1729, ...
Множество кармайкловых чисел бесконечно, и их плотность стремится к нулю на бесконечности. Несложно доказать следующее утверждение.
Предложение 5. Пусть
 m = p1e1 * p2e2 * ... * pkek
это представление целого числа m в виде произведения степеней простых. Число m является кармайкловым тогда и только тогда, когда:
:1) для всякого i показатель степени ei = 1;
:2) k >= 3;
:3) для всякого i число pi - 1 делит m - 1.
=== Доказательство. ===
Докажем только обратную, наиболее интересную импликацию. Пусть число m удовлетворяет условиям 1-3.
Рассмотрим произвольное b, взаимно простое с m. По Китайской теореме об остатках, кольцо Zm представляется в виде прямой суммы
 Zm = Zp1 + Zp2 + ... + Zpk
При этом изоморфизме элемент b представляется в виде строки
 b = (b1, b2, ..., bk)
Тогда
 bm-1 = (b1m-1, b2m-1, ..., bkm-1
По малой теореме Ферма, для всякого i
 bim-1 = 1 (mod pi),
поскольку (m - 1) делится на (pi - 1).
Поэтому
 bm-1 = (1, 1, ..., 1)
т.е.
 bm-1 = 1 (mod m).
=== Пример. ===
Покажем, что число 561 является кармайкловым. Действительно, 561 = 3 * 11 * 17. Имеем
 (3 - 1) | 560, (11 - 1) | 560, (17 - 1) | 560.
Следовательно, число 561 удовлетворяет условиям предложения 5.
Итак, для кармайкловых чисел тест простоты, основанный на теореме Ферма, не работает. Тем не менее его модификация, предложенная Рабином, применима к любым целым числам.
Тест Рабина является вероятностным. Это означает, что он использует датчик случайных чисел и, таким образом, работает не детерминированно. Для входного целого числа m тест Рабина может выдать один из следующих двух ответов.
:1. Число m является составным.
:2. Не знаю.
В случае первого ответа число m действительно является составным, тест Рабина предъявляет доказательство этого факта. Второй ответ может быть выдан как для простого, так и для составного числа m. Однако для любого составного числа m вероятность второго ответа не превышает 1/4. Ценность теста Рабина состоит именно в неравенстве, ограничивающем сверху вероятность второго ответа для произвольного составного числа m.
Таким образом, если мы применим 100 раз тест Рабина к числу m и получим 100 ответов "не знаю", то можно с большой вероятностью утверждать, что число m простое. Более точно, вероятность получения ста ответов "не знаю" для составного числа m не превышает (1/4)100, т.е. практически равна нулю. Тем не менее тест Рабина не предъявляет доказательства того, что число m простое.
Перейдем непосредственно к изложению теста Рабина. Мы проверяем простоту входного числа m. Допустим сразу, что число m нечетное. (Существует только одно четное простое число - 2.) Тогда число m - 1 четное. Представим его в виде
 m - 1 = 2t * s
где s - нечетное число. Выберем случайное число b такое, что
 b ≠ 0
 b ≠ 1 (mod m)
 1 < b < m
При выборе b используется датчик случайных чисел.
Используя алгоритм быстрого возведения в степень по модулю m, вычислим следующую последовательность элементов кольца Zm:
 x0 = b2 (mod m), (1)
 x1 = x0 * x0 (mod m),
 x2 = x1 * x1 (mod m),
 ...
 xt = xt-1 * xt-1 = bm-1 (mod m)
(На каждом шаге мы возводим в квадрат число, полученное на предыдущем шаге.)
Тест Рабина выдает ответ 'm - составное число' в случае, если
:1) xt ≠ 1 (mod m), или
:2) в последовательности x0, x1, x2, ..., xt имеется фрагмент вида ..., *, 1, ... где звездочкой обозначено число, отличное от единицы или минус единицы по модулю m.
В противном случае тест Рабина выдает ответ "не знаю". Последовательность x0, x1, ..., xt в этом "плохом" случае либо начинается с единицы, либо содержит (-1) где-нибудь не в конце.
Cуществует алгоритм, доказывающий простоту, со сложностью O(ln3 n), согласно которому необходимо провести тест Рабина со всеми числами 
 2 <= b < k ln2 m
а затем проверить, не является ли m степенью простого числа. Однако его правильность зависит от недоказанной в настоящее время гипотезы Римана.
Этот алгоритм, опираясь на недоказанный факт, в принципе может 'соврать' в отношении доказательства простоты, хотя если тест Рабина говорит, что число составное, значит так оно и есть. На практике он работает очень даже неплохо.
== Теорема (законность теста Рабина). ==
:1. Если тест Рабина выдает ответ 'm - составное число', то m действительно является составным.
:2. Вероятность ответа 'не знаю' для составного числа m не превосходит 1/4.
=== Доказательство. ===
Докажем только первое утверждение. Если xt ≠ 1 (mod m), то m не удовлетворяет малой теореме Ферма и, следовательно, не является простым. Если же последовательность (1) содержит фрагмент ..., a, 1, ..., где a ≠ ±1 (mod m), то имеем
 a2 = 1 (mod m)
 a ≠ 1
 a ≠ -1 (mod m)
Если бы m было простым, то кольцо Zm являлось бы полем.
Но в любом поле есть только два квадратных корня из единицы: это единица и минус единица. (По теореме Безу, число корней многочлена не превосходит его степени, квадратные корни из единицы - это корни многочлена x2 - 1.) Следовательно, число m не является простым.
== Пример программы ==
{ IsPrime.Pas ver. 2.0 (c) Max Alekseyev <relf@os2.ru>, 2:5015/60@FidoNet }
{ Реализация вероятностного алгоритма Миллера-Рабина с 20 раундами.
 Для примера выдает простые на отрезке [1000000000, 1000100000].
 Вероятность ошибки (то, что составное число будет названо простым) меньше
 4^(-Rounds). }
const
 Rounds = 20;
function MulMod(x, y, m: LongInt): LongInt; assembler;
asm
{$IFDEF USE32}
 mov eax,x
 mul y
 div m
 mov eax,edx
{$ELSE}
 db $66; mov ax,word ptr x
 db $66; mul word ptr y
 db $66; div word ptr m
 mov ax,dx
 db $66; shr dx,16
{$ENDIF}
end;
function PowMod(x, a, m: LongInt): LongInt;
var
 r: LongInt;
begin
 r := 1;
 while a > 0 do
 begin
 if Odd(a) then
 r := MulMod(r, x, m);
 a := a shr 1;
 x := MulMod(x, x, m);
 end;
 Result := r;
end;
function IsPrime(p: LongInt): Boolean;
var
 q, i, a: LongInt;
begin
 if Odd(p) and (p > 1) then
 begin
 Result := True;
 q := p - 1;
 repeat
 q := q shr 1;
 until Odd(q);
 for i := 1 to Rounds do
 begin
{$IFDEF USE32}
 a := Random(p-2)+2;
{$ELSE}
 a := 2 + Trunc(Random * (p-2));
{$ENDIF}
 if PowMod(a, p-1, p) <> 1 then
 begin
 Result := False;
 Break;
 end;
 a := PowMod(a, q, p);
 if a <> 1 then
 begin
 while (a <> 1) and (a <> p-1) do
 a := MulMod(a, a, p);
 if a = 1 then
 begin
 Result := False;
 Break;
 end;
 end;
 end;
 end
 else
 Result := (p = 2);
end;
var
 t: LongInt;
begin
 Randomize; { Don't forget to reset Random Generator! }
 for t := 1000000000 to 1000100000 do
 if IsPrime(t) then WriteLn(t);
end.
Source: http://algolist.manual.ru
ID: 04240
Category:Алгоритмы
Category:Математика