Liczby pierwsze - szybki algorytm
BcbMan
Poniższy kod znajduje liczby pierwsze, obsługuje liczby Int64, pokazuje czas szukania a jak jakas liczba jest bardzo dlugo szukana to pokazuje ile procent obliczen masz juz za soba. Jak lubicie cwiczyc swoja cierpliwosc to mozecie sprawdzic nawet czy 9000000288000001463 jest liczba pierwsza (ja Wam nie powiem :P - na wynik czekalem prawie 15 min - Celeron 2,4Ghz).
Objaśnienie specjalnie dla Ciebie, Brodny (pozdrowienia i szacunek - mam nadzieje, ze ponizszy opis rozwiaze problem):
wykorzystałem algorytm znajdowania liczb pierwszych (dowod znalazlem w necie, ale nie pamietam, gdzie - wazne ze dziala) zakładajacy, ze jesli sprawdzamy, czy N jest liczba pierwsza to wystarczy sprawdzic podzielnosc N mod i, gdzie i=2, i=3, oraz i=6k-1, i=6k+1 dla k ze zbioru liczb calkowitych dodatnich, np. dla N=1001 wystarczy sprawdzic liczby: 2, 3, 5 (61-1), 7 (61+1), 11 (62-1), 13 (62+1), 17 (63-1), 19 (63+1), 23 (64-1), 25 (64+1), 29 (65-1), 31 (65+1), bo 34>pierwiastek(1001)>33.
Funkcja ND(N) znajduje namniejszy wspolny dzielnik liczby N wiekszy niz 1. Jesli ND(N) = N to N jest liczba pierwsza.
Poza tym nie ma tu nic wiecej procz interfejsu.
Ponizszy kod wklejcie zywcem do nowego projektu dla konsoli w Delphi:
program ExCzyLiczbaPierwsza;
{$APPTYPE CONSOLE}
uses
Math, SysUtils;
function ND(N: Int64): Int64;
var
i: Int64;
begin
if N<2 then ND := 0 // liczby mniejsze niz 2 nie sa pierwsze
else
if N<4 then ND := N // liczby 2 i 3 sa pierwsze
else
// dla liczb wiekszych rownych 4 sprawdzamy najpierw czy sa podzielne
// przez 2 i 3
if (N mod 2=0) then ND := 2
else
if (N mod 3=0) then ND := 3
else
begin
ND := N; i := 1;
// dopiero pozniej sprawdzamy, czy sa podzielne przez 6*i-1 i 6*i+1
// podczas, gdy 6*i-1<=czesci calkowitej z pierwiastka kwadratowego
// badanej liczby
while 6*i-1<=Int(Power(N, 0.5)) do
begin
if N mod(6*i-1)=0 then
begin
ND := 6*i-1; Break;
end
else
if N mod(6*i+1)=0 then
begin
ND := 6*i+1; Break;
end;
// to jest element interfejsu - nie jest niezbedny
// kiedy dlugo trzeba czekac na wynik, to ten kod pokazuje ile
// procent dzielnikow juz sprawdzono
if i mod 1000000=0 then
Write( #13, i div 1000000, 'M (',
100.0*i/Int(Power(N, 0.5)): 0:1, '%)'#13 );
Inc(i);
end
end;
end;
var
tylko_pierwsze, tylko_ilosc, znak : Char;
i, ile_pierwszych, iND, M, N : Int64;
czas, czas_calk : TDateTime;
begin
Writeln( 'Program szuka liczb pierwszych wsrod liczb nieparzystych '
+ 'z podanego zakresu.' );
repeat
tylko_pierwsze := #0;
Write(#13#10#13#10'Podaj poczatek zakresu: '); Readln(M);
Write('Podaj koniec zakresu: '); Readln(N);
Write('Pokaz tylko ilosc liczb pierwszych z tego zakresu [t/n]:');
Readln(tylko_ilosc);
if UpCase(tylko_ilosc)<>'T' then
begin
Write('Pokaz tylko liczby pierwsze [t/n]:');
Readln(tylko_pierwsze);
end;
if M mod 2=0 then M := M + 1;
czas_calk := Now; Writeln; i:=M; ile_pierwszych := 0;
while i<N do
begin
// To ta linijka sprawdza namniejszy dzielnik liczby wiekszy od 1
// oraz oblicza, ile czas zabralo jego szukanie
czas := Now; iND := ND(i); czas := 86400*(Now - czas);
if iND=i then
begin
if UpCase(tylko_ilosc)<>'T' then
Writeln(i, #9'czas: ', czas:0:3, ' s <--- liczba pierwsza ---');
ile_pierwszych := ile_pierwszych + 1;
end
else
if (UpCase(tylko_pierwsze)<>'T')and(UpCase(tylko_ilosc)<>'T') then
Writeln(i, #9'czas: ', czas:0:3, ' s'#9'najmniejszy dzielnik: ', iND);
i := i + 2;
end;
czas_calk := 86400*(Now - czas_calk);
Writeln( #13#10'Ilosc znalezionych liczb pierwszych: ', ile_pierwszych,
#13#10'Calkowity czas wyszukiwania (z wyswietl.): ', czas_calk:0:3,
' s'#13#10'(UWAGA! wyswietlanie bardzo wydluza oczekiwanie na ',
'wyniki)' );
Write(#13#10#13#10'Czy rozpoczac nowe wyszukiwanie [T/N]? '); Readln(znak);
until UpCase(znak)<>'T';
end.
Pozdrawiam, BcbMan.
Kolego Zardax, o probabilistyce też tu było.
Do zadania:
"Czy liczba 9000000288000001463 jest liczbą pierwszą?" można podejść jak do zagadki.
Jeżeli tak długo się to liczy, to pewnie jest iloczynem dwu wielkich liczb pierwszych. I to bliskich sobie.
Wtedy bardzo dobrze działa metoda Fermata.
Kod w PARI/GP ale chyba zrozumiały dla wszystkich.
Dla podanej liczby pętla WHILE wykonuje się ZERO razy.
Gdzie tam 2 minuty, nawet nie wiem jak zmierzyć tak mały czas .
Mamy nie tylko dowód złożoności ale także czynniki.
-----------------------------------------------------------------------------------------------------------[ Xitami ]-
otoż panowie trzeba się trochę zgłębić w podstawowe twierdzenia matematyczne aby ustalic w 2 minuty ze podana liczba jest zlozona mozna ja rozlozyc na iloczyn 2 liczb pierwszych a mianowicie 9000000288000001463=3000000019·3000000077 czy ktos slyszal o testach probabilistycznych
Zaczęło się od tego, że BcBMan zapytał o to czy 9.000.000.288.000.001.463 jest liczbą pierwszą. Sprawdzenie tego zajęło na początku 15 minut.
Liczba operacji wykonywanych przez tą procedurę proporcjonalna jest do pierwiastka z badanej liczby. Pierwiastek z tej liczby to około 3e9, pierwiastek z liczby o połowę mniejszej to 2e9, z proporcji wychodzi, że taką liczbę sprawdzałoby się 10 minut.
W 150 wierszach czystego Pascala, po drodze licząc liczby Fibonacciego i Lucasa, mnożąc macierze, podnosząc do zawrotnych potęg (równych badanej liczbie), liczby o połowę mniejsze testuję w tempie ponad 16'000 na sekundę.
Mniejsze, bo mam kłopot z kilkulinijkową procedurą. Wydaj mi się, że mogę badać liczby nie większe niż około 6.522.000.000.000.000.001, mniejsze od 2^62 z pewnością tak (4,6E18).
Przyśpieszenie o ponad 9 milionów razy :).
Może kiedyś zbiorę to wszystko co tu powypisywałem i ułożę z tego osobny artykuł.
-----------------------------------------------------------------------------------------------------------[ Xitami ]-
Ważnym algorytmem związanym z liczbami pierwszymi jest sito Eratostenesa.
Nieco zoptymalizowana wersja:
Wśród liczb pierwszych jest tylko jedna parzysta liczba (2), więc można by zaoszczędzić połowę pamięci.
Do określenia czy liczba jest pierwsza, czy złożona wystarczy jeden bit, czyli można by zmniejszyć tablicę ośmiokrotnie.
W troszkę ponad 2 minuty znajduje wszystkie liczby pierwsze mniejsze od miliarda, tablica waży 62.5 MB. Odpowiednie ustawienie opcji kompilatora zmniejsza czas o połowę.
-----------------------------------------------------------------------------------------------------------[ Xitami ]-
Sam test jest stosunkowo prosty. W załączonym Pascal’owym programie zaznaczyłem 4 wiersze tego testu. Sprawa jest jak widać prosta kiedy posługujemy się „wbudowanymi” typami danych. Rzecz komplikuje się gdy „n*n” jest większe niż dostępny, „wbudowany” typ (wiersz //3 testu).
Załączam programik drukujący kilka liczb pierwszych Mersenne’a i odpowiadających im liczb doskonałych.
Przypomnę, że liczba Mersenne’a jest liczbą o postaci 2^p - 1.
Liczba 2^p - 1 może bycz liczbą pierwszą (Mersenne’a) wtety i tylko wtedy gdy "p" jest liczbą pierwszą.
Przed przystąpieniem do testu Lucasa-Lehmera warto sprawdzić pierwszość wykładnika. Dla małych liczb (jak tu) ten krok nie ma sensu.
-----------------------------------------------------------------------------------------------------------[ Xitami ]-
Chciałbym się jeszcze zapytać, czy ktoś z Was wie jak szuka się bardzo wielkich liczb pierwszych, np. "44 liczba Mersenne'a: 2^32582657?1 i liczy sobie 9808358 cyfr w zapisie dziesiętnym." - cytat z Wikipedii.
Wujek.wroc podsunął pomysła, dla liczb <=9'223'372'036'854'775'807 lepszym okazuje się typ "comp", zmiennoprzecinkowe 63 bity plus znak (choć przyjmuje tylko wartości całkowite).
Liczby z przedziału 232 <= N < 264 warto sprawdzać inaczej.
Czas wykonania tego przykładu to 127 sekund, czyli 2 min zamiast 15
var
n,i,k,t: comp ;
begin
n:= 9000000288000001517;
k:= sqrt(n);
t:= 0;
repeat
until true;
if t>k then writeln(' is prime') else writeln(' composite');
end.
62 bity, góra 2 minuty
-----------------------------------------------------------------------------------------------------[ Xitami ]---
No, wreszcie odezwał się ktoś kompetentny. Dzięki Xitami za komentarz i sensowną krytykę. W wolnej chwili poprawię powyższy kod.
Jeśli chodzi o tworzenie tablicy to ma ono sens dla względnie małych liczb. Zawsze znajdzie się taka liczba, dla której tworząc tablicę zabraknie pamięci, opisany przeze mnie sposób jest uniwersalny przez co nie jest rozwiązaniem optymalnym dla większości problemów (ale jest jednym z najlepszych). Jedno jest pewne: algorytmy w innych artykułach stworzonych przed moim są mniej efektywne. I stąd mój artykuł.
Jeśli ktoś stworzył lepszy algorytm od mojego to niech dołączy link do swojego artykułu, a wszyscy będziemy mu wdzięczni :)
bardziej przydałaby się funkcja sprawdzająca, czy liczba jest pierwsza, np. ta, którą sobie dzisiaj napisałem w Pascalu:
function pierwsza(x: integer):boolean;
var b: integer;
f: boolean;
begin
if (x=0) or (x=1) then f:=false
else if x=2 then f:=true
else
begin
for b:=2 to x-1 do
begin
if x mod b>0 then
begin
f:=true;
end
else
begin
f:=false;
Break;
end;
end;
end;
pierwsza:=f;
end;
A jako ciekawostkę dodam najwolniejszy algorytm testujący pierwszość (i szukający czynników).
Wujek.wroc a jakiego typu u Ciebie jest liczba: 9000000288000001463?
Od kiedy ?nieco ponad? jest szybsze od ?poniżej??
----------------------------------------------------------------------------------------- Xitami
Witam,
sprawdziłem, że dobrze napisany program szukający liczb pierwszych wg algorytmu naiwnego jest szybszy od Twojego. Pozbywając się ProgressBar-ów i innych upiększeń sprawdziłem czy liczba 9000000288000001463 jest pierwsza w nieco ponad 2 min. na Duronie 1000. Korzystając z pomysłu przytoczonego w Twoim artykule zszedłem poniżej 2 minut. Sprawdziłem, że funkcja mod jest nieefektywna i zamiast tego liczę w sposób następujący:
Dzielna := Round(Liczba / Dzielnik);
if Dzielnik * Dzielna = Liczba then Pierwsza:=False;
Znacznie skróciło to obliczenia.
BcbMan:
Sporo zależy od zakresu jaki Cię interesuje.
Przesadziłeś trochę z szacowaniem wielkości koniecznej pamięci. Żeby sprawdzić czy podana przez ciebie liczba jest pierwsza, potrzebuję listy wszystkich liczb pierwszych nie większych od jej pierwiastka, jest ich około 3E9/ln(3E9) = 137?476?711, razy 4 bajty (longword) raptem 550 MB, to nie tak wiele.
Może nawet nie trzeba 4 bajtów, wystarczy jeden z różnicą między kolejnymi.
Dryobates:
____ If n < 1,373,653 is a both 2 and 3-SPRP, then n is prime.
____ If n < 25,326,001 is a 2, 3 and 5-SPRP, then n is prime.
____ If n < 25,000,000,000 is a 2, 3, 5 and 7-SPRP, then either n = 3,215,031,751 or n is prime. (This is actually true for n < 118,670,087,467.)
____ If n < 2,152,302,898,747 is a 2, 3, 5, 7 and 11-SPRP, then n is prime.
____ If n < 3,474,749,660,383 is a 2, 3, 5, 7, 11 and 13-SPRP, then n is prime.
____ If n < 341,550,071,728,321 is a 2, 3, 5, 7, 11, 13 and 17-SPRP, then n is prime.
inaczej:
____ If n < 9,080,191 is a both 31 and 73-SPRP, then n is prime.
Mam fajnego gotowca, poczynając od 9000000288000001457 sprawdziłem 4400 kolejnych (nieparzystych) liczb, zajęło to jedną sekundę.
Napisane w assemblerze (dla FreePascal może nada się i do Delphi)
http://www.polarhome.com/~franco/
----------------------------------------------------------------------------------------- Xitami
Drogi TeWuX-ie, sprawdź za pomocą tego programu, czy np. liczba 9000000288000001457 jest liczbą pierwszą (powinno się udać za xxx lat - w miejsce x wstaw dowolną cyfrę większą od 0 - jak komputery będą miały terabajty terabajtów pamięci RAM :)).
A poważnie: taka tablica zajmuje zbyt wiele miejsca w pamięci i dlatego nie można jej wykorzystać przy większych liczbach. Trzeba szukać pewnych reguł, które pozwolą na najmnijeszą liczbę sprawdzeń na podzielność. Twój sposób jest najlepszy, ale tylko dla względnie małych liczb, tzn. dla takich dla których tworzona przez Ciebie tablica zmieści się w pamięci RAM, a wierz mi, że to miejsce bardzo szybko się zapełni niezależnie czy wsadzisz tam 256 MB, czy 2 GB (czyli max ile można wsadzić do 32-bitowej maszyny).
I muszę Cię rozczarować: nawet przesiadka na maszynę 64-bitową nie pozwoli na sprawdzenie w/w liczby. Dlatego szukam lepszego sposobu.
Ja tam i tak najbardziej lubię to:
Kapuściane materiały
dziala na nieco podobnej zasadzie co algoryt opisany przeze mnie:
Wyszukiwanie liczb pierwszych
Tylko, że mój zapisuje liczby pierwsze które znalazł do tablicy, i kolejną liczbe sprawdza porównując własnie z liczbami z tej tablicy.
Mój algolrym, potrzebuje mniej pamieci od Sita, a twój tyle co naiwny algorytm (czyli niewiele ;) ), przyczym jest pewnie szybszy od naiwnego.
Sdaju, gdybyś lepiej przyjrzał się temu, co napisał BcbMan to zobaczyłbyś, że funkcja ND() bada czy parametr z którym jest wołana jest liczbą pierwszą. Dokładniej, szuka podzielnika, ale robi to dość niekonsekwentnie, np. dla jeden, odpowiada, że zero jest podzielnikiem ;-). BcbMan zrobił to znacznie zgrabniej niż Ty.
Prostym (choć nieefektywnym) sposobem testowania pierwszości danej liczby jest szukanie podzielnika, czyli liczby która dzieli daną bez reszty. Można próbować dzielić daną liczbę ?N? przez wszystkie możliwe dzielniki, czyli liczby od dwa do N-1. Ale czy warto? Jeżeli testujemy akurat liczbę pierwszą to będziemy sprawdzać jej podzielność przez N-2 liczb.
Czy N-1 może być podzielnikiem N? Z pewnością nie. Największą liczbą, która być może jest dzielnikiem N jest N/2. A Ty szukasz aż do N-1 (pętla: for b:=2 to x-1 do). Więc można by poprawić: for b:=2 to x div 2 do. Jeżeli akurat sprawdzana liczba jest pierwsza, funkcja powie nam o tym w o połowę krótszym czasie.
Ale to nie wszystko. Czy musimy wypróbowywać liczby większe od pierwiastka z N? Jeżeli N=pq, oraz p<pierwiastek(N) wtedy q musi być większe od pierwiastka, więc - nie warto sprawdzać czy liczby większe od pierwiastka z N są podzielnikami N. Gdyby tak było to pętla już wcześniej by się zakończyła (BREAK), ale jeżeli testujemy akurat liczbę pierwszą to znowu będziemy sprawdzać coś, co już wcześniej było sprawdzone. Poprawmy, zamiast testować aż do N/2 sprawdzajmy tylko do pierwiastka z N. Pętla zamiast kręcić się jałowo miliard razy, zrobi tylko jakieś 65 tysięcy obrotów, a to chyba spora różnica. I to nie wszystko. Jeżeli sprawdziłem, że N nie dzieli się przez 2 (N jest liczbą nieparzystą) to czy warto sprawdzać czy dzieli się przez 4, 6, 8, itd.? No chyba nie warto. Więc zamiast pętli FOR zapiszmy to jako WHILE, a potencjalny dzielnik (X) niech wystartuje od 3 i zmienia się nie o jeden, ale o dwa. Znowu przyśpieszyliśmy dwukrotnie.
Kiedy ma zakończyć się pętla WHILE? No wiemy, wtedy gdy X przekroczy pierwiastek z N. Ale czy warto z każdym obrotem liczyć pierwiastek? Czy kompilator nam pomoże i tylko raz policzy pierwiastek? Nie koniecznie, pomóżmy mu jakąś pomocniczą zmienną. Raz policzmy pierwiastek z N, i niech WHILE kręci się aż do X>=PIERWIASTEK (BcbMan: popraw to). Algorytm zaproponowany przez BcbMan?a jest jeszcze lepszy, bo skoro sprawdziliśmy, że N nie dzieli się bez reszty przez 3 to nie warto sprawdzać czy aby nie dzieli się przez 6, 9, itd. BcbMan dość często liczy 6I-1, można zrobić to inaczej. Zmienna X nie musi przyrastać o dwa, ale na zmianę, o 2 albo o 4. Np. tak: przed pętlą: SKOK:=2; X:=3; a w pętli X:=X+SKOK; SKOK:=6-SKOK;. Jeżeli zaczniemy od X:=5, podzielność N będziemy testować tylko dla: 5,7, 11,13, 17,19, 23,25, 29,31... .
BcbMan?ie, mnie jakoś bardziej podoba się ten sposób, niż mnożenie przez 6 i odejmowanie.
Można by tak ulepszoną funkcję zapisać np. tak:
function pierwsza(n:cardinal):boolean;
var x, pierwiastek, skok:cardinal;
begin
result:= false; //zakładamy że N jest złożona
end;
. Taki fragment kodu można sobie oczywiście wygenerować maleńkim programikiem, albo poprosić mnie, to podeślę ;-).
Czyli dla 32. bitowych liczb, zamiast ewentualnie testować podzielność przez ponad 4 miliardy liczb, sprawdzamy tylko troszkę ponad 6 tysięcy. A to jest spora różnica!!! Prawie milion razy!!!
Dla liczb 64 bitowych to już sposób trochę kłopotliwy, tablica 1e17 liczb 32 bitowych, więc trochę dużo ;-).
Na moje oko tyle sam powinien wymyślić ?komputerowiec?. Bo ANALITYK to już zupełnie inna historia.
Pomysłowość (i złośliwość ;-) ) ludzka nie zna prawie granic.
Pokażę pewną implementację testu Rabin-Miller?a. Testu, który z szansą lepszą niż ofiaruje moneta, odpowie nam na pytanie czy N jest liczbą pierwszą?
Pozostanę przy 32 bitach, no nie do końca, ponieważ istotnym elementem jest obliczenie AB MOD N, gdzie liczby A, B i N są trzydziesto dwu bitowe, uwaga: pośredni wynik (AB) jest 64 bitowy. Dzięki promocji np. wartości ?A? do 64 bitów (typ QWORD we FreePascalu, long long unsigned w C) jesteśmy w stanie to policzyć, np. tak: W:= (QWORD(A)*B) mod N;.
Funkcja composite (z wielki prawdopodobieństwem) odpowie nam czy N jest liczbą złożoną. To prawdopodobieństwo da się zamienić w pewność, o tem potem ;-)
function composite(N, r, s, a:longword):boolean; //N-1 == r * 2^s
var y: longword;
begin
end;
Wcześniej pojawiła się funkcja ?modularnego potęgowania?. Temat sam w sobie ciekawy. Jak obliczyć A</sup>R mod N? Można np. R-razy liczyć B:=BA mod N (na początek B:=A). Tu pomoże matematyk, bo powie że takie ?R?-krotne mnożenie da ten sam wynik co bezpośrednie policzenie AR i wzięcie tego modulo. To pierwszy krok, ale R może być sporą liczbą (N-1 == R * 2s) powiesz. I na to jest sposób ? nie potrzebujemy R iteracji, wystarczy tylko tyle ile R ma bitów (czyli log2(R)):
function mpower(a, r, n:longword):longword; // a**r mod
code> Wcześniej pojawiła się funkcja ?modularnego potęgowania?. Temat sam w sobie ciekawy. Jak obliczyć A</sup>R mod N? Można np. R-razy liczyć B:=BA mod N (na początek B:=A). Tu pomoże matematyk, bo powie że takie ?R?-krotne mnożenie da ten sam wynik co bezpośrednie policzenie AR i wzięcie tego modulo. To pierwszy krok, ale R może być sporą liczbą (N-1 == R * 2s) powiesz. I na to jest sposób ? nie potrzebujemy R iteracji, wystarczy tylko tyle ile R ma bitów (czyli log2(R)): ` function mpower(a, r, n:longword):longword; // a**r mod nvar b:longword;
begin
if r>0 then begin
b:=1;
a:=a mod n;
end;
Sprawdzałem tylko liczby nieparzyste.
Poprawki w 3 algorytmie to nie tylko zapisanie modularnego mnożenia w asemblerze, to nie wszystko, troszkę poprawia sprawę wstępne, próbne dzielenie N przez kilka początkowych liczb pierwszych.
No to co Sdaju, wyszło, że mój słaby komputer i kompilator jest jednak tysiące razy szybszy od twojego ;-)
---------------------------------------------------------------------------------------------[ Xitami ]---
cakiem "maly" algorytm ale gratuluje ze ci sie go udalo napisac!!!!
Jest i opis...
Artykuł i żadnego wyjaśnienia? Jak ktoś będzie szukał opisu na sieci to znajdzie i implementację.
toż to nic innego, jak namiastka sita Eratostenesa - sprawdzanie podzielności przez wszystkie liczby pierwsze (jak się okazuje w tym algorytmie nie tylko) mniejsze od pierwiastka ze sprawdzanej liczby.
jest jednak "drobna" roznica. Ile bedziesz potrzebowal pamieci RAM, zeby utworzyc tablice liczb od 2 do 9000000288000001463 ?? !! :P (co najmniej 8*9000000288000001463 B = ok. 8.185.452,6 terabajtów pamieci RAM) nie mowiac juz ile czasu zajmie wyluskanie liczb pierwszych, podczas gdy ja chce wiedziec, czy jedna z nich jest pierwsza...
powyzszy program zajmuje w pamieci RAM ok. 100 B, ale przede wszystkim doczekasz sie na wynik (ja czekalem ok. 15 minut).
I co Ty na to? (mimo wszystko dalej licze na to, ze ktos ulepszy jeszcze ten algorytm)
acha, a propos Twojego komentarza "jak się okazuje w tym algorytmie nie tylko" to lepiej przeanalizuj kod, bo sie mylisz. Poza tym algorytm sita zaklada ze sprawdzasz wszystkie liczby od 1 do N, wiec ten algorytm jest i tak o niebo lepszy.