Rosetta Code

Erdős-Nicolas numbers: Difference between revisions

Page Discussion

Erdős-Nicolas numbers (view source)

Revision as of 14:59, 30 March 2023

3,077 bytes added ,  1 year ago
m
→‎{{header|Free Pascal}}: using fast prime decomposition less then 16 s til 1e8 on TIO.RUN.@home 7.2s
(Added C)
m (→‎{{header|Free Pascal}}: using fast prime decomposition less then 16 s til 1e8 on TIO.RUN.@home 7.2s)
Line 444:
=={{header|Pascal}}==
==={{header|Free Pascal}}===
using [[Factors_of_an_integer#using_Prime_decomposition]] , using a sort of sieve.
<syntaxhighlight lang="pascal">
program ErdoesNumb;
 
//formerly program FacOfInt;
// gets factors of consecutive integers fast
// limited to 1.2e11
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ENDIFELSE}
{$IFDEF WINDOWS}
{$APPTYPE CONSOLE}
{$ENDIF}
Line 465:
HCN_DivCnt = 4096;
type
tItem = Uint64;
tDivisors = array [0..HCN_DivCnt] of tItem;
tpDivisor = pUint64;
const
 
SizePrDeFe = 32768;//*56 <= 64kb level I or 2 Mb ~ level 2 cache or more
type
tdigits = array [0..31] of Uint32;
//the first number with 11 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31 = 2E11
Line 478 ⟶ 480:
pfDivCnt : Uint32;
pfMaxIdx : Uint32;
pfpotPrimIdx : array[0..9] of word;
pfpotMax : array[0..11] of byte;
pfpotPrimIdx : array[0..9] of word;
end;
tpPrimeFac = ^tprimeFac;
 
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;
 
var
{$ALIGN 8}
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
pdfIDX,pdfOfs: NativeInt;
 
procedure InitSmallPrimes;
Line 530 ⟶ 538:
end;
 
function smplPrimeDecompCnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):tprimeFacNativeInt;
//n must be multiple of base aka n mod base must be 0
var
q,r: Uint64;
pr,i,pot,fac,q :NativeUInt;
i : NativeInt;
Begin
fillchar(dgt,SizeOf(dgt),#0);
with result do
Begini := 0;
pfDivCntn := 1n div base;
pfSumOfDivsresult := 10;
repeat
pfRemain := n;
pfMaxIdxr := 0n;
pfpotPrimIdx[0]q := 1n div base;
pfpotMax[0]r : -= 0q*base;
n := q;
dgt[i] := r;
inc(i);
until (q = 0);
//searching lowest pot in base
result := 0;
while (result<i) AND (dgt[result] = 0) do
inc(result);
inc(result);
end;
 
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
i := 0;
var
while i < High(SmallPrimes) do
q :NativeInt;
Begin
result := 0;
q := dgt[result]+1;
if q = base then
repeat
dgt[result] := 0;
inc(result);
q := dgt[result]+1;
until q <> base;
dgt[result] := q;
result +=1;
end;
 
function SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : Uint64;
begin
n := pdfOfs;
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then
EXIT(FALSE);
//init
for i := 0 to SizePrDeFe-1 do
begin
with pdf[i] do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
end;
end;
//first factor 2. Make n+i even
i := (pdfIdx+n) AND 1;
IF (n = 0) AND (pdfIdx<2) then
i := 2;
 
repeat
with pdf[i] do
begin
prj := SmallPrimes[BsfQWord(n+i]);
qpfMaxIdx := n DIV pr1;
//ifpfpotPrimIdx[0] n:= < pr*pr0;
ifpfpotMax[0] pr:= > q thenj;
pfRemain := (n+i) BREAKshr j;
if npfSumOfDivs := pr*q(Uint64(1) shl then(j+1))-1;
pfDivCnt := j+1;
end;
i += 2;
until i >=SizePrDeFe;
//i now index in SmallPrimes
i := 0;
maxP := trunc(sqrt(n+SizePrDeFe))+1;
repeat
//search next prime that is in bounds of sieve
if n = 0 then
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if k < SizePrDeFe then
break;
until pr > MaxP;
end
else
begin
repeat
inc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if (k = pr) AND (n>0) then
k:= 0;
if k < SizePrDeFe then
break;
until pr > MaxP;
end;
 
//no need to use higher primes
if pr*pr > n+SizePrDeFe then
BREAK;
 
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
with pdf[k] do
Begin
pfpotPrimIdx[pfMaxIdx] := i;
potpfpotMax[pfMaxIdx] := 0j;
pfDivCnt *= j+1;
fac := pr;
repeat
npfRemain := qpfRemain DIV pr;
q := n div prdec(j);
pot+=1;
fac *= pr;
until n j<>= pr*q0;
pfpotMax[pfMaxIdx] := pot;
pfDivCnt *= pot+1;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
k += pr;
j := IncByBaseInBase(dgt,pr);
end;
until k >= SizePrDeFe;
until false;
 
//correct sum of & count of divisors
for i := 0 to High(pdf) do
Begin
with pdf[i] do
begin
j := pfRemain;
if j <> 1 then
begin
pfSumOFDivs *= (j+1);
pfDivCnt *=2;
end;
inc(i);
end;
pfRemain := n;
if n > 1 then
Begin
pfDivCnt *= 2;
pfSumOfDivs *= n+1
end;
end;
result := true;
end;
 
function NextSieve:boolean;
begin
dec(pdfIDX,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
result := SieveOneSieve(PrimeDecompField);
end;
 
function GetNextPrimeDecomp:tpPrimeFac;
begin
if pdfIDX >= SizePrDeFe then
if Not(NextSieve) then
EXIT(NIL);
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
end;
 
function Init_Sieve(n:NativeUint):boolean;
//Init Sieve pdfIdx,pdfOfs are Global
begin
pdfIdx := n MOD SizePrDeFe;
pdfOfs := n-pdfIdx;
result := SieveOneSieve(PrimeDecompField);
end;
 
Line 644 ⟶ 778:
 
var
pPrimeDecomp :tpPrimeFac;
Mypd : tPrimeFac;
Divs:tDivisors;
pDivs : tpDivisor;
T0:Int64;
i: NativeInt;
n,s : NativeUInt;
i : Int32;
Begin
T0 := GetTickCount64;
InitSmallPrimes;
Init_Sieve(0);
T0 := GetTickCount64;
//jump over 0
pPrimeDecomp:= GetNextPrimeDecomp;
n := 1;
pDivs := @Divs[0];
repeat
MypdpPrimeDecomp:= smplPrimeDecomp(n)GetNextPrimeDecomp;
ifs Mypd:= pPrimeDecomp^.pfSumOfDivs>2*n then;
if (s > 2*n) then
begin
GetDivisors(@Mypd,Divs)s -= n;
s// :=75% 1;of runtime
GetDivisors(pPrimeDecomp,Divs);
For i := 1 to Mypd.pfDivCnt-3 do
//calculate downwards. Not really an impact
For i := pPrimeDecomp^.pfDivCnt-2 downto 0 do
Begin
inc(s,pDivs -= Divs[i]);
if s<n then
break;
if s = n then
writeln(Format('%8d equals the sum of its first %4d divisors',
[n,i+1]));
if s>n then
break;
end;
end;
n += 1;
until n > 2100*1000*1000+1;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
end.
end.</syntaxhighlight>
</syntaxhighlight>
{{out|@TIO.RUN}}
<pre>
TIO.RUN
24 equals the sum of its first 6 divisors
2016 equals the sum of its first 31 divisors
Line 686 ⟶ 824:
392448 equals the sum of its first 68 divisors
714240 equals the sum of its first 113 divisors
1571328 equals the sum of its first 115 divisors
runtime 1.219 s
@home
...
1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors
runtime 4015.504760 s</pre>
 
Real time: 15.909 s User time: 15.721 s CPU share: 99.09 %
</pre>
 
=={{header|Perl}}==
{{libheader|ntheory}}
Horst
132

edits

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