Sphenic numbers: Difference between revisions

Content added Content deleted
m (→‎{{header|AppleScript}}: Duh! The sort's not needed if the sphenics are initially stored in the list slots they index.)
m (→‎{{header|Free Pascal}}: get rid of sorting. Tested to one billion)
Line 391: Line 391:
{{trans|Wren}} Output Format {{trans|AppleScript}}
{{trans|Wren}} Output Format {{trans|AppleScript}}
Most of the time, ~ 75% in this case, is spent with sort.
Most of the time, ~ 75% in this case, is spent with sort.
Now changed to use sieve of erathostenes and insert sphenic numbers in that array,So no sort is needed.
A little bit lengthy...
<syntaxhighlight lang=pascal>
<syntaxhighlight lang=pascal>
program sphenic;
program sphenic;
{$IFDEF FPC}{$MODE DELPHI}{$ENDIF}
{$IFDEF WINDOWS}{$APPTYPE CONSOLE}{$ENDIF}
const
const
Limit= 1000*1000;
Limit= 1*1000*1000;

LimitPrime = Limit div (2*3);
type
primeCount = 15224;
tPrimesSieve = array of boolean;
var
primes : array of Uint32;
tElement = Uint32;
sphenics : array of Uint32;
tarrElement = array of tElement;
tpPrimes = pBoolean;


procedure GetPrimes;
//trial division
var
var
PrimeSieve : tPrimesSieve;
p,p2,i,maxIdx,maxSqrtPrimeIdx : Uint32;
primes : tarrElement;
begin
sphenics : tarrElement;
setlength(Primes,primeCount);

primes[0] := 2;
function BuildWheel(pPrimes:tpPrimes;lmt:Uint32): longint;
maxIdx := 0;
var
maxSqrtPrimeIdx := 0;
wheelSize, wpno, pr, pw, i, k: NativeUint;
p2:= sqr(primes[maxSqrtPrimeIdx]);
wheelprimes: array[0..15] of byte;
For p := 3 to LimitPrime do
begin
begin
pr := 1;//the mother of all numbers 1 ;-)
if p2 <=p then
pPrimes[1] := True;
begin
WheelSize := 1;
inc(maxSqrtPrimeIdx);

p2 := sqr(primes[maxSqrtPrimeIdx]);
end;
wpno := 0;
i := 0;
repeat
repeat
Inc(pr);
if p mod primes[i] = 0 then
//pw = pr projected in wheel of wheelsize
BREAK;
inc(i);
pw := pr;
until i >= maxSqrtPrimeIdx;
if pw > wheelsize then
Dec(pw, wheelsize);
if i >= maxSqrtPrimeIdx then
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], Limit - 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
begin
inc(maxIdx);
Dec(wpno);
primes[maxIdx] := p;
pPrimes[wheelPrimes[wpno]] := True;
end;
end;
result := pr;
end;
end;
setlength(primes,maxidx+1);
end;


procedure QuickSort(var A: array of Uint32);
procedure Sieve(pPrimes:tpPrimes;lmt:Uint32);
procedure QSort(L, R: Int32);
var
var
I, J: Int32;
sieveprime, fakt, i: UInt32;
Tmp, Pivot: Uint32;
begin
begin
sieveprime := BuildWheel(pPrimes,lmt);
if R - L < 1 then exit;
I := L; J := R;
{$push}{$q-}{$r-}Pivot := A[(L + R) shr 1];{$pop}
repeat
repeat
repeat
while A[I] < Pivot do Inc(I);
while A[J] > Pivot do Dec(J);
Inc(sieveprime);
if I <= J then begin
until pPrimes[sieveprime];
Tmp := A[I];
fakt := Lmt div sieveprime;
A[I] := A[J];
while Not(pPrimes[fakt]) do
A[J] := Tmp;
dec(fakt);
Inc(I); Dec(J);
if fakt < sieveprime then
end;
BREAK;
until I > J;
i := (fakt + 1) mod 6;
if i = 0 then
QSort(L, J);
QSort(I, R);
i := 4;
repeat
pPrimes[sieveprime * fakt] := False;
repeat
Dec(fakt, i);
i := 6 - i;
until pPrimes[fakt];
if fakt < sieveprime then
BREAK;
until False;
until False;
pPrimes[1] := False;//remove 1
end;
end;

procedure InitAndGetPrimes;
var
prCnt,i,idx,lmt : UInt32;
begin
begin
setlength(PrimeSieve,Limit+1);
QSort(0, High(A));
lmt := Limit DIV (2*3);
if Lmt < 65536 then
setlength(Primes,6542)
else
setlength(Primes,trunc(lmt/(ln(lmt)-1.1)));
Sieve(@PrimeSieve[0],lmt);
prCnt := 0;
for i := 1 to Lmt do
Begin
if PrimeSieve[i] then
begin
primes[prCnt] := i;
inc(prCnt);
end;
end;
setlength(primes,prCnt);
end;
end;


function binary_search(value: Uint32): Int32;
function binary_search(value: Uint32;const A:tarrElement): Int32;
var
var
p : Uint32;
p : Uint32;
l, m, h: UInt32;
l, m, h: tElement;
begin
begin
l := Low(primes);
l := Low(primes);
Line 469: Line 533:
begin
begin
m := (l + h) div 2;
m := (l + h) div 2;
p := primes[m];
p := A[m];
if p > value then
if p > value then
begin
begin
Line 493: Line 557:
p1,p2,p,i : Uint32;
p1,p2,p,i : Uint32;
begin
begin
fillchar(PrimeSieve[0],Limit+1,#0);
i := 0;
i := 0;
idx1 := trunc(exp(1/3*ln(Limit)));
idx1 := trunc(exp(1/3*ln(Limit)));
idx1 := binary_search(idx1)-1;
idx1 := binary_search(idx1,Primes)-1;
idx3 := idx1+2;
idx3 := idx1+2;
For i1 := 0 to idx1 do
For i1 := 0 to idx1 do
Line 501: Line 566:
p1 := primes[i1];
p1 := primes[i1];
idx2 := trunc(sqrt(Limit DIV p1));
idx2 := trunc(sqrt(Limit DIV p1));
idx2:= binary_search(idx2)+1;
idx2:= binary_search(idx2,Primes)+1;
For i2 := i1+1 to idx2 do
For i2 := i1+1 to idx2 do
begin
begin
Line 510: Line 575:
if p > Limit then
if p > Limit then
break;
break;
sphenics[i] := p;
//mark as sphenic number
PrimeSieve[p]:= true;
inc(i);
inc(i);
end;
end;
end;
end;
end;
end;
//insert
setlength(sphenics,i);
setlength(sphenics,i);
p := 0;
QuickSort(sphenics);
For i := 0 to Limit do
begin
if PrimeSieve[i] then
begin
sphenics[p] := i;
inc(p);
end;
end;
end;
end;


Line 524: Line 599:
end;
end;


function CheckTriplets(i:Uint32):boolean;
function CheckTriplets(i:Uint32):boolean;inline;
begin
begin
CheckTriplets:= (sphenics[i]+1=sphenics[i+1]) AND
CheckTriplets:= PrimeSieve[i] AND PrimeSieve[i+1] AND PrimeSieve[i+2];
(sphenics[i+1]+1=sphenics[i+2]);
end;
end;


Line 533: Line 607:
i,j,t5000 : Uint32;
i,j,t5000 : Uint32;
begin
begin
InitAndGetPrimes;
GetPrimes;
setlength(sphenics,21*Limit div 100+100);//203,834,084 for Limit 1E9
CreateSphenics;
CreateSphenics;
writeln('Sphenic numbers < 1,000:');
writeln('Sphenic numbers < 1,000:');
i := 1;
i := 0;
repeat
repeat
if sphenics[i] > 1000 then
if sphenics[i] > 1000 then
break;
break;
write(sphenics[i]:4);
write(sphenics[i]:4);
inc(i);
if i Mod 15 = 0 then
if i Mod 15 = 0 then
writeln;
writeln;
inc(i);
until i>= High(sphenics);
until i>= High(sphenics);
writeln;
writeln;
Line 551: Line 624:
j := 0;
j := 0;
repeat
repeat
if CheckTriplets(i) then
if CheckTriplets(sphenics[i]) then
Begin
Begin
OutTriplet(i);
OutTriplet(i);
Line 570: Line 643:
writeln('There are ',length(sphenics),' sphenic numbers < ',limit);
writeln('There are ',length(sphenics),' sphenic numbers < ',limit);
repeat
repeat
if CheckTriplets(i) then
if CheckTriplets(sphenics[i]) then
Begin
Begin
inc(j);
inc(j);
Line 579: Line 652:
until i+2 >high(sphenics);
until i+2 >high(sphenics);
writeln('There are ',j,' sphenic triplets numbers < ',limit);
writeln('There are ',j,' sphenic triplets numbers < ',limit);
writeln('The 200,000th sphenic number is ',sphenics[200000]);
writeln('The 200,000th sphenic number is ',sphenics[200000-1]);
write('The 5,000th sphenic triplet is ');OutTriplet(T5000);
write('The 5,000th sphenic triplet is ');OutTriplet(T5000);
end.
end.
</syntaxhighlight>
</syntaxhighlight>
{{out}}
{{out|@TIO.RUN}}
<pre>
<pre>
Sphenic numbers < 1,000:
Sphenic numbers < 1,000:
Line 607: Line 680:
The 200,000th sphenic number is 966467
The 200,000th sphenic number is 966467
The 5,000th sphenic triplet is {918005,918006,918007}
The 5,000th sphenic triplet is {918005,918006,918007}
Real time: 0.096 s User time: 0.067 s Sys. time: 0.028 s

@home (4.4 Ghz Ryzen 5600G):
There are 2086746 sphenic numbers < 10000000
There are 20710806 sphenic numbers < 100000000
There are 203834084 sphenic numbers < 1000000000

There are 55576 sphenic triplets numbers < 10000000
There are 527138 sphenic triplets numbers < 100000000
There are 4824694 sphenic triplets numbers < 1000000000
real 0m5,530s
</pre>
</pre>

=={{header|Phix}}==
=={{header|Phix}}==
{{trans|Wren}}
{{trans|Wren}}
Retrieved from "https://rosettacode.org/wiki/Sphenic_numbers"