Rosetta Code

Ormiston pairs: Difference between revisions

Page Discussion

Ormiston pairs (view source)

Revision as of 13:42, 22 February 2023

4,812 bytes added ,  1 year ago
m
→‎{{header|Free Pascal}}: changed to segmented prime sieve. Generating and extract primes takes the most time.
m (→‎{{header|Free Pascal}}: checking every second number for prime....useful especially for bigger numbers)
m (→‎{{header|Free Pascal}}: changed to segmented prime sieve. Generating and extract primes takes the most time.)
Line 420:
uses
sysutils,strUtils;
const
Limit= 100*1000*1000;
 
//********* segmented sieve of erathostenes *********
{segmented sieve of Erathostenes using only odd numbers}
{using presieved sieve of small primes, to reduce the most time consuming}
const
smlPrimes :array [0..10] of Byte = (2,3,5,7,11,13,17,19,23,29,31);
maxPreSievePrimeNum = 8;
maxPreSievePrime = 19;//smlPrimes[maxPreSievePrimeNum];
cSieveSize = 2*16384;//<= High(Word)+1 // Level I Data Cache
type
tSievePrim = record
tPrimesSieve = array of boolean;
svdeltaPrime:word;//diff between actual and new prime
tpPrimes = pBoolean;
svSivOfs:word; //Offset in sieve
svSivNum:LongWord;//1 shl (1+16+32) = 5.6e14
end;
tpSievePrim = ^tSievePrim;
 
var
//sieved with primes 3..maxPreSievePrime.here about 255255 Byte
PrimeSieve : tPrimesSieve;
{$ALIGN 32}
//*********** Sieve of erathostenes
preSieve :array[0..3*5*7*11*13*17*19-1] of Byte;//must be > cSieveSize
procedure ClearAll;
{$ALIGN 32}
begin
Sieve :array[0..cSieveSize-1] of Byte;
setlength(PrimeSieve,0);
{$ALIGN 32}
//prime = FoundPrimesOffset + 2*FoundPrimes[0..FoundPrimesCnt]
FoundPrimes : array[0..12252] of Word;
{$ALIGN 32}
sievePrimes : array[0..1077863] of tSievePrim;
FoundPrimesOffset : Uint64;
FoundPrimesCnt,
FoundPrimesIdx,
FoundPrimesTotal,
SieveNum,
SieveMaxIdx,
preSieveOffset,
LastInsertedSievePrime :NativeUInt;
 
procedure CopyPreSieveInSieve; forward;
procedure CollectPrimes; forward;
procedure sieveOneSieve; forward;
procedure Init0Sieve; forward;
procedure SieveOneBlock; forward;
 
procedure preSieveInit;
var
i,pr,j,umf : NativeInt;
Begin
fillchar(preSieve[0],SizeOf(preSieve),#1);
i := 1;
pr := 3;// starts with pr = 3
umf := 1;
repeat
IF preSieve[i] =1 then
Begin
pr := 2*i+1;
j := i;
repeat
preSieve[j] := 0;
inc(j,pr);
until j> High(preSieve);
umf := umf*pr;
end;
inc(i);
until (pr = maxPreSievePrime)OR(umf>High(preSieve)) ;
preSieveOffset := 0;
end;
 
function InsertSievePrimes(PrimPos:NativeInt):NativeInt;
var
j :NativeUINt;
i,pr : NativeUInt;
begin
i := 0;
//ignore first primes already sieved with
if SieveNum = 0 then
i := maxPreSievePrimeNum;
pr :=0;
j := Uint64(SieveNum)*cSieveSize*2-LastInsertedSievePrime;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := pr+j;
j := pr;
end;
inc(PrimPos);
for i := i+1 to FoundPrimesCnt-1 do
Begin
IF PrimPos > High(sievePrimes) then
BREAK;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i]*2+1;
svdeltaPrime := (pr-j);
j := pr;
end;
inc(PrimPos);
end;
LastInsertedSievePrime :=Uint64(SieveNum)*cSieveSize*2+pr;
result := PrimPos;
end;
 
procedure CalcSievePrimOfs(lmt:NativeUint);
function BuildWheel(pPrimes:tpPrimes;lmt:Uint64): Uint64;
//lmt High(sievePrimes)
var
var
wheelSize, wpno, pr, pw, i, k: NativeUint;
i,pr : NativeUInt;
wheelprimes: array[0..15] of byte;
sq : Uint64;
begin
pr := 0;
i := 0;
repeat
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
IF sqr(pr) < (cSieveSize*2) then
Begin
svSivNum := 0;
svSivOfs := (pr*pr-1) DIV 2;
end
else
Begin
SieveMaxIdx := i;
pr := pr-svdeltaPrime;
BREAK;
end;
end;
inc(i);
until i > lmt;
 
for i := i to lmt do
begin
with sievePrimes[i] do
pr := 1;//the mother of all numbers 1 ;-)
Begin
pPrimes[1] := True;
WheelSize pr := 1pr+svdeltaPrime;
sq := sqr(pr);
svSivNum := sq DIV (2*cSieveSize);
svSivOfs := ( (sq - Uint64(svSivNum)*(2*cSieveSize))-1)DIV 2;
end;
end;
end;
 
procedure sievePrimesInit;
wpno := 0;
var
i,j,pr,PrimPos:NativeInt;
Begin
LastInsertedSievePrime := 0;
preSieveOffset := 0;
SieveNum :=0;
CopyPreSieveInSieve;
//normal sieving of first sieve
i := 1; // start with 3
repeat
while Sieve[i] = 0 do
inc(i);
pr := 2*i+1;
inc(i);
j := ((pr*pr)-1) DIV 2;
if j > High(Sieve) then
BREAK;
repeat
Inc(pr)Sieve[j] := 0;
inc(j,pr);
//pw = pr projected in wheel of wheelsize
until j pw> := prHigh(Sieve);
until false;
if pw > wheelsize then
Dec(pw, wheelsize);
if pPrimes[pw] then
begin
k := WheelSize + 1;
//turn the wheel (pr-1)-times
for i := 1 to pr - 1 do
begin
Inc(k, WheelSize);
if k < lmt then
move(pPrimes[1], pPrimes[k - WheelSize], WheelSize)
else
begin
move(pPrimes[1], pPrimes[k - WheelSize], Lmt - WheelSize * i);
break;
end;
end;
Dec(k);
if k > lmt then
k := lmt;
wheelPrimes[wpno] := pr;
pPrimes[pr] := False;
Inc(wpno);
 
CollectPrimes;
WheelSize := k;//the new wheelsize
PrimPos := InsertSievePrimes(0);
//sieve multiples of the new found prime
LastInsertedSievePrime := FoundPrimes[PrimPos]*2+1;
i := pr;
i := i * i;
while i <= k do
begin
pPrimes[i] := False;
Inc(i, pr);
end;
end;
until WheelSize >= lmt;
 
IF PrimPos < High(sievePrimes) then
//re-insert wheel-primes 1 still stays prime
Begin
while wpno > 0 do
beginInit0Sieve;
Dec(wpno)sieveOneBlock;
repeat
pPrimes[wheelPrimes[wpno]] := True;
end sieveOneBlock;
result := prdec(SieveNum);
PrimPos := InsertSievePrimes(PrimPos);
inc(SieveNum);
until PrimPos > High(sievePrimes);
end;
Init0Sieve;
end;
 
procedure Sieve(pPrimes:tpPrimes;lmt:Uint64)Init0Sieve;
begin
var
FoundPrimesTotal :=0;
sieveprime, fakt, i: UInt64;
preSieveOffset := 0;
SieveNum :=0;
CalcSievePrimOfs(High(sievePrimes));
end;
 
procedure CopyPreSieveInSieve;
var
lmt : NativeInt;
Begin
lmt := preSieveOffset+cSieveSize;
lmt := lmt-(High(preSieve)+1);
IF lmt<= 0 then
begin
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize);
sieveprime := BuildWheel(pPrimes,lmt);
if lmt <> 0 then
repeat
inc(preSieveOffset,cSieveSize)
repeat
else
Inc(sieveprime);
preSieveOffset := 0;
until pPrimes[sieveprime];
end
fakt := Lmt div sieveprime;
else
while Not(pPrimes[fakt]) do
begin
dec(fakt);
Move(preSieve[preSieveOffset],Sieve[0],cSieveSize-lmt);
if fakt < sieveprime then
Move(preSieve[0],Sieve[cSieveSize-lmt],lmt);
BREAK;
ipreSieveOffset := (fakt + 1) mod 6;lmt
end;
if i = 0 then
end;
i := 4;
 
repeat
procedure sieveOneSieve;
pPrimes[sieveprime * fakt] := False;
var
sp:tpSievePrim;
pSieve :pByte;
i,j,pr,sn,dSievNum :NativeUint;
Begin
pr := 0;
sn := sieveNum;
sp := @sievePrimes[0];
pSieve := @Sieve[0];
For i := SieveMaxIdx downto 0 do
with sp^ do
begin
pr := pr+svdeltaPrime;
IF svSivNum = sn then
Begin
j := svSivOfs;
repeat
Dec(fakt,pSieve[j] i):= 0;
i := 6 - iinc(j,pr);
until (faktj <> sieveprimeHigh(Sieve) OR pPrimes[fakt];
ifdSievNum fakt:= <j sieveprimeDIV thencSieveSize;
svSivOfs := BREAKj-dSievNum*cSieveSize;
until False svSivNum := sn+dSievNum;
until False end;
inc(sp);
pPrimes[1] := False;//remove 1
end;
i := SieveMaxIdx+1;
repeat
if i > High(SievePrimes) then
BREAK;
with sp^ do
begin
if svSivNum > sn then
Begin
SieveMaxIdx := I-1;
Break;
end;
pr := pr+svdeltaPrime;
j := svSivOfs;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
dSievNum := j DIV cSieveSize;
svSivOfs := j-dSievNum*cSieveSize;
svSivNum := sn+dSievNum;
end;
inc(i);
inc(sp);
until false;
end;
 
procedure CollectPrimes;
//extract primes to FoundPrimes
var
pSieve : pbyte;
pFound : pWord;
i,idx : NativeUint;
Begin
FoundPrimesOffset := SieveNum*2*cSieveSize;
FoundPrimesIdx := 0;
pFound :=@FoundPrimes[0];
i := 0;
idx := 0;
IF SieveNum = 0 then
//include small primes used to pre-sieve
Begin
repeat
pFound[idx]:= (smlPrimes[idx]-1) DIV 2;
inc(idx);
until smlPrimes[idx]>maxPreSievePrime;
i := (smlPrimes[idx] -1) DIV 2;
end;
//grabbing the primes without if then -> reduces time extremly
//primes are born to let branch-prediction fail.
pSieve:= @Sieve[Low(Sieve)];
repeat
//store every value until a prime aka 1 is found
pFound[idx]:= i;
inc(idx,pSieve[i]);
inc(i);
until i>High(Sieve);
FoundPrimesCnt:= idx;
inc(FoundPrimesTotal,Idx);
end;
 
procedure InitAndGetPrimesSieveOneBlock;
begin
CopyPreSieveInSieve;
setlength(PrimeSieve,Limit+2);
sieveOneSieve;
Sieve(@PrimeSieve[0],Limit);
CollectPrimes;
PrimeSieve[Limit+1] := true;
inc(SieveNum);
end;
//*********** End Sieve of erathostenes
 
function Nextprime:Uint64;
Begin
result := FoundPrimes[FoundPrimesIdx]*2+1+FoundPrimesOffset;
if (FoundPrimesIdx=0) AND (sievenum = 1) then
inc(result);
inc(FoundPrimesIdx);
If FoundPrimesIdx>= FoundPrimesCnt then
SieveOneBlock;
end;
 
function PosOfPrime: Uint64;inline;
Begin
result := FoundPrimesTotal-FoundPrimesCnt+FoundPrimesIdx;
end;
 
function TotalCount :Uint64;inline;
begin
result := FoundPrimesTotal;
end;
 
function SieveSize :LongInt;inline;
Begin
result := 2*cSieveSize;
end;
 
function SieveStart:Uint64;inline;
Begin
result := (SieveNum-1)*2*cSieveSize;
end;
 
procedure InitPrime;
Begin
Init0Sieve;
SieveOneBlock;
end;
//********* segmented sieve of erathostenes *********
 
const
Limit= 10*1000*1000*1000;
type
tDigits10 = array[0..15] of byte;
td10_UsedDgts2 = array[0..3] of Uint32;
td10_UsedDgts3 = array[0..1] of Uint64;
tpd10_UsedDgts2 = ^td10_UsedDgts2;
tpd10_UsedDgts3 = ^td10_UsedDgts3;
var
p_zero : tDigits10;
 
procedure OutIn(cnt,p1,p2:NativeInt);
Line 556 ⟶ 789:
r : NativeUint;
begin
// fillchar(outP,SizeOf(outP),#0);//takes longer
outP := p_zero;
td10_UsedDgts3(outP)[0]:=0;td10_UsedDgts3(outP)[1]:=0;
repeat
r := p DIV 10;
Line 564 ⟶ 798:
end;
 
function CheckOrmiston(const d1,d2:tpd10_UsedDgts2tpd10_UsedDgts3):boolean;inline;
begin
result := (d1^[1]=d2^[1]) AND (d1^[0]=d2^[0]) AND (d1^[21]=d2^[21]);
end;
 
var
pSieve : tpPrimes;
T1,T0: TDateTime;
{$align 16}
p1,p2 :tDigits10;
pr,pr1,prLimit :nativeInt;
delta,cnt : NativeUint;
Begin
preSieveInit;
T0 := now;
sievePrimesInit;
InitAndGetPrimes;
T1 := nowInitPrime;
 
Writeln('time for sieving ',FormatDateTime('NN:SS.ZZZ',T1-T0));
pSieve := @PrimeSieve[0];
pr := 3;
prLimit := 100*1000;
cnt := 0;
pr1 := nextprime;
repeat
deltapr := 0nextprime;
repeat
pr +=2;
delta += 2;
until pSieve[pr];
if pr > limit then
BREAK;
if (pr-pr1) mod 18 = 0 then
 
if delta mod 18 = 0 then
begin
Convert2Digits10(prpr1,p1);
Convert2Digits10(pr-delta,p2);
if CheckOrmiston(@p1,@p2) then
begin
inc(cnt);
IF cnt <= 30 then
OutIn(cnt,pr-deltapr1,pr);
end;
end;
if pr >=prLimit then
prlimit:= OutByPot10(cnt,prlimit);
pr1:= pr;
until false;
OutByPot10(cnt,prlimit);
end.
end.</syntaxhighlight>
</syntaxhighlight>
{{out|@TIO.RUN}}
<pre>
//only get all primes to 1E10 Real time: 13.294 s CPU share: 99.11 %
time for sieving 00:00.313
[ 1,913| 1,931][18,379|18,397][19,013|19,031][25,013|25,031][34,613|34,631]
[35,617|35,671][35,879|35,897][36,979|36,997][37,379|37,397][37,813|37,831]
Line 621 ⟶ 850:
382 Ormiston pairs before 1,000,000
3,722 Ormiston pairs before 10,000,000
34,901 Ormiston pairs before 100,000,000
Real time: 0.620 s User time: 0.528 s Sys. time: 0.085 s CPU share: 98.95 %
@home 1E10
// Limit= 10*1000*1000*1000;
time for sieving 00:17.708
....
34,901 Ormiston pairs before 100,000,000
326,926 Ormiston pairs before 1,000,000,000
3,037,903 Ormiston pairs before 10,000,000,000
Real time: 21.114 s User time: 20.862 s Sys. time: 0.057 s CPU share: 99.07 %
 
@home real 0m22 0m6,777s873s user 0m20 0m6,654s864s sys 0m2 0m0,115s008s
</pre>
 
Horst
132

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/337454"