Sieb des Eratosthenes

[Amateurprofi]

Zitat von Bjoerk:

Ich nehm mal an, daß Amateurprofi den Algo von sx2008 benutzt. Hätte mich auch mal interessiert, ob Amateurprofi diesen um Faktor 2 beschleunigt hat, denn das Sieb braucht die geraden Zahlen ja nicht.

Hast du mal beide verglichen?
Lass mal das von sx2008 laufen.
Und dann das von mir.

[p80286]

Ich hab hier ein altes TP-Sieb ausgegraben:

Delphi-Quellcode:

program Project2;

{$APPTYPE CONSOLE}

uses
  SysUtils;

const
  LOOPEND= $FFFFFFF;
var
  sieb : array [0..LOOPEND] of byte;
  i,j,lauf : longint;
  start,
  ende : tdatetime;
begin
  writeln (' Sieber');
  fillchar(sieb,sizeof(sieb),#0);
  start:=time;
  i:=2;
  while i < (LOOPEND shr 1) do begin
    j:=i *2;
    while j <= LOOPEND do begin;
      sieb[j] := 1;
      inc(j,i);
    end;
    repeat
      inc(i,1);
    until sieb[i]=0;
  end;
  ende:=time;
  j:=1;
  for i:=1 to loopend do if sieb[i]=0 then inc(j,1);
  writeln(' im Bereich von 0 bis ',formatfloat('0,',loopend/1),' habe ich ',j,' Primzahlen gefunden!');
 writeln('Ben”tigte Zeit:',formatdatetime('HH.MM:SS:HH:ZZZ',ende-start));
  readln;
  for i:=1 to 20 {LOOPEND} do begin
    if sieb[i]=0 then write(i:8);
  end;
  readln;
end.

Die reine Rechenzeit für ca 250 Mio liegt bei ~15-16sec

Gruß
K-H

[Amateurprofi]

Für 250 Mio braucht meins (bei mir - Intel I7) etwa 600ms, und das ist nicht nur die Rechenzeit.
Ich schrieb ja, das es etwas bummelig ist und dass ich mit dem Zeitverhalten nicht so ganz zufrieden bin.

[p80286]

Unbenommen, wenn man den Aufwand für die Oberfläche mitrechnet ist Deine Version bestimmt schneller, aber meine Intention war erst einmal eine einfache Implementierung zur Verfügung zu stellen.
Hier wäre der nächste Schritt bestimmt, die geraden Zahlen(2) aus der Berechnung zu eliminieren.

Gruß
K-H

[Horst_]

Hallo,

Was mich mal speziell interessieren würde, ob jemand die gleichen Probleme hat mit einem AMD Rechner hat.

"damals" gabe es doch mal eine parallele Diskussion zu Primzahlsieben:
Delphi-Forum: Optimierung einer Primzahlfunktion und hier: sehr-schneller-primzahl-finder
Zum Ende der Threads werden die Programme immer schneller und aufwendiger
Wenn man nicht alle Zahlen auf einmal siebt, sondern immer nur Cache-freundliche Abschnitte, wird die Sache richtig schnell.
Kim Walisch siebt, glaube ich, mit zwei Zahlen gleichzeitig, weil die Daten in der gleichen Cache-line sind. Dessen Quellcode ist ja offen, aber C++
http://code.google.com/p/primesieve/

Zitat:

primesieve generates the first 50,847,534 primes up to 10^9 in just 0.4 seconds on a single core of an Intel Core i7-920 2.66GHz

Ich habe mal eine Quelloffene Variante angehängt, die abschnittsweise siebt.
Es hält ja niemand einen davon ab, nach jedem Siebdurchgang, das Sieb abzuspeichern.
4 Sekunden bis Obergrenze = 1 shl 32 -1, also etwa 250% langsamer in dem Bereich als Kim Walisch Version. ( 1 Sekunde bis 1e9)

Gruß Horst

[Morphie]

Habe auch einen AMD und bekomme ebenfalls beim Starten die Fehlermeldung.
edit: Windows 8 Pro x64

[Amateurprofi]

Zitat von Morphie:

Habe auch einen AMD und bekomme ebenfalls beim Starten die Fehlermeldung.
edit: Windows 8 Pro x64

Ich habe meine Vermutungen, wo das Problem liegt, nämlich dass ich irgendwo eine ASM-Instruktion
benutze, die von den AMD-Prozessoren nicht unterstützt wird.

Ich habe Horst die Source-Codes zugeschickt, mit der Bitte, die in Frage kommenden Teile Mal mit F7 durchzuklappern.

[Amateurprofi]

Zitat von p80286:

Unbenommen, wenn man den Aufwand für die Oberfläche mitrechnet ist Deine Version bestimmt schneller, aber meine Intention war erst einmal eine einfache Implementierung zur Verfügung zu stellen.
Hier wäre der nächste Schritt bestimmt, die geraden Zahlen(2) aus der Berechnung zu eliminieren.

Gruß
K-H

Nee, die ist nicht nur dann schneller, wenn man den Aufwand für die Oberfläche mitrechnet.
Du schriebst 15-16 für 250M, ich maß 600ms für 250M also um den Faktor 25 schneller.
Aber das nur nebenbei - meine Version ist recht unoptimiert, und deshalb für meine Begriffe bummelig.
Horst hat ja einen Link reingestellt auf eine wirklich schnelle Version.

Und der nächste Schritt, die geraden Zahlen zu eleminieren ist eigentlich einfach.
In etwa so:

Mit CreateSieve(250000000); wird zum Beispiel ein Sieb für 1..250M erstellt.
Mit AIsPrime(N) kann man dann für Zahlen im Bereich des Siebs feststellen, ob N prim ist.

Delphi-Quellcode:

var
   PSieve:Pointer;
   LastBit:NativeUInt;

FUNCTION AClearMultiples(P:NativeUInt):NativeUInt;
asm         // In : EAX==P
            mov      ecx,LastBit;
            mov      edx,PSieve;
            push     ebx               // EBX retten
            mov      ebx,eax           // V:=P
            imul     ebx,eax           // V:=P^2
            shr      ebx,1             // Bit Nr. für P^2
@ClearLoop: // Alle Vielfachen von P = 0 setzen
            btr      [edx],ebx         // PData[V]:=0;
            add      ebx,eax           // V:=V+2P
            cmp      ebx,ecx           // V<=Max
            jbe      @ClearLoop        // ja, weiter
            pop      ebx               // EBX restaurieren
            // Nächstes 1 Bit suchen
            shr      eax,1             // Bit Nr für P
@NextLoop:  add      eax,1             // P:=P+2
            bt       [edx],eax         // PData[P]=1?
            jnc      @NextLoop         // nein weiter suchen
            lea      eax,[eax*2+1]     // Nächstes P
end;

PROCEDURE CreateSieve(Numbers:NativeUInt);
var Bits,Size,Last,P:NativeUInt;
begin
   Bits:=(Numbers div 2)+(Numbers and 1); // = Anzahl ungerader Zahlen 1..Numbers
   Size:=Bits div 8; // Anzahl Bytes für Bits
   if (Bits and 7)<>0 then Inc(Size); // Eventuell 1 Byte mehr
   LastBit:=Size shl 3-1;
   GetMem(PSieve,Size);
   FillChar(PSieve^,Size,$FF);  // alle Bits=1
   PByte(PSieve)^:=$FE;         // Bit 0 (entspricht der 1) = 0
   Last:=Trunc(Sqrt(Numbers));  // Letztes P deren Vielfache zu löschen sind
   P:=3;
   while P<=Last do P:=AClearMultiples(P);
   ShowMessage('Sieb erstellt');
end;

FUNCTION AIsPrime(N:NativeUInt):boolean;
asm         // In : EAX=N
            cmp      eax,2             // 0,1: CF=1, 2:CF=0
            jbe      @SetRes           // N<=2
            shr      eax,1             // Bitnr für N (N=even: CF=0)
            jnc      @ToggleCF         // N>2, even
            mov      edx,PSieve
            bt      [edx],eax
@ToggleCF:  cmc                        // CF:=not CF
@SetRes:    setnc    al                // N=2:true, N<2:false, N>2,even:false
end;