Factors of an integer: Difference between revisions

Content added Content deleted
(→‎using Prime decomposition: updated from zumkeller numbers. Fast for consecutive integers.)
m (→‎using Prime decomposition: using TIO.RUN .delete text accidently inserted by middle mouse click)
Line 4,281: Line 4,281:
//HCN(86) > 1.2E11 = 128,501,493,120 count of divs = 4096 7 3 1 1 1 1 1 1 1
//HCN(86) > 1.2E11 = 128,501,493,120 count of divs = 4096 7 3 1 1 1 1 1 1 1
HCN_DivCnt = 4096;
HCN_DivCnt = 4096;
RECCOUNTMAX = 100*1000*1000;
DELTAMAX = 1000*1000;
type
type
tItem = Uint64;
tItem = Uint64;
tDivisors = array [0..HCN_DivCnt-1] of tItem;
tDivisors = array [0..HCN_DivCnt] of tItem;
tpDivisor = pUint64;
tpDivisor = pUint64;
const
const
//used odd size for test only
SizePrDeFe = 12791;//*72 <= 2 Mb ~ level 2 cache -32kB for DIVS
SizePrDeFe = 32768;//*72 <= 64kb level I or 2 Mb ~ level 2 cache
type
type
tdigits = array [0..31] of Uint32;
tdigits = array [0..31] of Uint32;
//the first number with 11 different divisors =
//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
//56 byte
tprimeFac = packed record
tprimeFac = packed record
pfSumOfDivs,
pfSumOfDivs,
pfRemain : Uint64;
pfRemain : Uint64;
pfpotPrim : array[0..9] of Uint32;//+10*4 = 56 Byte
pfDivCnt : Uint32;
pfpotMax : array[0..9] of byte; //10 = 66
pfMaxIdx : Uint32;
pfMaxIdx : Uint16; //68
pfpotPrimIdx : array[0..9] of word;
pfDivCnt : Uint32; //72
pfpotMax : array[0..11] of byte;
end;
end;
tpPrimeFac = ^tprimeFac;
tpPrimeFac = ^tprimeFac;
Line 4,307: Line 4,307:


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


procedure InitSmallPrimes;
procedure InitSmallPrimes;
//get primes Number 0..65535.Sieving only odd numbers
//get primes. #0..65535.Sieving only odd numbers
const
const
MAXLIMIT = (821641-1) shr 1;
MAXLIMIT = (821641-1) shr 1;
Line 4,352: Line 4,354:
flipflop := 3-flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
until (p > MAXLIMIT) OR (j>High(SmallPrimes));

end;
end;



function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;
function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;
Line 4,372: Line 4,372:
if n>0 then
if n>0 then
result += '*';
result += '*';
p := pFpotPrim[n];
p := SmallPrimes[pfpotPrimIdx[n]];
chk *= p;
chk *= p;
str(p,s);
str(p,s);
Line 4,401: Line 4,401:
end;
end;
end;
end;

function smplPrimeDecomp(n:Uint64):tprimeFac;
function smplPrimeDecomp(n:Uint64):tprimeFac;
var
var
Line 4,411: Line 4,412:
pfRemain := n;
pfRemain := n;
pfMaxIdx := 0;
pfMaxIdx := 0;
pfpotPrim[0] := 1;
pfpotPrimIdx[0] := 1;
pfpotMax[0] := 0;
pfpotMax[0] := 0;


if Not(ODD(n)) then
i := 0;
⚫
begin
⚫
pot := BsfQWord(n);
⚫
pfMaxIdx := 1;
pfpotPrim[0] := 2;
pfpotMax[0] := pot;
⚫
n := n shr pot;
pfRemain := n;
pfSumOfDivs := (1 shl (pot+1))-1;
pfDivCnt := pot+1;
⚫
end;

⚫
i := 1;
while i < High(SmallPrimes) do
while i < High(SmallPrimes) do
begin
begin
pr := SmallPrimes[i];
pr := SmallPrimes[i];
q := n DIV pr;
q := n DIV pr;
//if n < pr*pr
if pr > q then
if pr > q then
BREAK;
BREAK;
if n = pr*q then
if n = pr*q then
Begin
Begin
pfpotPrim[pfMaxIdx] := pr;
pfpotPrimIdx[pfMaxIdx] := i;
pot := 0;
pot := 0;
fac := pr;
fac := pr;
Line 4,461: Line 4,451:


function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base
//n must be multiple of base aka n mod base must be 0
var
var
q,r: Uint64;
q,r: Uint64;
Line 4,468: Line 4,458:
fillchar(dgt,SizeOf(dgt),#0);
fillchar(dgt,SizeOf(dgt),#0);
i := 0;
i := 0;
// dgtNum:= n;
n := n div base;
n := n div base;
result := 0;
result := 0;
Line 4,479: Line 4,468:
inc(i);
inc(i);
until (q = 0);
until (q = 0);
//searching lowest pot in base

result := 0;
result := 0;
while (result<i) AND (dgt[result] = 0) do
while (result<i) AND (dgt[result] = 0) do
Line 4,502: Line 4,491:
end;
end;


procedure SieveOneSieve(var pdf:tPrimeDecompField);
function SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
var
dgt:tDigits;
dgt:tDigits;
i, j, k,pr,fac,n : NativeUInt;
i,j,k,pr,fac,n,MaxP : NativeUInt;
begin
begin
n := pdfOfs;
n := pdfOfs;
⚫
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then
EXIT(FALSE);
//init
//init
for i := 0 to SizePrDeFe-1 do
for i := 0 to SizePrDeFe-1 do
Line 4,517: Line 4,508:
pfRemain := n+i;
pfRemain := n+i;
pfMaxIdx := 0;
pfMaxIdx := 0;
pfpotPrim[0] := 1;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
pfpotMax[0] := 0;
end;
end;
end;
end;
//first factor 2. Make n+i even
//first factor 2. Make n+i even

i := (pdfIdx+n) AND 1;
i := (pdfIdx+n) AND 1;
IF (n = 0) AND (pdfIdx<2) then
IF (n = 0) AND (pdfIdx<2) then
Line 4,532: Line 4,522:
j := BsfQWord(n+i);
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfMaxIdx := 1;
pfpotPrim[0] := 2;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := j;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfRemain := (n+i) shr j;
Line 4,540: Line 4,530:
i += 2;
i += 2;
until i >=SizePrDeFe;
until i >=SizePrDeFe;
// i now index in SmallPrimes
//i now index in SmallPrimes
i := 0;
i := 0;
maxP := trunc(sqrt(n+SizePrDeFe))+1;
repeat
repeat
//search next prime that is in bounds of sieve
//search next prime that is in bounds of sieve
repeat
if n = 0 then
inc(i);
begin
repeat
⚫
if i >= High(SmallPrimes) then
BREAK;
inc(i);
pr := SmallPrimes[i];
pr := SmallPrimes[i];
k := pr-n MOD pr;
k := pr-n MOD pr;
if (k = pr) AND (n>0) then
if k < SizePrDeFe then
k:= 0;
break;
if k < SizePrDeFe then
until pr > MaxP;
break;
end
until false;
else
⚫
begin
⚫
if i >= High(SmallPrimes) then
BREAK;
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
//no need to use higher primes
if pr*pr > n+SizePrDeFe then
if pr*pr > n+SizePrDeFe then
BREAK;
BREAK;


// j is power of prime
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
j := CnvtoBASE(dgt,n+k,pr);
repeat
repeat
with pdf[k] do
with pdf[k] do
Begin
Begin
pfpotPrim[pfMaxIdx] := pr;
pfpotPrimIdx[pfMaxIdx] := i;
pfpotMax[pfMaxIdx] := j;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
pfDivCnt *= j+1;
Line 4,596: Line 4,597:
end;
end;
end;
end;
result := true;
end;
end;


procedure NextSieve;
function NextSieve:boolean;
begin
begin
dec(pdfIDX,SizePrDeFe);
dec(pdfIDX,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
SieveOneSieve(PrimeDecompField);
result := SieveOneSieve(PrimeDecompField);
end;
end;


Line 4,608: Line 4,610:
begin
begin
if pdfIDX >= SizePrDeFe then
if pdfIDX >= SizePrDeFe then
NextSieve;
if Not(NextSieve) then
EXIT(NIL);
result := @PrimeDecompField[pdfIDX];
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
inc(pdfIDX);
end;
end;


procedure Init_Sieve(n:NativeUint);
function Init_Sieve(n:NativeUint):boolean;
//Init Sieve pdfIdx,pdfOfs are Global
//Init Sieve pdfIdx,pdfOfs are Global
begin
begin
pdfIdx := n MOD SizePrDeFe;
pdfIdx := n MOD SizePrDeFe;
pdfOfs := n-pdfIdx;
pdfOfs := n-pdfIdx;
SieveOneSieve(PrimeDecompField);
result := SieveOneSieve(PrimeDecompField);
end;
end;


Line 4,645: Line 4,648:
i,len,j,l,p,k: Int32;
i,len,j,l,p,k: Int32;
Begin
Begin
i := pD^.pfDivCnt;
pDivs := @Divs[0];
pDivs := @Divs[0];
pDivs[0] := 1;
pDivs[0] := 1;
Line 4,657: Line 4,659:
//and append them to the list
//and append them to the list
k := pfpotMax[i];
k := pfpotMax[i];
p := pfpotPrim[i];
p := SmallPrimes[pfpotPrimIdx[i]];
pPot :=1;
pPot :=1;
repeat
repeat
Line 4,683: Line 4,685:
//Sort. Insertsort much faster than QuickSort in this special case
//Sort. Insertsort much faster than QuickSort in this special case
InsertSort(pDivs,0,len-1);
InsertSort(pDivs,0,len-1);
//end marker
pDivs[len] :=0;
end;
end;


procedure AllFacsOut(pD:tpPrimeFac;var Divs:tdivisors;proper:boolean=true);
procedure AllFacsOut(var Divs:tdivisors;proper:boolean=true);
var
var
k,j: Int32;
k,j: Int32;
Begin
Begin
k := pd.pfDivCnt-1;
k := 0;
IF proper then
j := 1;
dec(k);
if Proper then
IF k > 0 then
j:= 2;
repeat
For j := 0 to k-1 do
write(Divs[j],',');
IF Divs[j] = 0 then
⚫
BREAK;
write(Divs[k],',');
inc(j);
inc(k);
until false;
writeln(Divs[k]);
writeln(Divs[k]);
end;
end;
Line 4,709: Line 4,718:
T0 := GetTickCount64;
T0 := GetTickCount64;
n := 0;
n := 0;

Init_Sieve(0);
Init_Sieve(0);
repeat
repeat
pPrimeDecomp:= GetNextPrimeDecomp;
pPrimeDecomp:= GetNextPrimeDecomp;
// GetDivisors(pPrimeDecomp,Divs);
GetDivisors(pPrimeDecomp,Divs);
inc(n);
inc(n);
until n > 10*1000*1000+1;
until n > 10*1000*1000+1;
Line 4,719: Line 4,727:
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
GetDivisors(pPrimeDecomp,Divs);
GetDivisors(pPrimeDecomp,Divs);
AllFacsOut(pPrimeDecomp,Divs,true);
AllFacsOut(Divs,true);
AllFacsOut(pPrimeDecomp,Divs,false);
AllFacsOut(Divs,false);
writeln('simple version');
writeln('simple version');
T0 := GetTickCount64;
T0 := GetTickCount64;
n := 0;
n := 0;
repeat
repeatSpeedTest 0.296 secs for 1..1000000
mean count of divisors 13.970
**SpeedTest 5.707 secs for 1..10000000
**mean count of divisors 16.273
****SpeedTest 121.159 secs for 1..100000000
****mean count of divisors 18.575

Enter a number between 1 and 4294967295:
3491888400 is a nice choice :
123456789
Proper number version
123456789 = 3^2*3607*3803 got 11 proper divisors with sum : 54966027
1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263
including n version
123456789 = 3^2*3607*3803 got 12 divisors with sum : 178422816
1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789
Mypd:= smplPrimeDecomp(n);
Mypd:= smplPrimeDecomp(n);
// GetDivisors(@Mypd,Divs);
GetDivisors(@Mypd,Divs);
inc(n);
inc(n);
until n > 10*1000*1000+1;
until n > 10*1000*1000+1;
Line 4,747: Line 4,740:
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
GetDivisors(@Mypd,Divs);
GetDivisors(@Mypd,Divs);
AllFacsOut(@Mypd,Divs,true);
AllFacsOut(Divs,true);
AllFacsOut(@Mypd,Divs,false);
AllFacsOut(Divs,false);
end.
end.</lang>
⚫
</lang>
{{out}}
{{out}}
<pre>
<pre>
TIO.RUN
⚫
runtime 1.216 s
⚫
//out-commented GetDivisors, but still calculates sum of divisors and count of divisors
⚫
runtime 0.555 s
1,11,909091
1,11,909091
1,11,909091,10000001
1,11,909091,10000001
simple version
simple version
runtime 5.854 s
runtime 8.167 s
1,11,909091
1,11,909091
1,11,909091,10000001
1,11,909091,10000001
Real time: 8.868 s CPU share: 99.57 %

//with GetDivisors
⚫
//out-commented GetDivisors, but still calculates sum of divisors and count of divisors
⚫
runtime 1.815 s
// calculating explicit divisisors takes the same time ~ 0.9s for 1e7
⚫
runtime 0.311 s
1,11,909091
1,11,909091
1,11,909091,10000001
1,11,909091,10000001
simple version
simple version
runtime 4.922 s
runtime 11.057 s
1,11,909091
1,11,909091
1,11,909091,10000001</pre>
1,11,909091,10000001
Real time: 13.082 s CPU share: 99.16 %
⚫
</pre>


=={{header|Perl}}==
=={{header|Perl}}==
Retrieved from "https://rosettacode.org/wiki/Factors_of_an_integer"