Rosetta Code

Earliest difference between prime gaps: Difference between revisions

Page Discussion

Earliest difference between prime gaps (view source)

Revision as of 13:38, 23 November 2021

9,195 bytes added ,  2 years ago
→‎{{header|Julia}}: append Pascal
m (Minor edit to C++ code)
(→‎{{header|Julia}}: append Pascal)
Line 303:
Earliest difference > 1,000,000,000 between adjacent prime gap starting primes:
Gap 276 starts at 649,580,171, gap 278 starts at 4,260,928,601, difference is 3,611,348,430.
</pre>
=={{header|Pascal}}==
==={{header|Pascal}}===
<lang pascal>program primesieve;
// sieving small ranges of 65536
//{$O+,R+}
{$IFDEF FPC}
{$MODE DELPHI}{$OPTIMIZATION ON,ALL}{$CODEALIGN proc=32}
uses
sysutils;
{$ENDIF}
{$IFDEF WINDOWS}
{$APPTYPE CONSOLE}
{$ENDIF}
 
const
smlPrimes :array [0..10] of Byte = (2,3,5,7,11,13,17,19,23,29,31);
maxPreSievePrime = 17;
sieveSize = 1 shl 15;//32768*2 ->max count of FoundPrimes = 6542
type
tSievePrim = record
svdeltaPrime:word;//diff between actual and new prime
svSivOfs:word;//-> sieveSize< 1 shl 16
svSivNum:LongWord;// 1 shl (16+32) = 2.8e14
end;
var
{$Align 16}
//primes up to 1E6-> sieving to 1E12
sievePrimes : array[0..78497] of tSievePrim;
{$Align 16}
preSieve :array[0..3*5*7*11*13*17-1] of Byte;
{$Align 16}
Sieve :array[0..sieveSize-1] of Byte;
{$Align 16}
FoundPrimes : array[0..6542] of LongWord;
{$Align 16}
Gaps : array[0..255] Of Uint64;
{$Align 16}
Limit,OffSet : Uint64;
 
SieveMaxIdx,
preSieveOffset,
SieveNum,
FoundPrimesCnt,
PrimPos,
LastInsertedSievePrime :NativeUInt;
 
procedure CopyPreSieveInSieve;forward;
procedure CollectPrimes;forward;
procedure sieveOneSieve;forward;
 
procedure preSieveInit;
var
i,pr,j,umf : NativeInt;
Begin
fillchar(preSieve[0],SizeOf(preSieve),#1);
i := 1;// starts with pr = 3
umf := 1;
repeat
IF preSieve[i] <> 0 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 umf>High(preSieve);
preSieveOffset := 0;
end;
 
procedure CalcSievePrimOfs(lmt:NativeUint);
var
i,pr : NativeUInt;
sq : Uint64;
begin
pr := 0;
i := 0;
repeat
with sievePrimes[i] do
Begin
pr := pr+svdeltaPrime;
IF sqr(pr) < (SieveSize*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
Begin
pr := pr+svdeltaPrime;
sq := sqr(pr);
svSivNum := sq DIV (2*SieveSize);
svSivOfs := ( (sq - Uint64(svSivNum)*(2*SieveSize))-1)DIV 2;
end;
end;
end;
 
procedure InitSieve;
begin
preSieveOffset := 0;
SieveNum :=0;
CalcSievePrimOfs(PrimPos-1);
end;
 
procedure InsertSievePrimes;
var
j :NativeUINt;
i,pr : NativeUInt;
begin
i := 0;
//ignore first primes already sieved with
if SieveNum = 0 then
repeat
inc(i);
until FoundPrimes[i] > maxPreSievePrime;
 
pr :=0;
j := Uint64(SieveNum)*SieveSize*2-LastInsertedSievePrime;
with sievePrimes[PrimPos] do
Begin
pr := FoundPrimes[i];
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];
svdeltaPrime := (pr-j);
j := pr;
end;
inc(PrimPos);
end;
LastInsertedSievePrime :=Uint64(SieveNum)*(SieveSize*2)+pr;
end;
 
procedure sievePrimesInit;
var
i,j,pr:NativeInt;
Begin
LastInsertedSievePrime := 0;
 
PrimPos := 0;
preSieveOffset := 0;
SieveNum :=0;
CopyPreSieveInSieve;
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
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
until false;
 
CollectPrimes;
InsertSievePrimes;
IF PrimPos < High(sievePrimes) then
Begin
InitSieve;
//Erste Sieb nochmals, aber ohne Eintrag
CopyPreSieveInSieve;
sieveOneSieve;
repeat
inc(SieveNum);
CopyPreSieveInSieve;
sieveOneSieve;
CollectPrimes;
InsertSievePrimes;
until PrimPos > High(sievePrimes);
end;
end;
 
procedure OutGaps(g1,g2,delta:NativeUint);
begin
if g2= 0 then
writeln(2*g1:4,2*g1+2:4,delta:13,Gaps[g1]:13,Gaps[g1+1]:13)
else
writeln(2*g2-1:4,2*g1:4,delta:13,Gaps[g1-1]:13,Gaps[g1]:13);
end;
function GetDiffval(val1,val2: NativeInt): NativeInt;
begin
if val1*val2 = 0 then
EXIT(0);
dec(val1,val2);
if val1<0 then
val1 :=-val1;
GetDiffval := val1;
end;
 
procedure CheckGaps;
var
val1,val2,val3,i,DekaLimit : NativeInt;
Begin
writeln('Gap1 Gap2 difference first second prime');
dekaLimit := 10;
i := 1;
repeat
val1 := Gaps[i];
if val1 <> 0 then
Begin
val2 := GetDiffval(val1,Gaps[i-1]);
val3 := GetDiffval(val1,Gaps[i+1]);
while (val2>DekaLimit) or (val3>DekaLimit) do
begin
writeln(DekaLimit:21,'<');
if val2 = 0 then
begin
if val3 > 0 then
OutGaps(i,0,val3);
end
else
begin
if val3 = 0 then
OutGaps(i,1,val2)
else
Begin
if val3 > val2 then
OutGaps(i,0,val3)
else
OutGaps(i,1,val2);
end;
end;
DekaLimit := 10*DekaLimit;
end;
end;
inc(i);
until i>=254;
end;
 
procedure CopyPreSieveInSieve;
var
lmt : NativeInt;
Begin
lmt := preSieveOffset+sieveSize;
lmt := lmt-(High(preSieve)+1);
IF lmt<= 0 then
begin
Move(preSieve[preSieveOffset],Sieve[0],sieveSize);
if lmt <> 0 then
inc(preSieveOffset,sieveSize)
else
preSieveOffset := 0;
end
else
begin
Move(preSieve[preSieveOffset],Sieve[0],sieveSize-lmt);
Move(preSieve[0],Sieve[sieveSize-lmt],lmt);
preSieveOffset := lmt
end;
end;
 
procedure CollectPrimes;
var
pSieve : pbyte;
pFound : pLongWord;
i,idx : NativeInt;
Begin
pFound := @FoundPrimes[0];
idx := 0;
i := 0;
IF SieveNum = 0 then
Begin
repeat
pFound[idx] := smlPrimes[idx];
inc(idx);
until smlPrimes[idx]>maxPreSievePrime;
i := (smlPrimes[idx] -1) DIV 2;
end;
 
pSieve := @Sieve[0];
repeat
pFound[idx]:= 2*i+1;
inc(idx,pSieve[i]);
inc(i);
until i>High(Sieve);
FoundPrimesCnt:= idx;
end;
 
procedure sieveOneSieve;
var
i,j,pr,dSievNum :NativeUint;
Begin
pr := 0;
For i := 0 to SieveMaxIdx do
with sievePrimes[i] do
begin
pr := pr+svdeltaPrime;
IF svSivNum = sieveNum then
Begin
j := svSivOfs;
repeat
Sieve[j] := 0;
inc(j,pr);
until j > High(Sieve);
 
dSievNum := j DIV SieveSize;
svSivOfs := j-dSievNum*SieveSize;
inc(svSivNum,dSievNum);
end;
end;
 
i := SieveMaxIdx+1;
repeat
if i > High(SievePrimes) then
BREAK;
with sievePrimes[i] do
begin
if svSivNum > sieveNum 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 SieveSize;
svSivOfs := j-dSievNum*SieveSize;
inc(svSivNum,dSievNum);
inc(i);
end;
until false;
end;
 
var
T1,T0,CNT,ActPrime,LastPrime,delta : Int64;
i: Int32;
 
begin
T0 := GetTickCount64;
Limit := 10*1000*1000*1000;//158*1000*1000*1000;
preSieveInit;
sievePrimesInit;
 
InitSieve;
offset := 0;
Cnt := 1;//==2
LastPrime := 2;
 
 
repeat
CopyPreSieveInSieve;sieveOneSieve;CollectPrimes;
inc(Cnt,FoundPrimesCnt);
//collect first occurrence of gap
i := 0;
repeat
ActPrime := Offset+FoundPrimes[i];
delta := (ActPrime - LastPrime) shr 1;
If Gaps[delta] = 0 then
Gaps[delta] := LastPrime;
LastPrime := ActPrime;
inc(i);
until (i >= FoundPrimesCnt);
 
inc(SieveNum);
inc(offset,2*SieveSize);
until SieveNum > (Limit DIV (2*SieveSize));
CheckGaps;
T1 := GetTickCount64;
OffSet := Uint64(SieveNum-1)*(2*SieveSize);
 
 
i := FoundPrimesCnt;
repeat
dec(i);
dec(cnt);
until (i = 0) OR (OffSet+FoundPrimes[i]<Limit);
writeln;
writeln(cnt,' in ',Limit,' takes ',T1-T0,' ms');
{$IFDEF WINDOWS}
writeln('Press <Enter>');readln;
{$ENDIF}
end.</lang>
{{out|@TIO.RUN}}
<pre>
Gap1 Gap2 difference first second prime
10<
4 6 16 7 23
100<
14 16 1718 113 1831
1000<
14 16 1718 113 1831
10000<
36 38 21042 9551 30593
100000<
70 72 141962 173359 31397
1000000<
100 102 1047576 396733 1444309
10000000<
148 150 11615524 2010733 13626257
100000000<
198 200 332037210 46006769 378043979
 
50847534 in 1000000000 takes 886 ms
//@home primes til 158,000,000,000 like C++
276 278 3611348430 649580171 4260928601
10000000000<
332 334 24933958388 5893180121 30827138509
100000000000<
386 388 121560146636 35238645587 156798792223
 
6385991032 in 158000000000 takes 247532 ms
</pre>
 
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/317758"