Sieb des Eratosthenes

[Amateurprofi]

Dieser Beitrag http://www.delphipraxis.net/1213282-post18.html hat mich neugierig gemacht,
und dieser http://www.delphipraxis.net/1215236-post28.html hat mich dann dazu gebracht, mich auch mal ein wenig mit diesem Thema zu beschäftigen.

Tja, "ein wenig" ist untertrieben, ich kam vom hundertsten zum tausendsten, aber ich denke das Ergebnis ist ganz brauchbar,
auch wenn ich mit dem Zeitverhalten bei der Sieberstellung noch nicht so ganz zufrieden bin - etwas bummelig für meinen Geschmack.

Im Anhang sind eine 32 Bit Version und eine 64 Bit Version meines Programms sowie drei Help-Dateien, Help.docx, Help.pdf, Help.mht,
die über die Hilfe Funktion des Programms aufrufbar sind, solange sie sich im selben Verzeichnis befinden, wie die .exe.
Die Help Dateien sind inhaltlich identisch - aber jeder hat ja seine persönlichen Vorlieben.

Die 32 Bit Version kann Siebe bis max 2 ^ 32 - 1 erstellen, die 64 Bit Version bis 200 Giga - wenn genug Ram verfügbar ist.

Ich habe mich bemüht, alle Fehler zu finden - da ich aber weiß, wo die Juckepunkte sind, habe ich mich auf diese konzentriert und wahrscheinlich die simpelsten Ungereimtheiten übersehen.
Wer also Fehler findet, bitte Info in diesem Thread oder als PN.

Ich empfehle, die Help-Datei zumindest querzulesen.

[Horst_]

Hallo,

ich habe gerade die 32-Bit Version Primesieve32 ( Win7 ) getestet und erhalte:
Schon beim Start "External Exception C0000001D"
Bei MAX = 1e9 passiert das scheinbar mit jeder neuen Siebzahl, die beim Sieben eingeblendet wird, aber auch wenn ich das Tausender-Trennzeichen umstelle. Dubios das.
http://www.delphigroups.info/2/56/526165.html
Ohne Sourcen mußt Du es wohl selbst entdecken ;-)

Gruß Horst

[Furtbichler]

Keine Quellen? Wozu dann das Ganze?

[Amateurprofi]

Zitat:

Zitat von Horst_ (Beitrag 1217302)

Hallo,

ich habe gerade die 32-Bit Version Primesieve32 ( Win7 ) getestet und erhalte:
Schon beim Start "External Exception C0000001D"
Bei MAX = 1e9 passiert das scheinbar mit jeder neuen Siebzahl, die beim Sieben eingeblendet wird, aber auch wenn ich das Tausender-Trennzeichen umstelle. Dubios das.
http://www.delphigroups.info/2/56/526165.html
Ohne Sourcen mußt Du es wohl selbst entdecken ;-)

Gruß Horst

Hallo Horst,
Kann ich bei mir nicht nachvollziehen.
Daher ein paar Fragen:

Wenn du sagst "Schon beim Start "External Exception ..." , meinst du dann gleich wenn du das Programm startest, oder wenn du die Erstellung des Siebs startest?

Wenn du sagst "Bei MAX= 1e9 ...", meinst du dann, die Exceptions kommen bei der Erstellung des Siebs oder bei der Anzeige der Zahlen in der Liste ganz rechts?
Und: wie ist das kleineren Werten?

Läuft das bei dir als Single-Thread oder als Multi-Thread?
Anders gefragt: Welche Zahl steht in der Statusbar im dritten Panel von rechts?

Last, not least:
Welche CPU hat dein Rechner. (External Exception C0000001D heißt Illegal Instruction, es könnte also gut sein, dass ich einen Assembler Befehl nutze, den die CPU nicht "kann" - eher unwarscheinlich, aber möglich).

Wieviel Ram ist verfügbar.

Danke für Deine Mithilfe.

[Horst_]

Hallo,

AMD Phenom X4 955, Win7.
Schon zum Start kommt die Fehlermeldung.
Bei Ausführen dann extrem oft.Das II.te Bild ist bei 4E9 single Thread.
Heute morgen war es 4-fach thread.
Aus http://www.delphigroups.info/2/56/526165.html:

Zitat:

This is most probably memory corruption at work, some faulty code in your
program may have managed to overwrite a return address on the stack. When
trying to return to the mangled address the CPU ends up on some data bytes it
tries to interpret as instructions, and it fails with this error. In fact it
probably ends up in code, but in the middle of a multibyte instruction. Data
usually does not have the EXECUTE privilege, so you get an access violation,
not an invalid instruction.

It is quite easy to overwrite stack memory by incorrect usage of functions
that have untyped VAR parameters (FillChar, Move, the stream Read methods).
If you pass a pointer variable to them and forget to dereference it -> BOOM!.

Gruß Horst

[Furtbichler]

Gibt es einen Grund, uns den Quellcode vorzuenthalten?

Meinst Du, das sei ein Mehrwert für die Menschheit, den man zu Geld machen kann?

[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. ;)

[Amateurprofi]

Zitat:

Zitat von Bjoerk (Beitrag 1217420)

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:

Zitat von p80286 (Beitrag 1217432)

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;

[Amateurprofi]

Zitat:

Zitat von Morphie (Beitrag 1217445)

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.

[Furtbichler]

Na ja. Das Atkins-Sieb macht das in < 1 Sec.
Aber so, als Spielerei? Nett.

[Horst_]

Hallo,

lazarus kann das Programm nicht umwandeln.
Allein die Unit _Main hat 3800 Zeilen , Sieve um die 2500 .....
Das Trect plötzlich width und centerpoint hat mag ja noch angehen, aber TListbox.count läßt sich nichts zuweisen und das kommt andauernd vor.
Die gesamte Form wird selbst erzeugt und gestaltet. Holla, was für ein Aufwand.

Im Assembler Teil verschluckte sich Freepascal an so etwas:

Delphi-Quellcode:

ASM
mov      eax,eax.TSieve.fPData     ;hoffentlich ist das:
mov      eax,TSieve([EAX]).fPData  ;der passende Ersatz

Nach unmengen auskommentieren und umschriebseln bricht das Programm beim
"Error reading lbData.Style invalid value for property ab." und kommt deshalb wahrscheinlich auch nicht zur ursprünglichen Abbruchstelle.

Laut Debugger knallt es hier

Delphi-Quellcode:

asm
 pshufb   xmm1,xmm0

aus:
unit PrimeSieve_Progress;
Bereich Graphics:

Delphi-Quellcode:

PROCEDURE InitXMM(Range:Integer; Color1,Color2:TRGB);
FUNCTION NextColorXMM:TRGB;

AMD Phenom kennt es nicht:
http://en.wikipedia.org/wiki/SSSE3#CPUs_with_SSSE3

Also ging es um etwas nebensächliches....

Gruß Horst

[Amateurprofi]

Hallo Horst,
ja, ich treibe viel Aufwand um die Form korrekt zu gestalten.
Ich lege Wert darauf, dass der Anwender die FontSizes etc. nach seinen Wünschen einstellen kann, und das bringt dann mit sich, dass man die Größen und Positionen aller Controls selbst berechnen muss.
Da ich das aber bei all meinen Projekten so mache, ist das zum größten Teil Copy und Paste, der Aufwand hält sich also in Grenzen.

Ich werde in den nächsten Tagen alles was nicht "Standardregister" anfasst umschreiben, dann sollte das auch auf AMD Rechnern laufen. Ich hab das eh' nur gemacht um mich ein wenig in diese Technologie einzuarbeiten.
Der Performance-Gewinn an diesen Stellen ist zwar erheblich, jedoch nicht wirklich nötig.
Mit Kanonen auf Spatzen schießen, nennt man das wohl.

Das mit dem lbData.Style verstehe ich im Moment nicht, werde ich mir aber auch noch mal anschauen.

Nochmal Danke für Deine Mühe.

[Horst_]

Hallo,

ich glaube Lazarus hat wegen der Plattformunabhängigkeit nicht alles implementiert oder ist nicht immer auf dem neuesten Delphi-Stand. Siehe Trect.width etc pp.
Ich werde jetzt nicht in die Tiefen eintauchen ;-)
Oberflächen haben mich noch nie interessiert, deshalb meist Konsolenprogramme, die geschickte Darstellung mit den Hervorhebungen von Primzahlen / Zwillingen/ Vierlingen in diversen Farben ist dennoch sehr ansprechend und sehr schnell.
Jetzt weiß ich gar nicht, wie der Effekt im Original aussieht, den Du dort komponiert hast.

Gruß Horst
P.S.
@Furtbichler:
<1 Sekunde bezogen auf 1 shl 32-1 oder 1e9
alzaimars Version primgendelphi_117 gibt das auf meinen Rechner aus:
50864821 primes up to 1000359390 in 670 tics
Mein SiebmitStartzahl braucht dafür 838 Ticks aka ms