Algorytmy i Struktury Danych - Test Millera-Rabina

Opisany w tym rozdziale algorytm testowania pierwszości (ang. primality testing algorithm) liczby naturalnej p został opracowany w roku 1975 przez Michaela O. Rabina na podstawie prac Gary'ego L. Millera. Udostępnia on szybką metodę sprawdzania pierwszości liczby z możliwością kontrolowania poziomu prawdopodobieństwa popełnienia błędu – jest to zatem metoda probabilistyczna, zwana testem Millera-Rabina (ang. the Rabin-Miller Probabilistic Primalty Algorithm).

Test Millera-Rabina oparty jest na następującym twierdzeniu:

Niech p będzie nieparzystą liczbą pierwszą zapisaną jako p = 1+2^sd, gdzie d jest nieparzyste. Wtedy dla dowolnej liczby naturalnej a ∈ <2, p-2> ciąg Millera-Rabina:

a^d, a^2d, a^4d, …, a^{2^(r-1)}, a^{2^r}(mod p)

kończy się liczbą 1. Co więcej, jeśli a^d nie przystaje modulo p do 1, to wyraz ciągu Millera-Rabina bezpośrednio poprzedzający 1 jest równy p-1.

Jeśli liczba p przejdzie test, to jest albo pierwsza, albo silnie pseudopierwsza przy podstawie a. Test Millera-Rabina daje złe wyniki (p złożona) dla co najwyżej ¹/₄ baz a < p. Zatem dla jednego przebiegu prawdopodobieństwo błędu wynosi ¹/₄. Dla n przebiegów przy różnych podstawach a prawdopodobieństwo błędu spada do (¹/₄)ⁿ. Już dwadzieścia testów daje szansę 1 do 1099511627776, iż liczba złożona p zostanie potraktowana jako pierwsza. Widzimy zatem wyraźnie, iż liczbę p możemy sprawdzić na pierwszość z bardzo dużym prawdopodobieństwem wykonując odpowiednią liczbę testów Millera-Rabina.

W algorytmie testu Millera-Rabina wykorzystujemy procedury mnożenia i potęgowania modulo, które opisaliśmy w rozdziale o chińskim teście pierwszości. Test ten wykorzystują obecnie prawie wszystkie systemy kryptografii publicznej do testowania pierwszości dużych liczb potrzebnych przy generacji kluczy szyfrujących/deszyfrujących.

Algorytm sprawdzania pierwszości nieparzystej liczby naturalnej testem Millera-Rabina

Wejście:

p : liczba badana na pierwszość, p ∈ N, p > 2, p jest nieparzyste.
n : ilość powtórzeń testu Millera-Rabina, n ∈ N.

Wyjście:

TAK, jeśli p jest pierwsze lub silnie pseudopierwsze z prawdopodobieństwem (¹/₄)ⁿ.
NIE , jeśli p jest liczbą złożoną.

Elementy pomocnicze:

s : wykładnik potęgi 2 w dzielniku p-1. s ∈ N.
d : mnożnik potęgi 2 w dzielniku p-1. d ∈ N.
i : zlicza wykonane testy Millera-Rabina, i ∈ N.
a : baza. a ∈ N. a ∈ <2, p-2>.
x : wyraz ciągu Millera-Rabina, x ∈ N.
j : zlicza wyrazy ciągu Millera-Rabina, j ∈ N.

Lista kroków:

K01: d ← p-1 ; obliczamy s i d
K02: s ← 0
K03: Dopóki d mod 2 = 0: ; usuwamy z p-1 dzielniki 2, 
     wykonuj kroki K04…K05 ; zliczając je w s
K04:   s ← s+1
K05:   d ← d div 2
K06: Dla i = 1, 2, …, n: ; wykonujemy n testów Millera-Rabina
     wykonuj kroki K07…K15
K07:   a ← Losuj(2, p-2) ; losujemy bazę a
K08:   x ← a^d mod p ; wyliczamy pierwszy wyraz ciągu Millera-Rabina
K09:   Jeśli (x = 1) ∨ (x = p-1), ; jeśli x nie spełnia warunku, 
       to następny obieg pętli K06 ; wybieramy inne a
K10:   j ← 1
K11:   Dopóki (j < s) ∧ (x ≠ p-1): ; rozpoczynamy generację kolejnych
       wykonuj kroki K12…K14       ; wyrazów ciągu Millera-Rabina
K12:   x ← x² mod p ; obliczamy kolejny wyraz
K13:   Jeśli x = 1, ; tylko ostatni wyraz ciągu Millera-Rabina
       to idź do kroku K17 ; może mieć wartość 1!
K14:   j ← j+1
K15:   Jeśli x ≠ p-1, ; przedostatni wyraz ciągu Millera-Rabina
       to idź do kroku K17 ; musi być równy p-1
K16: Pisz "TAK" i zakończ  ; pętla wykonała n testów
                           ; i zakończyła się naturalnie.
K17: Pisz "NIE" i zakończ  ; liczba p nie przeszła
                           ; testów Millera-Rabina

Przykładowe programy

Uwaga:

Zanim uruchomisz program, przeczytaj wstęp do tego artykułu, w którym wyjaśniamy funkcje tych programów oraz sposób korzystania z nich.

Program odczytuje z pierwszego wiersza liczbę p oraz ilość testów Millera-Rabina n, a w drugim wierszu wypisuje:

TAK – jeśli liczba p jest pierwsza z prawdopodobieństwem 1-(¹/₄)ⁿNIE – jeśli liczba p jest liczbą złożoną

Pascal

// Test Millera-Rabina
// Data: 6.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

program prg;

// 64 bitowy generator pseudolosowy
//---------------------------------
function Losuj(a, b : qword) : qword;

var w : qword;
    i : integer;
begin
  for i := 1 to 8 do
  begin
    w := w shl 8;
    w := w or random(256);
  end;
  Losuj := a+(w mod (b-a));
end;

// Funkcja mnoży a i b mod n
//--------------------------
function MnozModulo(a, b, n : qword) : qword;

var m, w : qword;

begin
  w := 0; m := 1;
  while m <> 0 do
  begin
    if b and m <> 0 then
      w := (w+a) mod n;
    a := (a shl 1) mod n;
    m := m shl 1;
  end;
  MnozModulo := w;
end;

// Funkcja oblicza a^e mod n
//--------------------------
function PotegujModulo(a, e, n : qword) : qword;

var m, p, w : qword;

begin
  p := a; w := 1; m := 1;
  while m <> 0 do
  begin
    if e and m <> 0 then
      w := MnozModulo(w, p, n);
    p := MnozModulo(p, p, n);
    m := m shl 1;
  end;
  PotegujModulo := w;
end;

var a, d, p, x : qword;
    i, j, s, n : integer;
    t : boolean;

begin
  randomize;
  readln(p, n);
  d := p-1;
  s := 0;
  while d mod 2 = 0 do
  begin
    inc (s); d := d div 2;
  end;
  t := true;
  for i := 1 to n do
  begin
    a := Losuj(2, p-2);
    x := PotegujModulo(a, d, p);
    if(x = 1) or (x = p-1) then continue;
    j := 1;
    while(j < s) and (x <> p-1) do
    begin
      x := MnozModulo(x, x, p);
      if x = 1 then
      begin
        t := false; break;
      end;
      inc(j);
    end;
    if not t then break;
    if x <> p-1 then
    begin
      t := false; break;
    end;
  end;
  if t then writeln('TAK')
  else writeln('NIE');
end.

C++

// Test Millera-Rabina
// Data: 6.04.2008
// (C)2020 mgr Jerzy Wałaszek
//---------------------------

#include <iostream>
#include <cstdlib>
#include <ctime>

using namespace std;

typedef unsigned long long ulong;

// 64 bitowy generator pseudolosowy
//---------------------------------
ulong Losuj(ulong a, 
            ulong b)
{
  ulong w;
  int i;

  for(i = 1; i <= 8; i++)
  {
    w <<= 8;
    w |= rand() % 256;
  }
  return a+(w % (b-a));
}

// Funkcja mnoży a i b mod n
//--------------------------
ulong MnozModulo(ulong a, 
                 ulong b, 
                 ulong n)
{
  ulong m, w;

  w = 0;
  for(m = 1; m; m <<= 1)
  {
    if(b & m) w = (w+a) % n;
    a = (a << 1) % n;
  }
  return w;
}

// Funkcja oblicza a^e mod n
//--------------------------
ulong PotegujModulo(ulong a, 
                    ulong e, 
                    ulong n)
{
  ulong m, p, w;

  p = a; w = 1;
  for(m = 1; m; m <<= 1)
  {
    if(e & m) w = MnozModulo(w, p, n);
    p = MnozModulo(p, p, n);
  }
  return w;
}

int main()
{
  ulong a, d, p, x;
  int i, j, s, n;
  bool t;

  srand(time(NULL));
  cin >> p >> n;
  s = 0;
  for(d = p-1; d % 2 == 0; s++) d /= 2;
  t = true;
  for(i = 1; i <= n; i++)
  {
    a = Losuj(2, p-2);
    x = PotegujModulo(a, d, p);
    if((x == 1) || (x == p-1)) continue;
    for(j = 1; (j < s) && (x != p-1); j++)
    {
      x = MnozModulo(x, x, p);
      if(x == 1)
      {
        t = false; break;
      }
    }
    if(!t) break;
    if(x != p-1)
    {
      t = false; break;
    }
  }
  cout << (t ? "TAK": "NIE") << endl;
  return 0;
}

Basic

' Test Millera-Rabina
' Data: 6.04.2008
' (C)2020 mgr Jerzy Wałaszek
'---------------------------

' 64 bitowy generator pseudolosowy
'---------------------------------
Function Losuj(Byval a As Ulongint, _
               Byval b As Ulongint) _
               As Ulongint

  Dim w As Ulongint, i As Integer

  For i = 1 To 8
    w = w Shl 8
    w = w Or Culngint(Rnd(1)*256)
  Next
  Losuj = a+(w Mod (b-a))
End Function

' Funkcja mnoży a i b mod n
'--------------------------
Function MnozModulo(Byval a As Ulongint, _
                    Byval b As Ulongint, _
                    Byval n As Ulongint) _
                    As Ulongint

  Dim As Ulongint m, w

  w = 0: m = 1
  While m
    If b And m Then w = (w+a) Mod n
    a = (a Shl 1) Mod n
    m = m Shl 1
  Wend
  MnozModulo = w
End Function

' Funkcja oblicza a^e mod n
'--------------------------
Function PotegujModulo(Byval a As Ulongint, _
                       Byval e As Ulongint, _
                       Byval n As Ulongint) _
                       As Ulongint

  Dim As Ulongint m, p, w

  p = a: w = 1: m = 1
  While m
    If e And m Then _
      w = MnozModulo(w, p, n)
    p = MnozModulo(p, p, n)
    m = m Shl 1
  Wend
  PotegujModulo = w
End Function

Dim As Ulongint a, d, p, x
Dim As Integer i, j, s, n
Dim t As Byte

Randomize
Input p, n
d = p-1
s = 0
While d Mod 2 = 0
  s += 1: d \= 2
Wend
t = 1
For i = 1 To n
  a = Losuj(2, p-2)
  x = PotegujModulo(a, d, p)
  if(x = 1) Or (x = p-1) Then _
    Continue For
  j = 1
  while(j < s) And (x <> p-1)
    x = MnozModulo(x, x, p)
    If x = 1 Then
      t = 0: Exit While
    End If
    j += 1
  Wend
  If t = 0 Then Exit For
  If x <> p-1 Then
    t = 0: Exit For
  End If
Next
If t = 1 Then
  Print "TAK"
Else
  Print "NIE"
End If
End

Python (dodatek)

# Test Millera-Rabina
# Data: 14.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ---------------------------

import random

G = 2**64


# 64 bitowy generator pseudolosowy
# ---------------------------------
def losuj(a, b):
    w = 0
    for i in range(8):
        w <<= 8
        w |= int(random.random()*256)
    return a+(w % (b-a))


# Funkcja mnoży a i b % n
# -----------------------
def mnoz_modulo(a, b, n):
    w, m = 0, 1
    while m < G:
        if b & m: w = (w+a) % n
        a = (a << 1) % n
        m <<= 1
    return w

# Funkcja oblicza a**e % n
# ------------------------
def poteguj_modulo(a, e, n):
    p, w, m = a, 1, 1
    while m < G:
        if e & m:
            w = mnoz_modulo(w, p, n)
        p = mnoz_modulo(p, p, n)
        m <<= 1
    return w

random.seed(None, 2)

arr = input().split(" ")
p = int(arr[0])
n = int(arr[1])

d, s = p-1, 0
while not (d % 2):
    s += 1
    d //= 2
t = True
for i in range(n):
    a = losuj(2, p-2)
    x = poteguj_modulo(a, d, p)
    if (x == 1) or (x == p-1):
        continue
    j = 1
    while (j < s) and (x != p-1):
        x = mnoz_modulo(x, x, p)
        if x == 1:
            t = False
            break
        j += 1
    if not t: break
    if x != p-1:
        t = False
        break
if t:
    print("TAK")
else:
    print("NIE")

Wynik:

9223372036854775783 20
TAK

9223372036854775787 1
NIE

do podrozdziału do strony

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone pod warunkiem podania źródła oraz niepobierania za to pieniędzy.
Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl
Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.

Informacje dodatkowe.

Pierwszość liczby naturalnej – test Millera-Rabina