Rosetta Code

Factors of an integer: Difference between revisions

Page Discussion

Factors of an integer (view source)

Revision as of 06:19, 17 July 2021

386 bytes removed ,  2 years ago
m
→‎using Prime decomposition: using TIO.RUN .delete text accidently inserted by middle mouse click
(→‎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:
//HCN(86) > 1.2E11 = 128,501,493,120 count of divs = 4096 7 3 1 1 1 1 1 1 1
HCN_DivCnt = 4096;
RECCOUNTMAX = 100*1000*1000;
DELTAMAX = 1000*1000;
type
tItem = Uint64;
tDivisors = array [0..HCN_DivCnt-1] of tItem;
tpDivisor = pUint64;
const
//used odd size for test only
SizePrDeFe = 1279132768;//*72 <= 64kb level I or 2 Mb ~ level 2 cache -32kB for DIVS
type
tdigits = array [0..31] of Uint32;
//the first number with 11 different divisorsprime factors =
// 2*3*5*7*11*13*17*19*23*29*31 = 2E11
//56 byte
tprimeFac = packed record
pfSumOfDivs,
pfRemain : Uint64;
pfpotPrimpfDivCnt : array[0..9] of Uint32;//+10*4 = 56 Byte
pfpotMaxpfMaxIdx : array[0..9] of byteUint32; //10 = 66
pfMaxIdx pfpotPrimIdx : Uint16;array[0..9] //68of word;
pfDivCntpfpotMax : Uint32;array[0..11] //72of byte;
end;
tpPrimeFac = ^tprimeFac;
Line 4,307:
 
var
{$ALIGN 8}
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
pdfIDX,pdfOfs: NativeInt;
 
procedure InitSmallPrimes;
//get primes. Number #0..65535.Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
Line 4,352 ⟶ 4,354:
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
 
end;
 
 
function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring;
Line 4,372:
if n>0 then
result += '*';
p := pFpotPrimSmallPrimes[pfpotPrimIdx[n]];
chk *= p;
str(p,s);
Line 4,401:
end;
end;
 
function smplPrimeDecomp(n:Uint64):tprimeFac;
var
Line 4,411 ⟶ 4,412:
pfRemain := n;
pfMaxIdx := 0;
pfpotPrimpfpotPrimIdx[0] := 1;
pfpotMax[0] := 0;
 
ifi Not(ODD(n)):= then0;
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
begin
pr := SmallPrimes[i];
q := n DIV pr;
//if n < pr*pr
if pr > q then
BREAK;
if n = pr*q then
Begin
pfpotPrimpfpotPrimIdx[pfMaxIdx] := pri;
pot := 0;
fac := pr;
Line 4,461 ⟶ 4,451:
 
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
//n must be multiple of base aka n mod base must be 0
var
q,r: Uint64;
Line 4,468 ⟶ 4,458:
fillchar(dgt,SizeOf(dgt),#0);
i := 0;
// dgtNum:= n;
n := n div base;
result := 0;
Line 4,479 ⟶ 4,468:
inc(i);
until (q = 0);
//searching lowest pot in base
 
result := 0;
while (result<i) AND (dgt[result] = 0) do
Line 4,502 ⟶ 4,491:
end;
 
procedurefunction SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
dgt:tDigits;
i, j, k,pr,fac,n,MaxP : NativeUInt;
begin
n := pdfOfs;
if in+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then▼
EXIT(FALSE);
//init
for i := 0 to SizePrDeFe-1 do
Line 4,517 ⟶ 4,508:
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimpfpotPrimIdx[0] := 10;
pfpotMax[0] := 0;
end;
end;
//first factor 2. Make n+i even
 
i := (pdfIdx+n) AND 1;
IF (n = 0) AND (pdfIdx<2) then
Line 4,532 ⟶ 4,522:
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfpotPrimpfpotPrimIdx[0] := 20;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
Line 4,540 ⟶ 4,530:
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
repeatif n = 0 then
inc(i);begin
repeat
if i >= High(SmallPrimes) then▼
BREAKinc(i);
pr := SmallPrimes[i];
k := pr-n MOD pr;
if (k =if pr)k AND< (n>0)SizePrDeFe then
k:= 0 break;
ifuntil kpr <> SizePrDeFe thenMaxP;
break;end
until false;else
▲ begin
▲ if i >= High(SmallPrimes) then
BREAK;repeat
▲ pot := BsfQWordinc(ni);
▲ if i >pr := High(SmallPrimes) then[i];
▲ n k := pr-n shrMOD potpr;
if (k = pr) AND (n>0) then
▲ pfMaxIdx k:= 10;
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
pfpotPrimpfpotPrimIdx[pfMaxIdx] := pri;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
Line 4,596 ⟶ 4,597:
end;
end;
result := true;
end;
 
procedurefunction NextSieve:boolean;
begin
dec(pdfIDX,SizePrDeFe);
inc(pdfOfs,SizePrDeFe);
result := SieveOneSieve(PrimeDecompField);
end;
 
Line 4,608 ⟶ 4,610:
begin
if pdfIDX >= SizePrDeFe then
if Not(NextSieve;) then
EXIT(NIL);
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
end;
 
procedurefunction Init_Sieve(n:NativeUint):boolean;
//Init Sieve pdfIdx,pdfOfs are Global
begin
pdfIdx := n MOD SizePrDeFe;
pdfOfs := n-pdfIdx;
result := SieveOneSieve(PrimeDecompField);
end;
 
Line 4,645 ⟶ 4,648:
i,len,j,l,p,k: Int32;
Begin
i := pD^.pfDivCnt;
pDivs := @Divs[0];
pDivs[0] := 1;
Line 4,657 ⟶ 4,659:
//and append them to the list
k := pfpotMax[i];
p := pfpotPrimSmallPrimes[pfpotPrimIdx[i]];
pPot :=1;
repeat
Line 4,683 ⟶ 4,685:
//Sort. Insertsort much faster than QuickSort in this special case
InsertSort(pDivs,0,len-1);
//end marker
pDivs[len] :=0;
end;
 
procedure AllFacsOut(pD:tpPrimeFac;var Divs:tdivisors;proper:boolean=true);
var
k,j: Int32;
Begin
k := pd.pfDivCnt-10;
IFj proper:= then1;
if Proper dec(k);then
IF k >j:= 0 then2;
repeat
For j := 0 to k-1 do
IF write(Divs[j],','); = 0 then
▲ i := 1BREAK;
write(Divs[k],',');
inc(j);
inc(k);
until false;
writeln(Divs[k]);
end;
Line 4,709 ⟶ 4,718:
T0 := GetTickCount64;
n := 0;
 
Init_Sieve(0);
repeat
pPrimeDecomp:= GetNextPrimeDecomp;
// GetDivisors(pPrimeDecomp,Divs);
inc(n);
until n > 10*1000*1000+1;
Line 4,719 ⟶ 4,727:
writeln('runtime ',T0/1000:0:3,' s');
GetDivisors(pPrimeDecomp,Divs);
AllFacsOut(pPrimeDecomp,Divs,true);
AllFacsOut(pPrimeDecomp,Divs,false);
writeln('simple version');
T0 := GetTickCount64;
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);
// GetDivisors(@Mypd,Divs);
inc(n);
until n > 10*1000*1000+1;
Line 4,747 ⟶ 4,740:
writeln('runtime ',T0/1000:0:3,' s');
GetDivisors(@Mypd,Divs);
AllFacsOut(@Mypd,Divs,true);
AllFacsOut(@Mypd,Divs,false);
end.</lang>
</lang>▼
{{out}}
<pre>
TIO.RUN
runtime 1.216 s▼
//out-commented GetDivisors, but still calculates sum of divisors and count of divisors▼
runtime 0.311555 s▼
1,11,909091
1,11,909091,10000001
simple version
runtime 58.854167 s
1,11,909091
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.216815 s
// calculating explicit divisisors takes the same time ~ 0.9s for 1e7
▲runtime 0.311 s
1,11,909091
1,11,909091,10000001
simple version
runtime 411.922057 s
1,11,909091
1,11,909091,10000001</pre>
Real time: 13.082 s CPU share: 99.16 %
▲</langpre>
 
=={{header|Perl}}==
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/126499"