Algorytmy i Struktury Danych - Małe Twierdzenie Fermata

W roku 1640 Pierre de Ferm at odkrył ponownie chińskie twierdzenie o pierwszości:

Jeśli p jest liczbą pierwszą, to 2^p mod p = 2.

Równanie jest spełniane przez wszystkie bez wyjątku liczby pierwsze. Jeśli liczba p nie spełnia równania 2^pmod p = 2, to nie jest liczbą pierwszą – to jest 100% pewne. Pozwala to eliminować w prosty sposób liczby złożone. Jednakże spełnianie tego równania przez p nie gwarantuje, iż p jest pierwsze, ponieważ istnieją liczby złożone, które spełniają równanie. Nazywamy je liczbami pseudopierwszymi przy podstawie 2 (ang. base-2 pseudoprimes).

Przykład:

p = 341 = 11 × 31
2³⁴¹ mod 341 = 2

Pomyślisz w tym miejscu, iż w takim razie twierdzenie jest niepoprawne. Nic z tych rzeczy. W matematyce relacja wynikania nie jest zwrotna:

Z A wynika B – nie jest tym samym, co z B wynika A – wtedy mielibyśmy relację równoważności.

Twierdzenie mówi jedynie, iż jeśli liczba p jest pierwsza, to 2^p jest przystające modulo p do 2. Czyli jest spełnione dla każdej liczby pierwszej (z A wynika B). Jeśli jednak weźmiemy drugą część twierdzenia, czyli równanie 2^p mod p = 2, i znajdziemy liczbę p, która je spełnia, to wcale nie mamy gwarancji, iż p będzie liczbą pierwszą (z B nie wynika A):

Wszyscy Polacy mówią po polsku (z A wynika B), ale jeśli ktoś mówi po polsku, to wcale nie musi być Polakiem (z B nie wynika A).

Fermat dokładnie przestudiował chińskie twierdzenie o pierwszości oraz zbadał je dla innych podstaw niż 2. W efekcie znacznie poprawił dokładność chińskiego testu. Wyniki pracy Fermata są dzisiaj znane pod nazwą Małego Twierdzenia Fermata (ang. Fermat's Little Theorem – FLT). Można je przedstawić następująco:

Małe Twierdzenie Fermata

Niech p będzie dowolną liczbą pierwszą i a dowolną liczbą naturalną.

Jeśli NWD(p, a) = 1, to a^p-1 mod p = 1

Małe Twierdzenie Fermata jest bardziej ogólne od twierdzenia chińskiego. Jeśli w FLT podstawę a zastąpimy przez 2 i obie strony równania pomnożymy przez 2, to otrzymamy chińskie twierdzenie o pierwszości:

a = 2
a^p-1 mod p = 1 / ×2
(2^p-1)×2 mod p = 1×2
2^p-1+1 mod p = 2
2^p mod p = 2

Wzór Fermata, a^p-1 mod p = 1, pozwala przetestować liczbę p dla kilku różnych podstaw, co w olbrzymim stopniu poprawia rzetelność testu na pierwszość. Oczywiście Małe Twierdzenie Fermata również jest nawiedzane przez liczby pseudopierwsze przy podstawie a. Jednakże dla różnych a są to zwykle różne liczby. Zatem testy dla kilku podstaw a pozwalają je zwykle znakomicie wyeliminować. Jednakże należy tutaj wspomnieć o istnieniu pewnych liczb złożonych, zwanych liczbami Carmichaela, które są pseudopierwsze dla każdej podstawy a – test Fermata ich nie wyeliminuje. Na szczęście liczby Carmichaela są bardzo rzadkie i odległe od siebie – mniej więcej jedna co miliard liczb pierwszych – przykładowo poniżej 10¹⁶ jest tylko 246'683 liczb Carmichaela, natomiast liczb pierwszych jest 279'238'341'033'925.

Algorytm Fermata testujący pierwszość liczby p wybiera liczbę losową a w zakresie od 2 do p-1. Następnie sprawdza, czy jest ona względnie pierwsza z p, czyli czy NWD(p, a) = 1. Jeśli tak, to testowany jest warunek: a^p-1 mod p = 1. Jeśli liczba p nie spełnia tego warunku, to na pewno nie jest liczbą pierwszą. Jeśli spełnia go, to jest albo liczbą pierwszą, albo pseudopierwszą. Aby wykluczyć większość liczb pseudopierwszych można test Fermata wykonać kilkakrotnie, np. 10 razy. Jeśli liczba p nie zostanie odrzucona, to z bardzo dużym prawdopodobieństwem jest liczbą pierwszą lub liczbą pseudopierwszą Carmichaela. Te ostatnie można wyeliminować w ponad 95% sprawdzając podzielność liczby p przez liczby pierwsze z przedziału od 2 do 1000. Jeśli liczba p przejdzie pomyślnie te testy, to prawie na pewno jest liczbą pierwszą.

Do wyznaczania a^p-1 zastosujemy algorytmy mnożenia i potęgowania modulo n z poprzedniego rozdziału. Sposób obliczania NWD podaliśmy w rozdziale opisującym algorytm Euklidesa. Generację liczb pseudolosowych opisujemy w następnych rozdziałach.

Algorytm testu pierwszości Fermata

Wejście:

p : testowana liczba, p ∈ N, p > 2.

Wyjście:

TAK : liczba p jest pierwsza lub z bardzo małym prawdopodobieństwem pseudopierwsza Carmichaela.
NIE : liczba p nie jest pierwsza.

Elementy pomocnicze:

a : podstawa potęgi, a ∈ N.
i : zlicza ilość wykonanych testów Fermata, i ∈ N.
NWD(a, b) : oblicza NWD liczb a i b.
Losowa(a, b) : zwraca liczbę pseudolosową z zakresu od a do b-1.

Lista kroków:

K01: Sprawdź, czy p jest podzielne przez    ; wstępna eliminacja liczb złożonych oraz
     liczbę pierwszą z przedziału <2, 1000>. ; ewentualnych liczb pseudopierwszych Carmichaela
     Jeśli tak, to idź do kroku K07         ; Dla liczb złożonych zwracamy NIE i kończymy
K02: Dla i = 1, 2, …, 10: ; w pętli sprawdzamy warunek Fermata
     wykonuj kroki K03…K05
K03:   a ← Losowa(2, p-1) ; wybieramy losową podstawę a
K04:   Jeśli NWD(p, a) ≠ 1,  ; podstawa musi być względnie pierwsza z p
       to idź do kroku K07
K05:   Jeśli a^p-1 mod p ≠1,  ; test na liczbę złożoną
       to idź do kroku K07
K06: Pisz "TAK" i zakończ  ; liczba p przeszła pozytywnie wszystkie testy
K07: Pisz "NIE" i zakończ  ; liczba p nie przeszła testów, jest zatem liczbą złożoną

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, a w drugim wierszu wypisuje:

TAK – jeśli liczba p jest pierwsza lub z bardzo małym prawdopodobieństwem pseudopierwsza Carmichaela
NIE – jeśli liczba p jest liczbą złożoną

Pascal

// Test pierwszości Fermata
// Data: 5.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
  w := 0;
  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 oblicza NWD
//--------------------
function NWD(a, b : qword) : qword;

var t : qword;

begin
  while b <> 0 do
  begin
    t := b;
    b := a mod b;
    a := t;
  end;
  NWD := 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;

// Tablica początkowych 169 liczb pierwszych
//------------------------------------------
const lp : array [0..168] of integer = (
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009);

var p, a : qword;
    i : integer;
    t : boolean;
begin
  randomize;
  readln(p);
  t := true;
  for i := 0 to 168 do
    if(p <> lp[i]) and (p mod lp[i] = 0) then
    begin
      t := false; break;
    end;
  if t and (p > 1009) then
  begin
    for i := 1 to 10 do
    begin
      a := Losuj(2, p-1);
      if(NWD(p, a) <> 1) or (PotegujModulo(a, p-1, p) <> 1) then
      begin
        t := false; break;
      end;
    end;
  end;
  if t then writeln('TAK')
  else      writeln('NIE');
end.

C++

// Test pierwszości Fermata
// Data: 5.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 = 0;
  int i;

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

// Funkcja oblicza NWD
//--------------------
ulong NWD(ulong a, ulong b)
{
  ulong t;

  while(b)
  {
    t = b; b = a % b; a = t;
  }
  return 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;
}

// Tablica początkowych 169 liczb pierwszych
//------------------------------------------
const ulong lp[] = {
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, 
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009};

int main()
{
  ulong p, a;
  int   i;
  bool  t;

  srand((unsigned)time(NULL));
  cin >> p;
  t = true;
  for(i = 0; i < 169; i++)
    if((p != lp[i]) && (p % lp[i] == 0))
    {
      t = false; break;
    }
  if(t && (p > 1009))
  {
    for(i = 1; i <= 10; i++)
    {
      a = Losuj(2, p-1);
      if((NWD(p, a) != 1) || (PotegujModulo(a, p-1, p) != 1))
      {
        t = false; break;
      }
    }
  }
  cout << (t ? "TAK": "NIE") << endl;
  return 0;
}

Basic

' Test pierwszości Fermata
' Data: 5.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

  w = 0
  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 oblicza NWD
'--------------------
Function NWD(Byval a As Ulongint, _
             Byval b As Ulongint) _
             As Ulongint
  Dim t As Ulongint

  While b
    t = b: b = a Mod b: a = t
  Wend
  NWD = 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

' Tablica początkowych 169 liczb pierwszych
'------------------------------------------
Dim lp(168) As Integer => {_
  2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, _
 71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, _
167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, _
271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, _
389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, _
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, _
631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, _
757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, _
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009}

Dim As Ulongint p, a
Dim i As Integer, t As Byte

Randomize

Open Cons For Input As #1
Input #1, p
Close #1

t = 1

For i = 0 To 168
  if(p <> lp(i)) And (p Mod lp(i) = 0) Then
    t = 0: Exit For
  End If
Next

if (t = 1) And (p > 1009) Then
  For i = 1 To 10
    a = Losuj(2, p-1)
    if(NWD(p, a) <> 1) OrElse (PotegujModulo(a, p-1, p) <> 1) Then
      t = 0: Exit For
    End If
  Next
End If

If t = 1 Then Print "TAK": Else Print "NIE"

End

Python (dodatek)

# Test pierwszości Fermata
# 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 oblicza NWD
# -------------------
def nwd(a, b):
    while b:
        t = b
        b = a % b
        a = t
    return a


# Funkcja mnoży a×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


# Tablica początkowych 169 liczb pierwszych
# ------------------------------------------
lp = [ \
      2,   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,  53,  59,  61,  67, \
     71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, \
    167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, \
    271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, \
    389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, \
    503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, \
    631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, \
    757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, \
    883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009]

random.seed(None, 2)

p = int(input())
t = True
for i in range(169):
    if (p != lp[i]) and not (p%lp[i]):
        t = False
        break
if t and (p > 1009):
    for i in range(10):
        a = losuj(2, p-1)
        if (nwd(p, a) != 1) or (poteguj_modulo(a, p-1, p) != 1):
            t = False
            break
if t:
    print("TAK")
else:
    print("NIE")

Wynik:

9223372036854775783
TAK

9223372036854775787
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 – Małe Twierdzenie Fermata