Ormiston pairs: Difference between revisions

=={{header|Pascal}}==

==={{header|Free Pascal}}===

update the digits by adding difference.

<syntaxhighlight lang=jq>

program Ormiston;

{$IFDEF FPC}{$MODE DELPHI} {$OPTIMIZATION ON,ALL}{$ENDIF}

{$IFDEF WINDOWS}{$APPLICATION CONSOLE}{$ENDIF}

uses

sysutils,strUtils;

const

Limit= 100*1000*1000;

type

tPrimesSieve = array of boolean;

tpPrimes = pBoolean;

var

PrimeSieve : tPrimesSieve;

//*********** Sieve of erathostenes

procedure ClearAll;

begin

setlength(PrimeSieve,0);

end;

function BuildWheel(pPrimes:tpPrimes;lmt:Uint64): Uint64;

var

wheelSize, wpno, pr, pw, i, k: NativeUint;

wheelprimes: array[0..15] of byte;

begin

pr := 1;//the mother of all numbers 1 ;-)

pPrimes[1] := True;

WheelSize := 1;

wpno := 0;

repeat

Inc(pr);

//pw = pr projected in wheel of wheelsize

pw := pr;

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

WheelSize := k;//the new wheelsize

//sieve multiples of the new found prime

i := pr;

i := i * i;

while i <= k do

begin

pPrimes[i] := False;

Inc(i, pr);

end;

end;

until WheelSize >= lmt;

//re-insert wheel-primes 1 still stays prime

while wpno > 0 do

begin

Dec(wpno);

pPrimes[wheelPrimes[wpno]] := True;

end;

result := pr;

end;

procedure Sieve(pPrimes:tpPrimes;lmt:Uint64);

var

sieveprime, fakt, i: UInt64;

begin

sieveprime := BuildWheel(pPrimes,lmt);

repeat

repeat

Inc(sieveprime);

until pPrimes[sieveprime];

fakt := Lmt div sieveprime;

while Not(pPrimes[fakt]) do

dec(fakt);

if fakt < sieveprime then

BREAK;

i := (fakt + 1) mod 6;

if i = 0 then

i := 4;

repeat

pPrimes[sieveprime * fakt] := False;

repeat

Dec(fakt, i);

i := 6 - i;

until (fakt < sieveprime) OR pPrimes[fakt];

if fakt < sieveprime then

BREAK;

until False;

until False;

pPrimes[1] := False;//remove 1

end;

procedure InitAndGetPrimes;

begin

setlength(PrimeSieve,Limit+2);

Sieve(@PrimeSieve[0],Limit);

PrimeSieve[Limit+1] := true;

end;

//*********** End Sieve of erathostenes

type

tDigits10 = record

d10_UsedDgts : array[0..11] of byte;

d10_dgts :array[0..15] of byte;

d10_maxIdx : Uint32;

end;

procedure Out_Digits10(const p0:tDigits10);

var

i : Int32;

begin

with p0 do

begin

i := d10_maxIdx;

while i >= 0 do

begin

write(d10_dgts[i]);

dec(i);

end;

write(' ');

For i := 0 to 9 do

write(d10_UsedDgts[i]);

writeln;

end;

end;

procedure OutIn(cnt,p1,p2:NativeInt);

Begin

write('[',Numb2USA(IntToStr(p1)):6,'|',Numb2USA(IntToStr(p2)):6,']');

if cnt MOD 5 = 0 then

writeln;

end;

function OutByPot10(cnt,prLimit:NativeInt):NativeInt;

Begin

writeln(Numb2USA(IntToStr(cnt)):12,' Ormiston pairs before ',Numb2USA(IntToStr(prLimit)):14);

result := 10*prLimit;

end;

function Convert2Digits10(p:Uint32):tDigits10;

var

i,r,m10: Uint32;

begin

fillchar(result,SizeOf(result),#0);

with result do

begin

i := 0;

repeat

r := p DIV 10;

m10 := p-10*r;

d10_dgts[i] := m10;

inc(d10_UsedDgts[m10]);

inc(i);

p := r;

until r = 0;

end;

end;

function Add2Digits10(delta:Uint32;const p0:tDigits10):tDigits10;

var

i,r,m10,dgt: Uint32;

begin

result := p0;

with result do

begin

i := 0;

repeat

m10 := d10_dgts[i];

//remove old digit

dec(d10_UsedDgts[m10]);

r := delta DIV 10;

dgt := m10 + delta -10*r;

if dgt >= 10 then

begin

dgt -= 10;

r += 1;

end;

d10_dgts[i] := dgt;

//insert new digit

inc(d10_UsedDgts[dgt]);

delta := r;

inc(i);

until (i> d10_maxIdx) OR (delta = 0);

if (i> d10_maxIdx) AND (delta>0) then

begin

d10_maxIdx += 1;

r := delta DIV 10;

dgt := delta -10*r;

d10_dgts[i] := dgt;

inc(d10_UsedDgts[dgt]);

end;

end;

end;

var

pSieve : tpPrimes;

T1,T0: TDateTime;

p1,p2 :tDigits10;

pr,pr_before,

i,cnt,prLimit :nativeInt;

Begin

T0 := now;

InitAndGetPrimes;

T1 := now;

Writeln('time for sieving ',FormatDateTime('NN:SS.ZZZ',T1-T0));

pSieve := @PrimeSieve[0];

pr_before := 2;

p1 := Convert2Digits10(pr_before);

pr := pr_before;

prLimit := 100*1000;

cnt := 0;

repeat

repeat

inc(pr);

until pSieve[pr];// pr[limit+1] set to true

if pr > limit then

BREAK;

p2 := Add2Digits10(pr-pr_before,p1);

i := 9;

repeat

if p2.d10_UsedDgts[i] <> p1.d10_UsedDgts[i] then

BREAK;

dec(i);

until i < 0;

if i < 0 then

begin

inc(cnt);

IF cnt <= 30 then

OutIn(cnt,pr_before,pr);

// Out_Digits10(p1); write(pr:10,' ');Out_Digits10(p2);

end;

p1 := p2;

pr_before := pr;

if pr >=prLimit then

prlimit:= OutByPot10(cnt,prlimit);

until false;

OutByPot10(cnt,prlimit);

end.</syntaxhighlight>

{{out|@TIO.RUN}}

<pre>

time for sieving 00:00.333

[ 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]

[40,013|40,031][40,213|40,231][40,639|40,693][45,613|45,631][48,091|48,109]

[49,279|49,297][51,613|51,631][55,313|55,331][56,179|56,197][56,713|56,731]

[58,613|58,631][63,079|63,097][63,179|63,197][64,091|64,109][65,479|65,497]

[66,413|66,431][74,779|74,797][75,913|75,931][76,213|76,231][76,579|76,597]

40 Ormiston pairs before 100,000

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.986 s User time: 0.866 s Sys. time: 0.106 s CPU share: 98.64 %

@home 1E10

time for sieving 00:17.926

...

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 0m30,479s user 0m28,191s sys 0m2,279s

</pre>