Sehr schneller Primzahl-Finder

[GTA-Place]

Das DF (besonders Phantom1) hat mal weiter probiert und Phantom1 hat diese Funktion geschrieben:

Delphi-Quellcode:

procedure SavePrimes(MaxPrime: Cardinal; const FileName: string = '');  
const  
  CACHE = 64*1024;  
  STEMPEL: array[0..7] of Byte = (1, 7, 11, 13, 17, 19, 23, 29);  
  MODS: array[0..29] of Byte = (0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4, 0, 8, 0, 0,  
                                0, 16, 0, 32, 0, 0, 0, 64, 0, 0, 0, 0, 0, 128);  
var  
  Primes, PrimesLUT: array of Byte;  
  i, j, k, PrimeLen, PrimeBits, Num, Num2, m, mbit, s: Cardinal;  
  f: TextFile;  
begin  
  if FileName<>'' then begin  
    AssignFile(f, FileName);  
    ReWrite(f);  
    WriteLn(f, '2'+#10#13+'3'+#10#13+'5');  
  end;  
 
  SetLength(PrimesLUT, Trunc(Sqrt(MaxPrime)/30)); // max 2184 Byte für 2^32 ;-)  
  PrimesLUT[0]:=$01;  
  PrimeLen:=Length(PrimesLUT);  
  PrimeBits:=PrimeLen*30;  
  for i:=0 to Trunc(Sqrt(PrimeBits)/30) do  
    for j:=0 to 7 do  
      if PrimesLUT[i] and (1 shl j)=0 then begin  
        s:=STEMPEL[j];  
        Num:=i*30+s;  
        Num2:=Num*Num;  
        mbit:=Num2 mod 30;  
        m:=(Num2-mbit) div 30;  
        while m<PrimeLen do begin  
          PrimesLUT[m]:=PrimesLUT[m] or MODS[mbit];  
          Inc(m, i);  
          Inc(mbit, s);  
          if mbit>29 then begin  
            Dec(mbit, 30);  
            Inc(m);  
          end;  
        end;  
      end;  

  SetLength(Primes, CACHE);  
  PrimeLen:=Length(Primes);  
  PrimeBits:=PrimeLen*30;  
  for k:=0 to MaxPrime div PrimeBits do begin  
    FillChar(Primes[0], PrimeLen, 0);  
    for i:=0 to Trunc(Sqrt((k+1)*PrimeBits)/30) do  
      for j:=0 to 7 do  
        if PrimesLUT[i] and (1 shl j)=0 then begin  
          s:=STEMPEL[j];  
          Num:=i*30+s;  
          if k=0 then  
            Num2:=Num*Num  
          else  
            Num2:=Trunc(k*PrimeBits/Num)*Num+Num;  
          mbit:=Num2 mod 30;  
          m:=(Num2-mbit) div 30-k*PrimeLen;  
          while m<PrimeLen do begin  
            primes[m]:=Primes[m] or MODS[mbit];  
            Inc(m, i);  
            Inc(mbit, s);  
            if mbit>29 then begin  
              Dec(mbit, 30);  
              Inc(m);  
            end;  
          end;  
        end;  
    if FileName<>'' then  
      for i:=0 To PrimeLen-1 do  
        for j:=0 to 7 do begin  
          if k*PrimeBits+i*30+STEMPEL[j]>MaxPrime then  
            Break;  
          if not ((i=0) and (j=0) and (k=0)) and (Primes[i] and (1 shl j)=0) then  
            WriteLn(f, IntToStr(k*PrimeBits+i*30+STEMPEL[j]));  
        end;  
  end;  
  if FileName<>'' then  
    CloseFile(f);  
end;

Um 500.000.000 Zahlen zu überprüfen, braucht diese Funktion nur ca. 4 Sekunden.


EDIT: Code wurde von Phantom1 verbessert: 500 mio Zahlen in ca. 2.65 Sekunden.