Rzecz będzie o matematyce. Proszę się nie martwić, nie zamierzam dziś udowadniać, że P = NP ani hipotezy Riemanna. Nie zamierzam nawet tłumaczyć szczegółów tytułowego testu AKS. Wpomnę tylko o nim, a następnie rzucimy się w wir kodowania.
Test AKS jest testem pierwszości liczb. O testach pierwszości już wcześniej pisałem, o tutaj - wtedy jednak było o testach probabilistycznych (a więc takich, które pozwalają stwierdzić, czy dana liczba jest liczbą pierwszą, z dowolnym - jednak mniejszym od jedności - prawdopodobieństwem). Natomiast test AKS jest najlepszym obecnie znanym testem pierwszości, który wynik daje ze stuprocentową pewnością. W dodatku jest on jedynym znanym obecnie testem pierwszości, którego złożoność jest wielomianowa względem ilości cyfr sprawdzanej liczby (a więc, w odróżnieniu od innych testów, nie jest tak, że czas potrzebny na sprawdzenie pierwszości rośnie wykładniczo wraz ze wzrostem ilości cyfr liczby testowanej).
Sam test jest algorytmicznie dość skomplikowany. Poprzednie zdanie stoi pozornie w sprzeczności z końcówką poprzedniego akapitu - pisałem przecież, że złożoność testu nie jest wykładnicza, a zaraz potem piszę, że jest skomplikowany. Otóż chodzi o to, że wyjaśnienie zasady działania testu AKS jest raczej niestrawne dla kogoś, kto nie pasjonuje się od dzieciństwa matematyką, a hipotezy Fermata nie wyssał z mlekiem matki. Z drugiej jednak strony, jak się już ten test zaimplementuje w postaci algorytmu, czas potrzebny na jego zakończenie dla dowolnie dużej (jednak skończonej) liczby jest krótszy niż w przypadku innych algorytmów.
Wyjaśniwszy ten pozorny dualizm semantyczny, spróbuję wreszcie przejść do sedna.
Bez wnikania w szczegóły implementacji testu AKS, najważniejszym jego elementem jest to, że opiera się on na kongruencji liczb (x - a)^n oraz (x^n - a) (względem n) dla dowolnych naturalnych x oraz a tylko wtedy, jeżeli n jest liczbą pierwszą.
Co to znaczy, że dwie liczby są kongruentne względem jakiejś liczby N, zapyta Czujny Czytelnik?
Znaczy to, że obydwie dają tę samą resztę przy dzieleniu przez N. Na przykład 5 i 8 są kongruentne względem 3, bo reszta z dzielenia 5 przez 3 jest 2, tak samo jak reszta z dzielenia 8 przez 3.
No więc test AKS opera się na tej właśnie kongruencji. Nawiasem mówiąc, kongruencja ta jest uogólnieniem małego twierdzenia Fermata (zainteresowanych szczegółami odsyłam na Wikipedię), co akurat ma niewielkie znaczenie dla niniejszego wpisu.
Dzisiejsza pchełka skupi się na iteracyjnym sprawdzeniu tej kongruencji, dla kilku niedużych wartości A, N oraz X. Od razu mówię, że wbrew pozorom, ilość kombinacji do sprawdzenia jest całkiem pokaźna, nawet "tylko" w zakresie liczb 32-bitowych (większe liczby nie wchodzą w rachubę, ponieważ VBA sobie z nimi nie za bardzo radzi).
Powtórzę raz jeszcze, co będziemy sprawdzać: będziemy generować trójki liczb naturalnych A, N oraz X, takie, że A < X oraz N jest pierwsze, a następnie będziemy porównywać reszty z dzielenia (X^N - A) przez N oraz (X - A)^N przez N. Za każdym razem obydwie te reszty powinny być takie same.
Zaczniemy - jak zawsze - od pełnego kodu. A więc, Alt-F11, prawomysz na drzewku, wstaw nowy moduł, piszemy:
Option Explicit
Private Function DoPotegi(ByVal x As Long, ByVal n As Long) As Double
Dim i As Long, retval As Double
retval = 1
For i = 1 To n
retval = retval * x
Next i
DoPotegi = retval
End Function
Private Function NtaPierwsza(n As Long) As Long
NtaPierwsza = Switch(n = 1, 2, n = 2, 3, n = 3, 5, n = 4, 7, n = 5, 11, n = 6, 13, n = 7, 17, n = 8, 19, n = 9, 23, n = 10, 29, n = 11, 31, n = 12, 37, n = 13, 41, n = 14, 43, n = 15, 47, True, -1)
End Function
Sub Pchelka(ByVal x As Long, ByVal np As Long)
Dim cx As Long, ca As Long, cp As Long, crow As Long
crow = 1
Range("A:B").Clear
For cx = 2 To x
For ca = 1 To cx - 1
For cp = 1 To np
Range("A1").Offset(crow).Value = "sprawdzamy czy (" & cx & " - " & ca & ")^" & NtaPierwsza(cp) & " mod " & NtaPierwsza(cp) & " = (" & cx & "^" & NtaPierwsza(cp) & " - " & ca & ") mod " & NtaPierwsza(cp)
If DoPotegi(cx - ca, NtaPierwsza(cp)) Mod NtaPierwsza(cp) = (DoPotegi(cx, NtaPierwsza(cp)) - ca) Mod NtaPierwsza(cp) Then Range("B1").Offset(crow).Value = "TAK" Else Range("B1").Offset(crow).Value = "NIE"
crow = crow + 1
Next cp
Next ca
Next cx
End Sub
Teraz zrobimy rozbiór powyższego kodu na czynniki (nomen-omen) pierwsze, a więc kilka słów wyjaśnienia dla tych, którzy jeszcze nie przysnęli.
Option Explicit - ta linijka powinna znaleźć się na początku każdego modułu - jak już kilkakrotnie tłumaczyłem, wymusza ona na programiście deklarowanie wszystkich zmiennych użytych w kodzie, dzięki czemu łatwiej unikniemy literówek. Taki bacik na samego siebie, jednak bardzo przydatny.
Private Function DoPotegi(ByVal x As Long, ByVal n As Long) As Double
Dim i As Long, retval As Double
retval = 1
For i = 1 To n
retval = retval * x
Next i
DoPotegi = retval
End Function
Powyższa funkcja będzie nam podnosić x do potęgi n. VBA niestety nie posiada takiej funkcji wbudowanej, dlatego musimy ją sobie sami napisać. Ponieważ wszystkie elementy tej funkcji były już wcześniej omówione (proszę sobie zajrzeć do poprzednich pchełek), nie będę jej tutaj omawiał szczegółowo.
Private Function NtaPierwsza(n As Long) As Long
NtaPierwsza = Switch(n = 1, 2, n = 2, 3, n = 3, 5, n = 4, 7, n = 5, 11, n = 6, 13, n = 7, 17, n = 8, 19, n = 9, 23, n = 10, 29, n = 11, 31, n = 12, 37, n = 13, 41, n = 14, 43, n = 15, 47, True, -1)
End Function
Ta z kolei funkcja zwraca nam N-tą z kolei liczbę pierwszą, przy czym działa tylko dla N <=15. Ponieważ sprawdzana kongruencja wymaga podnoszenia do potęgi będącej liczbą pierwszą, siłą rzeczy przy dużych wykładnikach wyjdziemy poza zakres liczb 32-bitowych generując błąd, więc nie ma sensu zwracać liczb pierwszych większych niż 47. Funkcja Switch była już omawiana wcześniej, przypomnę tylko, że sprawdza ona kolejno warunki podawane w parametrach na pozycjach nieparzystych i zwraca wartość parametru na pozycji parzystej zaraz za znalezionym warunkiem prawdziwym. Jeżeli więc n będzie równe 4, funkcja Switch sprawdzi, że siódmy parametr jest prawdziwy, i zwróci wartość ósmego parametru (czyli siódemkę, która jest czwartą z kolei liczbą pierwszą).
Sub Pchelka(ByVal x As Long, ByVal np As Long)
Dim cx As Long, ca As Long, cp As Long, crow As Long
crow = 1
Range("A:B").Clear
For cx = 2 To x
For ca = 1 To cx - 1
For cp = 1 To np
Range("A1").Offset(crow).Value = "sprawdzamy czy (" & cx & " - " & ca & ")^" & NtaPierwsza(cp) & " mod " & NtaPierwsza(cp) & " = (" & cx & "^" & NtaPierwsza(cp) & " - " & ca & ") mod " & NtaPierwsza(cp)
If DoPotegi(cx - ca, NtaPierwsza(cp)) Mod NtaPierwsza(cp) = (DoPotegi(cx, NtaPierwsza(cp)) - ca) Mod NtaPierwsza(cp) Then Range("B1").Offset(crow).Value = "TAK" Else Range("B1").Offset(crow).Value = "NIE"
crow = crow + 1
Next cp
Next ca
Next cx
End Sub
Tutaj konsumujemy samo gęste, czyli faktyczne testowanie kongruencji dla zadanego zakresu maksymalnych wartości A, N oraz X. Procedurę nazwałem przewrotnie Pchelka, przyjmuje ona dwa parametry: maksymalną wartość X oraz maksymalny numer kolejny liczby pierwszej (a więc dla np = 3 kongruencja będzie sprawdzona dla trzech kolejnych liczb pierwszych: 2, 3 oraz 5 i tak dalej).
Najpierw definiujemy sobie cztery zmienne numeryczne: cx (bieżąca wartość X), ca (bieżąca wartość a), cp (bieżąca wartość indeksu kolejnej liczby pierwszej) oraz crow - wirtualny "kursor", czyli numer bieżącego wiersza w arkuszu, do którego będziemy wstawiać wynik sprawdzenia.
Następnie ustawiamy nasz "kursor" na pierwszym wierszu a potem zagnieżdżamy kolejno trzy pętle. Pierwsza, najbardziej zewnętrzna pętla, ogranicza wartość zmiennej X, w zakresie od 2 do cx. A więc pierwsza iteracja będzie dla dwójki, druga dla trójki i tak dalej aż do ostatniej iteracji dla x=cx.
Druga, wewnętrzna pętla, steruje zmienną a, która przyjmuje wartości kolejno od 1 do cx - 1. Proszę zauważyć, że górna granica jest zmienna i zależy od wartości cx (z zewnętrznej pętli). A więc ilość iteracji drugiej pętli będzie rosła o jeden z każdą iteracją zewnętrznej (pierwszej) pętli.
Wreszcie trzecia, najbardziej wewnętrzna pętla, iteruje po kolejnych liczbach pierwszych, począwszy od pierwszej aż do NP-tej. Ta pętla zawsze wykona się tą samą (NP) ilość razy.
Ciało pętli stanowią trzy instrukcje: najpierw, w kolumnie A, wpisujemy informację co w danej iteracji jest sprawdzane. Następnie, po słówku If, sprawdzamy, czy kongruencja zachodzi dla bieżącej trójki liczb ca, cx oraz cn. Wreszcie, po słówku Then, wpisujemy w kolumnie B słowo "TAK" jeżeli reszty z dzielenia były równe, bądź "NIE", jeżeli nie były.
Proszę sobie teraz wykonać procedurę Pchelka, dla X=20 oraz NP=4:
Ctrl-G
Pchelka 20, 4 <enter>
Przełączamy się na arkusz i widzimy, że w kolumnie A pojawiło się 760 wpisów o sprawdzaniu kongruencji, a w kolumnie B dla każdego sprawdzenia mamy wynik TAK.
Proszę popróbować z innymi wartościami parametrów - w większości przypadków, zwłaszcza dla NP większych niż 4, dostaniemy błąd przepełnienia, ponieważ podnosimy do potęgi jedenastej, wychodząc tym samym poza zakres 32-bitowy. Jednak celem niniejszej pchełki nie było udowodnienie tej kongruencji, tylko nauczenie się kolejnego kawałka VBA.
Jeszcze odrobina samokrytyki na koniec:
- Kod jest całkowicie niezoptymalizowany pod kątem wydajności: niepotrzebnie wołamy funkcję DoPotegi oraz NtaPierwsza tyle razy wewnątrz pętli,
- Kod nie działa dla liczb większych niż 4294967295 (maksymalna wartość typu Long) - gdyby się uprzeć, można by napisać bibliotekę operującą na większych liczbach, jednak staram się tutaj pozostać w dziedzinie "pchełek" czyli niedużych porcji kodu
- Ilość parametrów wewnątrz funkcji Switch należałoby zredukować do maksymalnie 10-12.
- I tak dalej...
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.