Niniejszy wpis łączy w sobie kilka ulubionych przeze mnie tematów, a więc spodziewam się z niego całkiem sporej ilości frajdy.
Przynajmniej dla mnie...
Tytuł dzisiejszego wpisu pochodzi ze strony http://projecteuler.net - jest to strona z zadaniami matematyczno-logicznymi idealnymi do ćwiczenia różnych języków programowania.
Dzisiaj zajmiemy się zagadką numer 10, czyli będziemy sumować liczby pierwsze. Ha!
Żeby radości było jeszcze więcej, sumowanie owo wykonamy za pomocą języka Python, który próbuję od jakiegoś czasu intensywnie przyswajać, pomimo tego, iż nie używam go (póki co!) w pracy. No bo po cholerę mi na zmywaku jakiś Python?
Zadanie Euler 10 brzmi następująco:
Wyznaczyć sumę wszystkich liczb pierwszych mniejszych od dwóch milionów.
Jak widać, bez pomocy komputera może być raczej ciężko. Liczby pierwsze mniejsze od stu, czy nawet od tysiąca, może jeszcze jakoś dałoby się zsumować "na piechotę", zużywając dużo ołówków i papieru. Ale dwa miliony... nie, tu trzeba końputra.
Zaczniemy od algorytmu naiwnego, który wygląda, z grubsza, o tak:
1. Ustaw sumę S na 2
2. Ustaw X = 3
3. Sprawdź, czy X jest pierwsze, jeżeli tak, ustaw S = S + X
4. Zwiększ X o 2
5. Jeżeli X jest mniejsze od 2 mln, przejdź do punktu 3
6. Wynikiem jest uzyskana na koniec suma S
W powyższym algorytmie najbardziej czasochłonny jest punkt 3, czyli sprawdzenie pierwszości liczby X. Użyjemy tutaj metody polegającej na próbie podzielenia liczby X przez każdą liczbę nieparzystą od 3 do pierwiastka z X powiększonego o 1. Jeżeli w trakcie dzielenia trafimy na resztę zero, wówczas liczba jest złożona i przechodzimy do kolejnej liczby X = X + 2.
Kod w Pythonie wygląda następująco:
# -*- coding: utf-8 -*-
import time
# Funkcja sprawdzająca pierwszość liczby
# używając algorytmu naiwnego: jeżeli parzysta i większa od dwójki, odrzucamy,
# jeżeli nieparzysta, sprawdzamy podzielność przez wszystkie
# nieparzyste od trójki aż do pierwiastka z liczby sprawdzanej + 1
def czy_pierwsza(liczba):
if liczba <= 2: return False if liczba == 2: return True if liczba % 2 == 0: return False for i in range(3, round(liczba ** 0.5) + 1, 2): if liczba % i == 0: return False return True # Generator zwracający kolejne liczby pierwsze. def generator_pierwszych(): yield 2 liczba = 3 while True: yield liczba liczba += 2 while not czy_pierwsza(liczba): liczba += 2 # Tu się zaczyna wykonywanie kodu. # Zapamiętujemy czas rozpoczęcia wykonywania. poczatek = time.time() # Zerujemy sumę suma = 0 # Dodajemy do siebie wszystkie liczby pierwsze mniejsze od dwóch milionów. for pierwsza in generator_pierwszych(): if pierwsza >= 2000000:
break
suma += pierwsza
# Wyliczamy czas, jaki upłynął od początku programu do teraz.
czas_wykonania = (time.time() - poczatek)
# Wyświetlamy wynik.
print("wynik: %s, czas wykonania: %s sekund" % (suma, czas_wykonania))
#### wynik: 142913828922, czas wykonania: 20.322700023651123 sekund
#### wynik: 142913828922, czas wykonania: 13.306761026382446 sekund
Powyższy kod zapuściłem na dwóch komputerach: na jednym wykonał się w odrobinę powyżej 20 sekund, na drugim w około 13.3 s. Wynik (na szczęście!) jest ten sam, czyli 142 913 828 922 (i jest to, przynajmniej według Google, wynik poprawny).
Spróbujemy teraz nieco innego podejścia i rozwiążemy ten sam problem przy użyciu Sita Eratostenesa. Jest to najszybszy możliwy algorytm generowania liczb pierwszych, którego jedyną wadą jest niemożność sprawdzenia dowolnie wybranej liczby pod kątem pierwszości - musimy wyznaczać je po kolei, 2, 3, 5, 7 i tak dalej. Nie da się za pomocą Sita Eratostenesa stwierdzić, czy liczba 2345832746589327645986324538279453 jest pierwsza, czy nie - ale jeżeli będziemy chcieli znaleźć WSZYSTKIE liczby pierwsze od 2 do 2345832746589327645986324538279453, o ile tylko będziemy dysponować odpowiednio dużą ilością pamięci, Sitem zrobimy to najszybciej.
Idea Sita Eratostenesa jest w gruncie rzeczy dość trywialna. Załóżmy, że chcemy wyznaczyć wszystkie liczby pierwsze od 2 do 2000000. W tym celu wrzucamy do bardzo długiego pudełka klocki ponumerowane od 2 do 2000000, po kolei. Następnie usuwamy wszystkie klocki z numerami podzielnymi przez 2 (z wyjątkiem dwójki). Przesuwamy się do następnego klocka (czyli do trójki) i usuwamy usuwamy wszystkie klocki z numerami podzielnymi przez 3 (oprócz trójki). Potem przechodzimy do kolejnego klocka (piątka, bo czwórkę usunęliśmy w pierwszym kroku) i usuwamy wszystkie klocki podzielne przez pięć, oprócz piątki. Dalej jest siódemka (szóstka wyleciała w pierwszym kroku), potem dopiero jedenastka (ósemka i dziesiątka wyleciały w pierwszym kroku, dziewiątka zaś w drugim). I tak dalej i tak dalej - za każdym razem usuwamy wszystkie wielokrotności bieżącego klocka, z wyjątkiem niego samego.
Po dotarciu do klocka, którego numer jest równy lub większy pierwiastkowi kwadratowemu z dwóch milionów (czyli w naszym przypadku klocka z numerem 1423 - zauważmy, że klocek 1415, który jest najbliżej pierwiastka kwadratowego z 2mln, zostanie usunięty w trzecim kroku) możemy zatrzymać algorytm. Numery na klockach, które pozostały, są liczbami pierwszymi. Wystarczy je posumować, aby dostać wynik zagadki Euler 10. Oczywiście aby uniknąć powtórnego iterowania po tych początkowych klockach, możemy sumę obliczać na bieżąco w tej samej pętli, w której eliminujemy poszczególne klocki.
Zobaczmy jak to wygląda w Pythonie:
# -*- coding: utf-8 -*-
import time
# Zapamiętujemy czas rozpoczęcia wykonywania.
poczatek = time.time()
# Ustalamy górną granicę zgodnie z warunkami zadania
ile_liczb = 2000000
# Najpierw oznaczamy wszystkie liczby jako pierwsze, czyli
# budujemy tablicę 2000001 elementów (liczb od 0 do 2000000 jest 2000001) i wszystkie ustawiamy na True.
# UWAGA: tutaj liczbę reprezentuje indeks elementu tej tablicy, nie jego wartość.
pierwsze = [True] * ile_liczb
# Eliminujemy zero i jedynkę (nie są pierwsze)
pierwsze[0], pierwsze[1] = [False] * 2
# Zerujemy sumę
suma = 0
# Dla każdego elementu tablicy ...
for ind, val in enumerate(pierwsze):
# ...jeżeli jest on równy True ...
if val is True:
# ... to znaczy, że znajduje się on na pozycji będącej liczbą pierwszą,
# a więc dodajemy tę pozycję do wyniku ...
suma += ind
# ... po czym ustawiamy wszystkie elementy na pozycjach będących całkowitą
# wielokrotnością pozycji bieżącego elementu, na False, o ile nie dotarliśmy
# już do pierwiastka z ilości wszystkich liczb, zaokrąglonego w górę
# (w przeciwnym razie wyeliminowaliśmy już wszystkie liczby złożone
# i trzeba tylko posumować to, co zostało)
if ind <= ile_liczb ** 0.5 + 1:
pierwsze[ind * 2 :: ind] = [False] * (((ile_liczb - 1) // ind) - 1)
czas_wykonania = (time.time() - poczatek)
print("wynik: %s, czas wykonania: %s sekund" % (suma, czas_wykonania))
### wynik: 142913828922, czas wykonania: 0.6189999580383301 sekund
### wynik: 142913828922, czas wykonania: 0.3670210838317871 sekund
I tym razem zapuściłem kod na dwóch różnych komputerach - na jednym wykonał się w ciut ponad 0.6s, na drugim w około 0.36s. Jak widać, różnica jest ponad trzydziestokrotna, co dowodzi znacznej przewagi Sita Eratostenesa nad algorytmem naiwnym. Także i w tym przypadku otrzymaliśmy wynik 142 913 828 922 (uff).
Ciekawostką jest, że ten sam algorytm zaimplementowany w językach bardziej niskopoziomowych może działać (na komputerach podobnej klasy) w czasie poniżej jednej dziesiątej sekundy - a więc około osiem-dziesięć razy szybciej, niż w Pythonie. Jednak tu już wychodzimy poza zakres Pchełek (a także niebezpiecznie zbliżamy się do granicy lenistwa autora wpisu), a więc może innym razem 😉
Chciałbym jeszcze dodać kilka komentarzy na temat dzisiejszego kodu. Po pierwsze, jak być może uważny Czytelnik zauważył, na samym początku obydwu skryptów znajduje się tajemnicza linia:
# -*- coding: utf-8 -*-
Jest to odpowiednik HTML-owego meta charset="utf-8", czyli deklaracja strony kodowej, w której będzie napisany skrypt. Oznacza to, że na przykład komentarze z polskimi znaczkami diakrytycznymi ("ogonkami") czy też jakiekolwiek wartości tekstowe zawierające znaki spoza standardowego zakresu ASCII będą poprawnie obsługiwane. Użycie tej opcji nie jest obowiązkowe, ale przydatne, jeżeli chcemy używać znaków specjalnych w kodzie.
Kolejny interesujący trick pojawia się w linii dziesiątej drugiego skryptu:
pierwsze[0], pierwsze[1] = [False] * 2
Tutaj podstawiamy za pierwsze dwa elementy listy wartość False. Po lewej stronie znaku równości mamy dwie osobne wartości listy, z prawej zaś listę jednoelementową [False] - podwojoną. Równie dobrze zamiast [False] * 2 można by napisać [False, False]
Zerknijmy jeszcze na linię 27 drugiego skryptu:
if ind <= ile_liczb ** 0.5 + 1
Tutaj zamiast używać funkcji sqrt() z biblioteki math, podnosimy liczbę do potęgi 1/2. Podwójny znak mnożenia to pythonowy operator podnoszenia do potęgi.
Najwięcej magii odbywa się w następnej linii (28):
pierwsze[ind * 2 :: ind] = [False] * (((ile_liczb - 1) // ind) - 1)
Tutaj po pierwsze widzimy wyrażenie typu L[a:b:c], które oznacza pocięcie L na "plasterki" o grubości c, począwszy od wartości na pozycji a, aż do pozycji (uwaga) b-1. Na przykład:
L = [1, 2, 3, 4, 5, 6] LL = L[0:5:2]
W wyniku wykonania powyższego kodu lista LL będzie równa [1, 3, 5], czyli co drugi element z listy L, od zerowego do czwartego (przypominam, że listy w Pythonie są poindeksowane od zera, a więc piątka w liście L jest na pozycji 4, a nie 5).
Weźmy inny przykład:
L = [1, 2, 3, 4, 5, 6, 7, 8, 9] LL = L[1::3]
Tutaj drugi parametr został pominięty (co oznacza, że bierzemy wszystkie elementy od drugiego aż do końca listy), a więc LL będzie równe [2, 5, 8] (co trzeci element, począwszy od dwójki)
Po powyższych przykładach powinno już być zrozumiałe, co robimy w linii 28: począwszy od elementu ind * 2 (a więc drugiej krotności aktualnie iterowanej liczby pierwszej), aż do końca (::), co ind (czyli kolejne krotności), wszystkie te elementy wypełniamy wartością [False]. Trik polega na tym, że aby podstawić wiele wartości do wielu elementów listy, musimy po obydwu stronach znaku równości mieć dokładnie tę samą liczbę elementów, w przeciwnym bowiem razie dostaniemy błąd. Temu służy wyliczenie ((ile_liczb-1)//ind)-1.
Przy okazji, operator // oznacza dzielenie całkowite, natomiast operator / służy do dzielenia ułamkowego (innymi słowy, 5//2==2, natomiast 5/2==2.5). Tu tkwi mała pułapka, ponieważ zachowanie się operatora / ("zwykłe" dzielenie) uległo zmianie między wersją 2 a 3 Pythona - w wersji 2 wynik dzielenia był albo liczbą całkowitą albo ułamkową, w zależności od tego, czy któryś z operandów był ułamkiem, czy nie. Natomiast w wersji 3 operator / zawsze zwróci liczbę ułamkową (tj. typu float) niezależnie od tego, czy dzielimy przez siebie liczby float, czy int. Dzięki temu działanie kodu jest bardziej przewidywalne i łatwiejsze do kontrolowania.
Wróćmy do naszego Eulera - linia 28 to jest "samo gęste" algorytmu Sita Eratostenesa, tu bowiem "wyciągamy klocki" z naszego długachnego pudełka.
Ktoś mógłby jeszcze postawić pytanie (aczkolwiek wątpię, ponieważ po takiej dawce nudy nawet najwięksi twardziele pewnikiem już drzemią), dlaczego usuwanie wielokrotności wykonujemy tylko do pierwiastka z wartości maksymalnej. Na to pytanie odpowiedź można znaleźć w prezentowanej przeze mnie kiedyś zagadce o stu przełącznikach (proszę sobie poszukać odpowiedniej Pchełki SQL).
Jeżeli chcesz do komentarza wstawić kod, użyj składni:
[code]
tutaj wstaw swój kod
[/code]
Jeżeli zrobisz literówkę lub zmienisz zdanie, możesz edytować komentarz po jego zatwierdzeniu.