Erdős-Nicolas numbers: Difference between revisions

Content added Content deleted
(Added C)
m (→‎{{header|Free Pascal}}: using fast prime decomposition less then 16 s til 1e8 on TIO.RUN.@home 7.2s)
Line 444: Line 444:
=={{header|Pascal}}==
=={{header|Pascal}}==
==={{header|Free Pascal}}===
==={{header|Free Pascal}}===
using [[Factors_of_an_integer#using_Prime_decomposition]] , using a sort of sieve.
<syntaxhighlight lang="pascal">
<syntaxhighlight lang="pascal">
program ErdoesNumb;
program ErdoesNumb;

//formerly program FacOfInt;
// gets factors of consecutive integers fast
// gets factors of consecutive integers fast
// limited to 1.2e11
// limited to 1.2e11
{$IFDEF FPC}
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ENDIF}
{$ELSE}
{$IFDEF WINDOWS}
{$APPTYPE CONSOLE}
{$APPTYPE CONSOLE}
{$ENDIF}
{$ENDIF}
Line 465: Line 465:
HCN_DivCnt = 4096;
HCN_DivCnt = 4096;
type
type
tItem = Uint64;
tItem = Uint64;
tDivisors = array [0..HCN_DivCnt] of tItem;
tDivisors = array [0..HCN_DivCnt] of tItem;
tpDivisor = pUint64;
tpDivisor = pUint64;
const

SizePrDeFe = 32768;//*56 <= 64kb level I or 2 Mb ~ level 2 cache or more
type
type
tdigits = array [0..31] of Uint32;
//the first number with 11 different prime factors =
//the first number with 11 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31 = 2E11
//2*3*5*7*11*13*17*19*23*29*31 = 2E11
Line 478: Line 480:
pfDivCnt : Uint32;
pfDivCnt : Uint32;
pfMaxIdx : Uint32;
pfMaxIdx : Uint32;
pfpotPrimIdx : array[0..9] of word;
pfpotMax : array[0..11] of byte;
pfpotMax : array[0..11] of byte;
pfpotPrimIdx : array[0..9] of word;
end;
end;
tpPrimeFac = ^tprimeFac;
tpPrimeFac = ^tprimeFac;

tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;
tPrimes = array[0..65535] of Uint32;


var
var
{$ALIGN 8}
SmallPrimes: tPrimes;
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
pdfIDX,pdfOfs: NativeInt;


procedure InitSmallPrimes;
procedure InitSmallPrimes;
Line 530: Line 538:
end;
end;


function smplPrimeDecomp(n:Uint64):tprimeFac;
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base aka n mod base must be 0
var
var
q,r: Uint64;
pr,i,pot,fac,q :NativeUInt;
i : NativeInt;
Begin
Begin
fillchar(dgt,SizeOf(dgt),#0);
with result do
Begin
i := 0;
pfDivCnt := 1;
n := n div base;
pfSumOfDivs := 1;
result := 0;
repeat
pfRemain := n;
pfMaxIdx := 0;
r := n;
pfpotPrimIdx[0] := 1;
q := n div base;
pfpotMax[0] := 0;
r -= q*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
begin
pr := SmallPrimes[i];
j := BsfQWord(n+i);
q := n DIV pr;
pfMaxIdx := 1;
//if n < pr*pr
pfpotPrimIdx[0] := 0;
if pr > q then
pfpotMax[0] := j;
BREAK;
pfRemain := (n+i) shr j;
if n = pr*q then
pfSumOfDivs := (Uint64(1) shl (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
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pfpotPrimIdx[pfMaxIdx] := i;
pot := 0;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
fac := pr;
fac := pr;
repeat
repeat
n := q;
pfRemain := pfRemain DIV pr;
q := n div pr;
dec(j);
pot+=1;
fac *= pr;
fac *= pr;
until n <> pr*q;
until j<= 0;
pfpotMax[pfMaxIdx] := pot;
pfDivCnt *= pot+1;
pfSumOfDivs *= (fac-1)DIV(pr-1);
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
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;
end;
inc(i);
end;
pfRemain := n;
if n > 1 then
Begin
pfDivCnt *= 2;
pfSumOfDivs *= n+1
end;
end;
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;
end;


Line 644: Line 778:


var
var
pPrimeDecomp :tpPrimeFac;
Mypd : tPrimeFac;
Divs:tDivisors;
Divs:tDivisors;
pDivs : tpDivisor;
T0:Int64;
T0:Int64;
i: NativeInt;
n,s : NativeUInt;
n,s : NativeUInt;
i : Int32;
Begin
Begin
T0 := GetTickCount64;
InitSmallPrimes;
InitSmallPrimes;
Init_Sieve(0);
T0 := GetTickCount64;
//jump over 0
pPrimeDecomp:= GetNextPrimeDecomp;
n := 1;
n := 1;
pDivs := @Divs[0];
repeat
repeat
Mypd:= smplPrimeDecomp(n);
pPrimeDecomp:= GetNextPrimeDecomp;
if Mypd.pfSumOfDivs>2*n then
s := pPrimeDecomp^.pfSumOfDivs;
if (s > 2*n) then
begin
begin
GetDivisors(@Mypd,Divs);
s -= n;
s := 1;
// 75% 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
Begin
inc(s,pDivs[i]);
s -= Divs[i];
if s<n then
break;
if s = n then
if s = n then
writeln(Format('%8d equals the sum of its first %4d divisors',
writeln(Format('%8d equals the sum of its first %4d divisors',
[n,i+1]));
[n,i]));
if s>n then
break;
end;
end;
end;
end;
n += 1;
n += 1;
until n > 2*1000*1000+1;
until n > 100*1000*1000+1;
T0 := GetTickCount64-T0;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
end.
end.</syntaxhighlight>
</syntaxhighlight>
{{out|@TIO.RUN}}
{{out|@TIO.RUN}}
<pre>
<pre>
TIO.RUN
24 equals the sum of its first 6 divisors
24 equals the sum of its first 6 divisors
2016 equals the sum of its first 31 divisors
2016 equals the sum of its first 31 divisors
Line 686: Line 824:
392448 equals the sum of its first 68 divisors
392448 equals the sum of its first 68 divisors
714240 equals the sum of its first 113 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
1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors
91963648 equals the sum of its first 142 divisors
runtime 40.504 s</pre>
runtime 15.760 s

Real time: 15.909 s User time: 15.721 s CPU share: 99.09 %
</pre>

=={{header|Perl}}==
=={{header|Perl}}==
{{libheader|ntheory}}
{{libheader|ntheory}}
Retrieved from "https://rosettacode.org/wiki/Erdős-Nicolas_numbers"