Algorytmy i Struktury Danych - Chiński test pierwszości

Problem

Sprawdzić, czy liczba naturalna p jest pierwsza.

Zainteresowanie testami pierwszości (ang. primality testing) wzrosło w ostatnich latach z uwagi na powszechne wprowadzenie do użytku różnych publicznych systemów kryptograficznych. Systemy te wymagają szybkich algorytmów weryfikacji, czy dana liczba naturalna jest liczbą pierwszą – liczby pierwsze są potrzebne do tworzenia tzw. kluczy szyfrujących. Z uwagi na charakter procesu szyfrowania informacji testowane liczby są bardzo duże. Zatem zwykłe metody testowania pierwszości, np. przez próbne dzielenia liczby przez czynniki, nie zdają egzaminu z powodu swojej powolności. Dla przykładu rozważmy test pierwszości liczby 200 cyfrowej. Możemy ją przybliżyć następująco:

p ≈ 10²⁰⁰

Aby sprawdzić, czy liczba p jest pierwsza, dzielimy ją próbnie przez liczby z przedziału <2, √p>. Liczb tych jest:

d ≈ 10¹⁰⁰

Jeśli wyeliminujemy z tego przedziału liczby podzielne przez 2 lub przez 3, to pozostanie ich ¹/₃ w stosunku do pierwotnej liczby, zatem:

d' ≈ 10¹⁰⁰ / 3

Tyle musimy wykonać dzieleń próbnych w algorytmie naiwnym. Wyobraźmy sobie, iż dysponujemy superkomputerem, który w ciągu 1 sekundy jest w stanie wykonać sto miliardów (10¹¹) takich testów podzielności (zwykłe Pentium IV tego nie potrafi, być może jakiś super Cray przyszłości lub komputer kwantowy). Zatem wykonanie operacji zajmie:

t ≈ 10¹⁰⁰ / 10¹¹ =
    10¹⁰⁰ / (3 × 10¹¹) =
    10⁸⁹ / 3 [s] ≈
    3,33 × 10⁸⁸ [s] ≈
    10⁸¹ [lat]

Widzimy wyraźnie, iż czas obliczeń jest niewyobrażalnie duży – cały nasz Wszechświat istnieje "zaledwie" 1, 3 × 10¹⁰ lat. Musimy zatem sięgnąć do innych metod testowania pierwszości.

Ideą takich metod jest wyszukiwanie własności, które spełniają tylko liczby pierwsze lub własności spełnianych jedynie przez liczby złożone. Następnie sprawdzamy, czy badana liczba spełnia określoną własność i na tej podstawie wnioskujemy o jej pierwszości. Pierwsza wskazówka pojawiła się około 500 lat p.n.e. w Chinach, gdzie odkryto fakt, iż jeśli liczba p jest pierwsza, to wyrażenie 2^p – 2 jest podzielne przez p. Sprawdźmy dla kilku początkowych liczb pierwszych:

p =  2, 2^p - 2 = 2²  - 2 =    4 - 2 =    2 =   1 × p
p =  3, 2^p - 2 = 2³  - 2 =    8 - 2 =    6 =   2 × p
p =  5, 2^p - 2 = 2⁵  - 2 =   32 - 2 =   30 =   6 × p
p =  7, 2^p - 2 = 2⁷  - 2 =  128 - 2 =  126 =  18 × p
p = 11, 2^p - 2 = 2¹¹ - 2 = 2048 - 2 = 2046 = 186 × p …

Powyższą własność można zapisać w postaci twierdzenia:

Niech p będzie liczbą naturalną większą od 1. Jeśli liczba 2^p nie jest przystająca do liczby 2 modulo p, to p jest liczbą złożoną.

Innymi słowy, jeśli 2^p mod p ≠ 2, to p jest złożone.

Jeśli twierdzenie jest spełnione, to mamy pewność, iż p jest liczbą złożoną. Niestety, jeśli jednak twierdzenie nie jest spełnione, czyli 2^p mod p = 2, to liczba p jest albo pierwszą (wszystkie liczby pierwsze p dzielą 2^p - 2), albo pseudopierwszą przy podstawie 2 (ang. base 2 pseudoprime), czyli złożoną dzielącą 2^p - 2. Na szczęście liczby pseudopierwsze przy podstawie 2 występują bardzo rzadko. Jeśli zatem liczba p nie spełnia testu, to istnieje tylko prawdopodobieństwo około 0,002%, iż mimo wszystko p jest liczbą złożoną.

Podstawowym problemem chińskiego testu pierwszości jest obliczanie dużych potęg liczby 2. Na pierwszy rzut oka możemy odnieść wrażenie, iż zadanie nie jest wykonalne w zakresie liczb oferowanym przez podstawowe typy danych. Na przykład liczby 64-bitowe posiadają jedynie zakres od 0 do 2⁶⁴ - 1. Jednakże nas nie interesuje konkretna wartość potęgi liczby 2, a jedynie reszta z dzielenia tej potęgi przez p. Możemy zatem wykorzystać zasady arytmetyki modularnej (ang. modular arithmetic), w której wynik każdej operacji jest resztą z dzielenia przez moduł.

Wynik operacji arytmetycznej modulo n wpada w przedział od 0 do n - 1. Zatem będzie tego samego rzędu, co n. Dzięki tej własności możemy obliczać reszty potęg dowolnej wielkości – obliczanie potęgi rozbijamy na operacje mnożenia modulo n. W efekcie potęgowanie sprowadza się do wymnażania odpowiednich reszt. Idea opiera się na prostej zależności:

a^x+^y mod n = (a^x mod n) × (a^y mod n) mod n

Oznacza to, iż zamiast wyliczać resztę potęgi a^x+^y mod n, możemy wyznaczyć reszty modulo n potęg a^x oraz a^y i wymnożyć je modulo n. Dla potęg a^x i a^y można dokonać rekurencyjnie takiego samego rozkładu.

Przykład:

Obliczyć 2³⁸⁷ mod 13.

Rozbijamy liczbę 387 na sumę potęg 2:

387 = 256 + 128 + 2 + 1 2³⁸⁷ mod 13 =
(2²⁵⁶ mod 13) × (2¹²⁸ mod 13) × (2² mod 13) × (2¹ mod 13) mod 13

Teraz obliczamy reszty potęg 2:

2¹   mod 13 = 2 mod 13 = 2
2²   mod 13 = (2¹   mod 13) × (2¹   mod 13) mod 13 = (2 × 2) mod 13 =  4 mod 13 = 4
2⁴   mod 13 = (2²   mod 13) × (2²   mod 13) mod 13 = (4 × 4) mod 13 = 16 mod 13 = 3
2⁸   mod 13 = (2⁴   mod 13) × (2⁴   mod 13) mod 13 = (3 × 3) mod 13 =  9 mod 13 = 9
2¹⁶  mod 13 = (2⁸   mod 13) × (2⁸   mod 13) mod 13 = (9 × 9) mod 13 = 81 mod 13 = 3
2³²  mod 13 = (2¹⁶  mod 13) × (2¹⁶  mod 13) mod 13 = (3 × 3) mod 13 =  9 mod 13 = 9
2⁶⁴  mod 13 = (2³²  mod 13) × (2³²  mod 13) mod 13 = (9 × 9) mod 13 = 81 mod 13 = 3
2¹²⁸ mod 13 = (2⁶⁴  mod 13) × (2⁶⁴  mod 13) mod 13 = (3 × 3) mod 13 =  9 mod 13 = 9
2²⁵⁶ mod 13 = (2¹²⁸ mod 13) × (2¹²⁸ mod 13) mod 13 = (9 × 9) mod 13 = 81 mod 13 = 3

Mając reszty potęg (zaznaczone na czerwono), wymnażamy je modulo 13:

2³⁸⁷ mod 13 = (2²⁵⁶ mod 13) × (2¹²⁸ mod 13) × (2² mod 13) × (2¹ mod 13) mod 13
2³⁸⁷ mod 13 = (3 × 9 × 4 × 2) mod 13 = 216 mod 13 = 8

Zwróć uwagę, iż w obliczeniu końcowym biorą udział potęgi liczby 2, które odpowiadają wagom bitów ustawionych na 1 w reprezentacji dwójkowej wykładnika potęgowego:

387₍₁₀₎ = 110000011₍₂₎

Spostrzeżenie to wykorzystuje poniższy algorytm szybkiego wyznaczania reszty modulo n potęgi liczby 2. Bada on stan bitów wykładnika i na tej podstawie tworzy wynik. Mnożenie modulo n można rozbić na dodawanie modulo n wielokrotności mnożnej. Dzięki temu wynik częściowy nie wykroczy wiele poza moduł i operację będzie można przeprowadzać dla dużych liczb.

Algorytm obliczania a × b mod n

Wejście:

a, b : mnożone liczby; a, b ∈ N.
n : moduł operacji mnożenia; n ∈ N.

Wyjście:

w = a × b mod n.

Zmienne pomocnicze:

m : maska bitowa służąca do testowania ustawionych w b bitów.

Lista kroków:

K01: w ← 0  ; wynik mnożenia
K02: m ← 1  ; maska bitowa
K03: Dopóki m ≠ 0:
     wykonuj kroki K04…K06
K04:   Jeśli b and m ≠ 0, ; sprawdzamy ustawione bity w mnożniku
       to w ← w + a mod n
K05:   a ← ( a shl 1) mod n ; przesuwamy mnożną o 1 bit w lewo
K06:   m ← m shl 1 ; przesuwamy bit w masce
K07: Pisz w
K08: Zakończ

Algorytm obliczania a^e mod n

Wejście:

a : podstawa potęgi, a ∈ N.
e : wykładnik dla potęgi a. e ∈ N.
n : moduł. n ∈ N. n > 1.

Wyjście:

w = 2^e mod n.

Zmienne pomocnicze:

m : maska bitowa służąca do testowania ustawionych w e bitów.
p : reszta z kolejnych potęg a modulo n. p ∈ N.

Lista kroków:

K01: p ← a                ; wartość początkowa potęgi a¹K02: w ← 1                ; wynik
K03: m ← 1                ; maska bitowa – ustawiony bit b₀K04: Dopóki m ≠ 0:        ; obliczamy 2^e mod n
     wykonuj kroki K05…K07
K05:   Jeśli e and m ≠ 0, ; sprawdzamy, czy e posiada ustawiony bit na pozycji z m
       to w ← w × p mod n ; jeśli tak, to wynik wymnażamy przez potęgę
K06:   p ← p × p mod n    ; wyliczamy następną potęgę
K07:   m ← m shl 1        ; przesuwamy w lewo bit maski
K08: Pisz w
K09: Zakończ

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 dwie liczby: e – wykładnik potęgi liczby 2 oraz n – moduł. W drugim wierszu wypisuje wynik 2^e mod n.

Pascal

// Reszty potęg 2 modulo n
// Data   : 4.04.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

program prg;

// 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;

var e, n, m, p, w : qword;

begin
  readln(e, n);
  p := 2; 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;
  writeln(w);
end.

C++

// Reszty potęg 2 modulo n
// Data   : 4.04.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>

using namespace std;

typedef unsigned long long ulong;

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

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

int main()
{
  ulong e, n, m, p, w;

  cin >> e >> n;
  p = 2; w = m = 1;
  while(m)
  {
    if(e & m) w = MnozModulo (w, p, n);
    p = MnozModulo (p, p, n);
    m <<= 1;
  }
  cout << w << endl;
  return 0;
}

Basic

' Reszty potęg 2 modulo n
' Data   : 4.04.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

' 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) <> 0 Then w = (w + a) Mod n
    a = (a Shl 1) Mod n
    m = m Shl 1
  Wend
  MnozModulo = w
End Function

Dim As Ulongint e, n, m, p, w

Input e, n
p = 2: w = 1: m = 1
While m
  if(e And m) <> 0 Then w = MnozModulo (w, p, n)
  p = MnozModulo (p, p, n)
  m = m Shl 1
Wend
Print w
End

Python (dodatek)

# Reszty potęg 2 modulo n
# Data   : 14.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ----------------------------

# sztuczna granica dla int
# ----------------------------
G = 2 ** 64


# Funkcja mnoży a i b mod 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


# ----------------------------

arr = input().split(" ")
e = int(arr[0])
n = int(arr[1])
p, w, m = 2, 1, 1
while m < G:
    if e & m: w = mnoz_modulo(w, p, n)
    p = mnoz_modulo(p, p, n)
    m <<= 1
print(w)

W Pythonie nie ma ograniczenia na wielkość liczb całkowitych, powyższy program wprowadza za pomocą stałej G (granica) sztuczne ograniczenie do liczb 64 bitowych.

Wynik:

283125637241 551842729812
502162450484

Na początek: podrozdziału strony

Po podaniu algorytmów mnożenia modulo oraz potęgowania modulo przy podstawie 2 chiński test pierwszości (ang. Chinese Primality Test) staje się bardzo prosty. Odczytujemy liczbę p. Obliczamy 2^p mod p. Jeśli reszta jest różna od 2, to liczba p na pewno jest liczbą złożoną. Jeśli reszta jest równa 2, to z dużym prawdopodobieństwem p jest liczbą pierwszą.

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 w pierwszym wierszu czyta liczbę p i w następnym wierszu wypisuje słowo NIE, jeśli p jest liczbą złożoną lub TAK, jeśli liczba p jest pierwsza lub pseudopierwsza przy podstawie 2. Zwróć uwagę, iż testowanie jest błyskawiczne.

Pascal

// Chiński test pierwszości
// Data   : 4.04.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

program prg;

// 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 2^e mod n
//--------------------------
function PotegujModulo(e, n : qword) : qword;

var m, p, w : qword;

begin
  p := 2; 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 p : qword;

begin
  readln(p);
  if PotegujModulo (p, p) = 2 then writeln('TAK') else writeln('NIE');
end.

C++

// Chiński test pierwszości
// Data   : 4.04.2008
// (C)2020 mgr Jerzy Wałaszek
//----------------------------

#include <iostream>

using namespace std;

typedef unsigned long long ulong;

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

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

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

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

int main()
{
  ulong p;

  cin >> p;
  cout << (PotegujModulo (p, p) == 2 ? "TAK": "NIE") << endl;
  return 0;
}

Basic

' Chiński test pierwszości
' Data   : 4.04.2008
' (C)2020 mgr Jerzy Wałaszek
'----------------------------

' 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 2^e mod n
'--------------------------
Function PotegujModulo (Byval e As Ulongint,
                        Byval n As Ulongint)
                        As Ulongint
  Dim As Ulongint m, p, w

  p = 2: 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 p As Ulongint

Input p
If PotegujModulo (p, p) = 2 Then Print "TAK": Else Print "NIE"
End

Python (dodatek)

# Chiński test pierwszości
# Data   : 14.02.2024
# (C)2024 mgr Jerzy Wałaszek
# ----------------------------

G = 2 ** 64  # sztuczna granica int


# Funkcja mnoży a i b mod 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 2**e mod n
# ---------------------------
def poteguj_modulo(e, n):
    p, w, m = 2, 1, 1
    while m < G:
        if e & m: w = mnoz_modulo(w, p, n)
        p = mnoz_modulo(p, p, n)
        m <<= 1
    return w


# --------------------------

p = int(input())
if poteguj_modulo(p, p) == 2:
    print("TAK")
else:
    print("NIE")

Wynik:

9223372036854775783
TAK

9223372036854775787
NIE

Na początek: podrozdziału 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 – Chiński Test Pierwszości

Problem

Wyznaczanie dużych potęg liczby 2 modulo n

Algorytm obliczania a × b mod n

Wejście:

Wyjście:

Zmienne pomocnicze:

Lista kroków:

Algorytm obliczania a^e mod n

Wejście:

Wyjście:

Zmienne pomocnicze:

Lista kroków:

Przykładowe programy

Chiński test pierwszości

Przykładowe programy

Pierwszość liczby naturalnej – Chiński Test Pierwszości

Problem

Wyznaczanie dużych potęg liczby 2 modulo n

Algorytm obliczania a × b mod n

Wejście:

Wyjście:

Zmienne pomocnicze:

Lista kroków:

Algorytm obliczania ae mod n

Wejście:

Wyjście:

Zmienne pomocnicze:

Lista kroków:

Przykładowe programy

Chiński test pierwszości

Przykładowe programy

Algorytm obliczania a^e mod n