Delphi » Gotowce

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.

Gotowce

21 komentarzy

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.

Fermat(N)={ local(a,b2);
	a= sqrtint(N)+1;       // pierwiastek całkowity
	b2= a*a - n;
	while( !issquare(b2),  // czy b2 jest kwadratem
		a++;
		b2 = a*a - N;
	);
	return(a-sqrtint(b2));
}

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:

for i:=3 to N do
    t[i]:= i and 1;      //parzyste oznaczone jako złożone
t[2]:=1;
i:=3; p:=trunc(sqrt(N));
while i<=p do begin
    if t[i]<>0 then begin
        j:=i*i;
        while j<=n do begin
            t[j]:=0;
            j:=j + 2*i 
        end
    end;
    i+=2
end;

Wśród liczb pierwszych jest tylko jedna parzysta liczba (2), więc można by zaoszczędzić połowę pamięci.

for i:=0 to N/2 do
     t[i]:=255;
for i:=1 to trunc(sqrt(N)) div 2 do begin
    if t[i] <> 0 then begin
        j:=(i*i+i)*2;
        while j<=N/2 do begin
            t[j]:=0;
            j+=2*i+1
        end
    end
end;
// jeżeli t[k] > 0 wtedy k*2+1 jest liczbą pierwszą
// jeżeli t[k div 2]>0 wtedy "k" jest liczbą pierwszą (dla nieparzystych "k")

Do określenia czy liczba jest pierwsza, czy złożona wystarczy jeden bit, czyli można by zmniejszyć tablicę ośmiokrotnie.

 //w powyższym programie zmieniamy trzy wiersze
//for i:=0 to N/2 do
for i:=0 to N/16 do

//if t[i]<>0 then begin
if t[i div 8] and (1 shl (i mod 8))<>0 then begin

//t[j]:=0;
t[j div 8]:= t[j div 8] and (255 xor (1 shl (j mod 8)));

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.

var
  p, i, n: longword;   //32 bity
  s, n1: int64;      //64 bity

begin
  p:=2;
  n:=3;  // n  = 2^p - 1
  n1:=2; // n1 = 2^(p-1)
         // n1*n może być liczbą doskonałą
  repeat

    s:=4;                   // 1. Lucas-Lehmer test
    for i:=3 to p do        // 2. 
      s:=(s*s - 2) mod n;   // 3. 
    if (s=0) or (p=2) then  // 4. 

       writeln('Mersenne=2^', p:2, '-1=', n:12, ',  perfect number = ', n1*n);
    p+=1;
    n1:=n+1;
    n:=2*n+1;
  until p>31;
end.

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 15var
n,i,k,t: comp ;
begin

n:= 9000000288000001517;

k:= sqrt(n);
t:= 0;

repeat

i:= n/2; if i*2 = n then break;
i:= n/3; if i*3 = n then break;

t:=5;
while t<=k do begin
  i:= n/t; if i*t = n then break;
  t:= t + 2;
  i:= n/t; if i*t = n then break;
  t:= t + 4;
end

until true;
if t>k then writeln(' is prime') else writeln(' composite');

end.

Pierwiastek z 9000000288000001517 to 3'000'000'047, sprawdzamy 1/3 czyli około miliarda podzielników, ale do 3'000'000'047 włącznie jest tylko 144,449,539 liczb pierwszych, dysponując listą można by ten czas zmniejszyć prawie 7 krotnie. 
Taką listę da się zmieścić w 5 MB ! 

------------------------------------------------------------------------------------------------------------------
Jeszcze jedna ciekawostka. Interesujące, że większość pracy wykonuje operując liczbami 32 bitowymi, unika także dzielenia. Działa dla liczb mniejszych od 4'611'685'936'823'009'641. Jest szybki. Autorem jest chyba Edsger Dijkstra, znalazłem to w jego książce ?Umiejętność programowania?. <code>var
    y,x,n,p:int64;
    k, f,ari, i, j:integer;
    ar:array[0..20] of integer;

var
    b:string;
begin
    //moj kompilator nie lubi stalych calkowitych > 2**31
    b:='1234567890123456789';
    n:=0;
    while b<>'' do begin
        n:=n*10+ord(b[1])-48;
        delete(b,1,1);
    end;

    //szukamy P, najmniejszego pierwszego podzielnika liczby N
    ari:=-1;
    x:=n; y:=2;
    while x<>0 do begin
        ari+=1;        ar[ari]:=x mod y;
        x:=x div y;
        y+=1;
    end;

    f:=2;
    while (ar[0]<>0) and (ari>1) do begin
        f+=1;
        if f>2147483629 then begin
            ar[0]:=1;
            break;
        end;
        i:=0;
        while i<>ari do begin
            j:=i+1;
            k:=ar[i] - j*ar[j];
            while (k<0) do begin
                k+= f + i;
                ar[j]-=1;
            end;
            ar[i]:=k;
            i:=j
        end;
        while ar[ari]=0 do
            ari-=1;
    end;
    if ar[0]=0 then p:=f     //zlozona, P jest podzielnikiem N
    else p:=n;                  //pierwsza
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).

Program paciorki
  var a,b,k,n:longword;
begin
  n:=1283;

  a:=n div 2;
  b:=2;
  k:=n and 1;  //n == a*b+k

  while (a>b) and (k<>0) do begin 
    a -= 1;
    k += b;
    if k >= a then begin
      b += 1;
      k -= a
    end
  end;
  // n == a*b+k
  if k <> 0 then writeln(n,? jest pierwsze?)
end.

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.

If n < 4,759,123,141 is a 2, 7 and 61-SPRP, then n is prime.
____ If n < 1,000,000,000,000 is a 2, 13, 23, and 1662803-SPRP, then n is prime.
Dla 32 bitów starczy więc 2, 7 i 61, a wiesz coś może o "n < 2^64" ?

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

if N<2 then exit; // 0 i 1 nie są pierwsze

if N<4 then begin // 2 i 3 to liczby pierwsze
    result:=true; exit end;

if (n and 1)=0 then exit; //wszystkie parzyste (>=4) są złożone

if n mod 3=0 then exit; // >=4 i podzielne przez 3 są złożone

pierwiastek:=trunc(sqrt(N)); //raz liczymy pierwiastek
x:=5;
skok:=2;
while x<=pierwiastek do begin
    if n mod x = 0 then exit; //jeżeli x dzieli bez reszty N, wtedy N jest złożone
    x+=skok;
    skok:=6-skok
end;
result:=true; //nie znaleźliśmy podzielników N, więc N jest pierwsze

end;

 <b>Sdaj</b>u, sprawdź jakiś milion liczb Twoim i ?moim? sposobem, zauważyłeś różnicę w jakim działają te dwie funkcje?
I tak mniej więcej działa funkcja zaproponowana przez <b>BcbMan</b>?a. Na moje (i nie tylko) oko jest znaaaacznie lepsza. 

Ale to dalej nie wszystko ;-).
Sprawdzamy czy <b>5</b>,7, 11,13, 17,19, 23,<b>25</b>, 29,31... nie są podzielnikami N. Skoro sprawdziliśmy, że 5 nie dzieli N, to czy warto sprawdzać czy 25 jest dzielnikiem? Gdyby podrążyć tę sprawę, to nawet nie matematyk dojdzie do wniosku, że jedynymi podzielnikami, które warto sprawdzać są liczby pierwsze. Więc gdybym wcześniej przygotował sobie tablicę liczb pierwszych i sprawdzał tylko podzielność N przez elementy tej tablicy (mniejsze od pierwiastka z N oczywiście), wtedy jeszcze bym sporo zyskał. Dla 32. bitowego INTEGERa czy raczej CADINALa to tylko pierwiastek(2<sup>32)/ln(pierwiastek(2</sup>32)) liczb, a dokładniej to tylko 6?542 liczb pierwszych mniejszych od pierwiastka z 2^32. Można to w programie policzyć raz, a potem korzystać, można też sobie przygotować tablicę stałych, np.: <code>const
  nprimes=6542;
  primetab:array[0..nprimes] of cardinal = (
      2, // dwa to dziwna liczba ;-)
      3,      5,      7,     11,     13,     17,     19,     23,     29,     31, 
      ........
  65407,  65413,  65419,  65423,  65437,  65447,  65449,  65479,  65497,  65519,
  65521,  65537);

. 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

result:=true;

y:=mpower(a, r, n); // y:=a^r mod n

if (Y<>1) and (Y<>N-1) then

begin
    while (s>1) and (y<>n-1) do begin
        y:=(qword(y)*y) mod n;        // y:=y*y mod n
        if y=1 then exit;
        s-=1
    end;
    if y<>n-1 then exit;
end;
result:=false;

end;

Trochę więcej parametrów, ponieważ zawołamy ją kilka razy dla ustalonego N. Jak zaznaczyłem, wołamy z parametrami R i S takimi, że N-1 == R * 2^S (bardzo łatwo znaleźć R i S, pomyśl o binarnej reprezentacji N-1, . to będzie coś tam (czyli R) z ilomaś tam (S) zerami po prawej stronie.
Parametr A jest pewną liczbą, która być może jest DOWODEM na to, że N jest liczbą ZŁOŻONĄ. Właśnie na tym polega problem, uzyskujemy TYLKO dowód złożoności N (dla bazy A), a nie dowód na to, że N jest pierwsze. Ale dobra wiadomość, wołając to kilkukrotnie tą funkcję, z losowym A (mniejszym od N) prawdopodobieństwo faktu uznania N za liczbę pierwszą dramatycznie maleje. Dramatycznie to znaczy: WYSTARCZAJĄCO. Ale to jeszcze nie wszystko... ;-)
Dla N<=2<sup>32 wystarczy sprawdzić czy A równe 2, 7, 61, nie są dowodami złożoności N. Można to zapisać np. tak:<code> function Rabin_Miller(N:longword):boolean;
var r, s:longword;
begin
    r:=N-1; s:=0;
    while not odd(r) do begin
            r:=r shr 1; s+=1 end;
    //n-1 == 2**s * r
    result:=false;
    if composite(N, r, s,  2) then exit;
    if composite(N, r, s,  7) then exit;
    if composite(N, r, s, 61) then exit;
    result:=true
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 modcode> 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 n
var b:longword;
begin
if r>0 then begin
b:=1;
a:=a mod n;

    while true do begin
        if (r and 1)=1 then
            b:= (qword(b) * a) mod n;

        r:=r shr 1;
        if r = 0 then BREAK;

        a:=(qword(a)*a) mod n;
    end;
    result:=b
end else
    result:=1

end;

Zapisz, sprawdź i co się okazuje? Ano, że nasz wymyślny Rabin z Millerem są jakoś dwa razy słabsi od prostego algorytmu opisanego wcześniej!!! Ale nie pisałem tego bez powodu. 
Okazuje się, że 64 bitowe mnożenie jest głównym sprawcą opóźnień. Zapisanie tego mnożenia (w funkcji mpower() i composite()) w asemblerze (kilkulinijkowe inline) przyśpiesza sprawę kilkudziesięciokrotnie!!! 
Porównajmy trzy algorytmy
1 ? dzielenie N przez wszystkie naturalne <sqrt(N) i nie podzielne przez 2 i 3
2 ? Rabin-Miller bez assemblera
3 ? Inline?owa wstawka mnożąca 
(Sprawdzania aż do N-1 nie bierzemy pod uwagę)
Dla każdego sprawdzę 1?000?000 (milion) liczb poczynając od miliarda, dwu, trzech i czterech:<code>1 ? próbne dzielenie:
        od            do       czas[sek]
           63      1000000    1.10
   1000000001   1001000000   25.71
   2000000001   2001000000   35.37
   3000000001   3001000000   41.77
   4000000001   4001000000   47.51
(zauważcie że czas nie rośnie liniowo z badanym przedziałem)
1a ? sprawdzam podzielność tylko przez wcześniej przygotowaną listę liczb pierwszych
        od            do       czas[sek]
           63      1000000    1.18
   1000000001   1001000000   15.60
   2000000001   2001000000   20.63
   3000000001   3001000000   23.93
   4000000001   4001000000   26.61

2 ? Rabin-Miller 
        od            do       czas[sek]
           63      1000000   38.65
   1000000001   1001000000   60.54
   2000000001   2001000000   62.59
   3000000001   3001000000   62.52
   4000000001   4001000000   64.72
(niby jest źle ;-) ale czas jest prawie stały!!!)
3- Rabin ? Miller (z poprawkami)
        od            do       czas[sek]
           63      1000000    1.51
   1000000001   1001000000    1.88
   2000000001   2001000000    1.91
   3000000001   3001000000    1.94
   4000000001   4001000000    1.96

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.