Sito Eratostenesa

27.11.2009 - Anna Piekarska

Już starożytni Grecy zastanawiali się nad zagadnieniem liczb pierwszych. Aktualnie problem zdaje się wciąż nie być dostatecznie zbadany. Jednak dziś wiemy nieco więcej niż kiedyś. Do rozwoju teorii przyczynił się zdecydowanie Eratostenes - wymyślił tzw. Sito Eratostenesa czyli dość prosty algorytm znajdowania liczb pierwszych. Prosty oczywiście dla małych liczb. Wciąż nie znamy algorytmu wielomianowego na faktoryzację. Moim zdaniem na jednymi z ważniejszych późniejszych odkryć są szczególne rodzaje liczb pierwszych - np. Fermata czy Mersenne'a oraz probabilistyczne testy pierwszości.

Czego się nauczymy?

Dowiemy się jak sprawdzić, czy liczba jest pierwsza. Znajdziemy dzielniki pierwsze liczby, czyli dokonamy faktoryzacji. Poznamy powody, dla których liczby pierwsze są tak ważne. Zapoznamy się ze szkicem rozwiązania zadania "Najdzielniejszy dzielnik" z tegorocznej Olimpiady Informatycznej.

Sposób naiwny

Chcemy dowiedzieć się, czy liczba n jest pierwsza. Czemu w takim razie nie sprawdzić po prostu dla dużej ilości liczb czy są dzielnikami n? Tylko które liczby sprawdzić? Czy 15 może być dzielnikiem 10? Dlaczego?
Zauważmy, że dzielnik liczby n nie może być większy od n. Sprawdźmy więc liczby z przedziału $[2,n-1]$ :

bool isPrime(int n)
{
    for(int d = 2; d < n; d++) // sprawdzamy każdy dzielnik
    {
	 if(n % d == 0) // sprawdzamy, czy n dzieli się przez d
	 {
		cout << d << "jest dzielnikiem " << n;
                cout << ", więc " << n << " nie jest pierwsze" << endl;
		return 0; // jeśli tak, to n nie jest pierwsza
	 }
    }
    return 1; // n jest pierwsza
}

Przyspieszamy

Oszacujmy złożoność tego algorytmu. Pętlę wykonujemy n razy. W pętli wykonujemy stałą liczbę operacji, więc złożoność wynosi $O(n)$ . Trochę wolno. Czy możemy to przyspieszyć?
Po usunięciu komendy return 0; z powyższego kodu nasz program wypisze wszystkie dzielniki. Popatrzmy na wyniki dla $n = 24$ .

2 jest dzielnikiem 24, więc 24 nie jest pierwsze
3 jest dzielnikiem 24, więc 24 nie jest pierwsze
4 jest dzielnikiem 24, więc 24 nie jest pierwsze
6 jest dzielnikiem 24, więc 24 nie jest pierwsze
8 jest dzielnikiem 24, więc 24 nie jest pierwsze
12 jest dzielnikiem 24, więc 24 nie jest pierwsze

Czy na pewno musimy sprawdzać liczbę 12? Przecież skoro 2 dzieli 24 to $\frac{24}{2}=12$ też.

Uogólniając zauważamy, że jeśli d jest dzielnikiem n to $\frac{n}{d}$ też. W takim razie każdemu dzielnikowi z przedziału $[\sqrt{n},n]$ możemy przyporządkować dzielnik z przedziału $[2,\sqrt{n}]$ taki, że ich iloczyn wyniesie n. Poprawmy więc nasz algorytm:

bool isPrime(int n)
{
    for(int d = 2; d*d <= n; d++) 
    // sprawdzamy każdy dzielnik do pierwiastka
    {
	if(n % d == 0)
        // sprawdzamy, czy n dzieli się przez d
	{
		cout << d << "jest dzielnikiem " << n;
                cout << ", więc " << n << " nie jest pierwsze" << endl;
		return 0; // jeśli tak, to n nie jest pierwsza
	}
     }
     return 1; // n jest pierwsza
}

Jaka złożoność jest teraz? Tym razem pętlę wykonujemy $\sqrt{n}$ razy, w każdej stałą liczbę operacji. Przyspieszyliśmy zatem do $O(\sqrt{n})$ . Na pojedyńcze sprawdzenie może wystarczy.

Więcej zapytań

A co jeśli dostaniemy więcej zapytań? Powiedzmy milion. Załóżmy, że nasze zapytania też są w okolicach miliona. Wykonując za każdym razem nasz algorytm, wykonamy około $10^6*\sqrt{10^6} = 10^9$ . To za wolne.
Spróbujmy na raz dla wszystkich liczb z przedziału $[2,n]$ sprawdzić, czy są pierwsze. Nazywa się to Sitem Eratostenesa. Dla każdej liczby pierwszej wykreślmy jej wielokrotności z listy - nie bądą one przecież pierwsze.

Przykładowo załóżmy $n=30$ .

Liczba 2 jest pierwsza, wykreślmy więc wielokrotności liczby 2.
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Liczba 3 jest pierwsza, wykreślmy więc wielokrotności liczby 3.
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Liczba 5 jest pierwsza, wykreślmy więc wielokrotności liczby 5.
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Liczba 7 jest pierwsza, wykreślmy więc wielokrotności liczby 7. Chwila! Nie mamy już nieskreślonych wielokrotności liczby 7. Dlaczego? Popatrzmy, jakie liczby skreślaliśmy w kolejnych krokach. Zauważmy, że w każdym kroku liczby, które skreślamy nie mają dzielników pierwszych mniejszych niż aktualnie przetwarzana liczba pierwsza. Wydaje się to być logiczne - wszystkie liczby z mniejszymi czynnikami pierwszymi zostały już skreślone! Czyli ostatnia przetwarzana przez nas liczba pierwsza nie powinna przekroczyć $\sqrt{n}$ bo większe i tak nam nic nowego nie wprowadzą.

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

const int MAX = 1000010; // tyle liczb chcemy sprawdzić
bool p[MAX]; // tablica, w której zaznaczamy, czy liczba jest pierwsza
void init()
{
   for(int i = 2; i < MAX; i++) p[i] = 1; 
   // na początku wszystkie liczby są pierwsze
   for(int i = 2; i*i < MAX; i++)
   {
	if(p[i])
	{
		for(int j = 2; i * j < MAX; j++) p[i*j] = 0; 
                // i*j jest wielokrotnością i, więc jest liczbą złożoną
	}
   }
}

bool isPrime(int n)
{
   if(n < MAX) return p[n];
   else cout << "Nie wiem. Liczba jest za duza" << endl;
}

Dlaczego to działa?

Oczywistym jest, że żadna liczba pierwsza nie została oznaczona jako złożona. Zastanówmy się, czy może istnieć taka liczba złożona, która nie została znaleziona. Jeśli udowodnimy, że nie ma takiej, będzie to oznaczało, że algorytm jest poprawny. Załóżmy nie wprost, że pewna liczba postaci $a*b \leq MAX$ pozostanie nieoznaczona przez nasz algorytm. Skoro $a*b \leq MAX$ oraz $a \geq 2$ to $b \leq MAX$ .

Jednak wszystkie liczby postaci $2b, 3b, 4b, ...$ oznaczyliśmy jako złożone, co prowadzi do sprzeczności a więc wszystkie liczby złożone faktycznie oznaczyliśmy jako złożone.

Czy to działa szybciej?

Musimy policzyć, jak szybko działa nasz algorytm. W końcu po co nam skomplikowany algorytm wolniejszy od tego prostego, który już znamy.

Nasz algorytm dla każdej liczby $2 \leq d \leq \sqrt{MAX}$ wykona $\frac{MAX}{d}$ operacji. W sumie będzie to $MAX * (\frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... \frac{1}{\sqrt{MAX}})$ . Niech $W_{n} = \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... \frac{1}{n}$ . Spróbujmy oszacować $W_{n}$ . Pogrupujmy wyrazy tak, żeby pierwszy (największy) wyraz w grupie wynosił $\frac{1}{2^k}$ dla kolejnych $k \in \mathbb{N}$ . Grup będzie $O(log n)$ .
$W_{n} = \frac{1}{2} + \frac{1}{3} + \frac{1}{4} + ... \frac{1}{n} = (\frac{1}{2} + \frac{1}{3}) + (\frac{1}{4} + \frac{1}{5} + \frac{1}{6} + \frac{1}{7}) + (\frac{1}{8} + \frac{1}{9} + \frac{1}{10} + ...) + ... \leq (\frac{1}{2} + \frac{1}{2}) + (\frac{1}{4} + \frac{1}{4} + \frac{1}{4} + \frac{1}{4}) + (\frac{1}{8} + \frac{1}{8} + \frac{1}{8} + ...) + ... = 2 * \frac{1}{2} + 4 * \frac{1}{4} + 8 * \frac{1}{8} + ... = O(log_{2}n) * 1 = O(log_{2}n)$