Factors of an integer: Difference between revisions
Content added Content deleted
m (→{{header|REXX}}: changed a comment.) |
(→using Prime decomposition: updated from zumkeller numbers. Fast for consecutive integers.) |
||
Line 4,263: | Line 4,263: | ||
Insertion sort was much faster, because mostly not so many factors need to be sorted.<BR> |
Insertion sort was much faster, because mostly not so many factors need to be sorted.<BR> |
||
"runtime overhead" +25% instead +100% for quicksort against no sort.<BR> |
"runtime overhead" +25% instead +100% for quicksort against no sort.<BR> |
||
Especially fast for consecutive integers. |
|||
[https://tio.run/##rVjdc9pIEn/nr5iHrQJs8WU7XlsE3yU23nDlGG6draurlG9LiAEGhKSVRibelP/183X3fGiEYfce1g8JTHf/uqe/h2S64qFspUEeBlFrnoavr2mWLLJgw26DcDwfxZIveNavff9hdHszvGW3k@uXWqfD2Pcffj72xscvbManxULEC5YWWZrkvIa0z@ObIbsZ3k0@jV7oYJxKsRG/B1IkMRvfex/u7hThOpnxD5FYxAwUh4OLF1A1vHsYIhXIHyaTL/@eDNn1@P5hfDck6v3N6PalVuQ8B578OS@kiPJ@rSafU9QuJ4lkA5bxMMlmiGL@0kROMrHxds4@B9@Y/4uI5elJvyTxeAaQgJaCCAdn7IXUIHPEzZnPgiwLnr922@3e6SNL5mRLf5/AdSy9fec34ukQ6aHYjJGe76XeFxv84LM3N9m5z8MmiCK0l@dwJWvw@buzHpn8L7giSsuR5AQJXPcQtyeOyEQBI0SeZFaeborsRE01GaipK/gUZEB21PsVY1CWfH0DXt6kful5sBtzg8@KjLNRnPNMPiSZbKAi36rrszs@lx77WSyWEvygVEMCs6bRPfLYPxwCapyIpwS5tflTDqkMx3MwX/gD1mPHBMtkooFnCZANm5YHRrTiq3hUbl/BgWAt1lNft0sRcdZYsasBgTVZEM8Ymf919ciuFEhTQbvgTMGujnuP/kCzm8CClxqrpvpGsd1hZymC4jFRVfhdN4qq85WH0syTPIcy8VbeJviGERCzb25SiZwkwGXTJIl4EPdrH7XBpQAD/d1@Ndpfu494fKLijB9P8WPGUx5Isl7kqYIeMJkVXPvSYlV4Ic/ITqS6SkoHjW41y5HmvEK1csljWxrTjAdryz9H@iaZGWRQ67J/dMJSsXUeRDnv7we1kQGJuAxYEUPLMoLq7nMLaXW6GlG69G/TgLpXL8kU/mwnNxABwntCssqAFJ3SY/kyYr1z1uq9TZOJLsFxAeWWJSnP/I8q7AMMka0s4a0hI0YqSUwCb4VcuiWtEry81TYTkjeodXl1cHddX2uNXqUW2dIVdIvlSImAhbhu9UytGCWqAUMBlgQVUdXesZArkbfK9UioH9WblsgjmmIHOP9T9xQqCekKd41YO0b8gQkVWKO8ovqQYq1UYMJiSFxckoniRp0toK@BhJ4orR64WLPPTPsms3MYHD5x2gHTopgomxyLDmED8p9CNnc60byIQ9oFECKBdGyESZzLygSwA6Dp6zmCTVvEkYh5mWUZz4tIUsqUsm1j2xt9xqK/XJ@96oEqUoxJLtCMRmw2DqgfhHSUlzUF5ToPQi8EF2tP/1Yk8juUwEbEuE5x5XgAzF/syFPT9mNpblttB2BxjFEAONtTFSx@pUoTlgDJFbP3TtI6ZYshVrAq4bs24d2OpfPWqnW2LavPaU@hYwYc1mzmOU0/fdPuzcSFZp9mR2l2FVdatmrF9A09R7agy2x3VHJEQ5JbSPqWZS1XB0Cl5MGHlc5TZa34wmiurp92xKk/iPpb3sr0U38xcqH11XXP3JX@d69rpxfOAreXuIqPdhWbYaH89P6qjCctIiqFQEpDHvccUGcAKrXA23w7HZEinMlkZ2qnE7CpWKhaYxGuYskTz/5Wo9jFNL9s1HZj5obhD@JVjda@vHUD1Tsw292LlT45qaS5qhogAU6oFlB1VNahaVuGYtuKTpLd3vITx3U8b2Af2dvMqMMgi2/3dttiUgVc7uxU@l7EY9j/Ii/1Uni/eGsPGnplbzYeE8Zi6uCOevIEBEhAgO54vJDLBk2BMlY5l1FJ8FTwU3PRv@MHaCz2UG@P5HyQKz/jJzjAz7lqcnRe3RiqnZrGoQn/7pLSfrtLuPnR6XyG7i/S6JlxyMNnM/lga4c3A1dwcEkW862ChZKS4HMrj7t/kKYQR@TboHnIH4nc1LbafyjhWjbbUn2IxtmzCdW55dnpEETGai7TGN2ysm4Bv7mr1J6@Ra6PaJ8EsCP1vdyxKZDg9eOBZXVJWBbRoWYwgwfM2hJV2a/f2zaogxxV6ufAxEVGMEP1C3watvUzEQIj2aYIl9BQckmbUhCzfxYiXCMbGAgnImd5ykMRRCwMaOTsvDFzr@uRq8Dami7AW7NPqPqjsqoO9Z2S8/VyrGe7@Wo8vm9F8HYL6i8sN92kTfdwsDy7rpWbzgGdu73oQxRBw8nxsRDbN2PVD339ipjue0WsvZV9RBjHUDFUnKwMhMrodNjvPEvYFKI2Mz6sPFaMX3b35DL3gLgGp3XfrjnVUqm@OWgx17Xg1T37DNALMlHWj82Dr@@HlPPZF3hmNu5gjZM6bawfvkB@5h6ktxkL6JTzM72wqUXbL3@GwfO42EzBreyXHfcRlHYhfqYwKigHX68feGUFRH0U7032mZtjPVsprxIUkvJsttPV7XvBXhfeBI2qIS0ysNnpdbtdv@ufwlsi52FOP8D02u26RwZU4TaQOCxEefzhyT496qXHOhG5FfBcSRMH5eQ/86Rx5PCn0T31BLH7a1WlpkjeU1VVhhevdYT/NGkNtVcYxtiNAuPtKZdbDg2vR78MnZ1cnl2e/3hy@c7Xb2Ird3p22bu4uDjrdhm0LZAXIWfhMsH//Lqu2mAGrArZ5LjWczXYl@kWfaLKRDPDeMOXhUlut7h3oq1XbIsDeRIVM/xNOP7/UTza@cqawV6@lDLN/U6Hx@2tWIuUz0TQTrJFB791PonFMnr@1fRM/qs2HJ4tn67vLQDIb7fbZQJRaxexaE0Fhz2SR7P2jHeCcCk2nSUhteU3SbN4LjLI1d4JeBmBhvc37dfX3snp2bvzHy8u/xvOo2CRv7bGp/8D Try it online!] |
|||
<lang pascal>program |
<lang pascal>program FacOfInt; |
||
// gets factors of consecutive integers fast |
|||
// limited to 1.2e11 |
|||
{$IFDEF FPC} |
{$IFDEF FPC} |
||
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON} |
|||
// {$R+,O+} debuging purpose |
|||
{$MODE DELPHI} |
|||
{$Optimization ON,ALL} |
|||
{$CodeAlign proc=8} |
|||
{$ELSE} |
{$ELSE} |
||
{$APPTYPE CONSOLE} |
|||
{$ENDIF} |
{$ENDIF} |
||
uses |
uses |
||
sysutils |
sysutils |
||
{$IFDEF WINDOWS},Windows{$ENDIF} |
|||
; |
|||
//###################################################################### |
|||
//prime decomposition |
|||
const |
|||
//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 |
type |
||
tItem = Uint64; |
|||
tDivisors = array [0..HCN_DivCnt-1] of tItem; |
|||
potPrim, |
|||
tpDivisor = pUint64; |
|||
potMax :Uint32; |
|||
const |
|||
end; |
|||
SizePrDeFe = 12791;//*72 <= 2 Mb ~ level 2 cache -32kB for DIVS |
|||
type |
|||
tprimeFac = record |
|||
tdigits = array [0..31] of Uint32; |
|||
//the first number with 11 different divisors = |
|||
pfCnt, |
|||
// 2*3*5*7*11*13*17*19*23*29*31 = 2E11 |
|||
pfDivCnt, |
|||
tprimeFac = packed record |
|||
pfSumOfDivs, |
pfSumOfDivs, |
||
pfRemain : Uint64; |
|||
pfpotPrim : array[0..9] of Uint32;//+10*4 = 56 Byte |
|||
pfpotMax : array[0..9] of byte; //10 = 66 |
|||
pfMaxIdx : Uint16; //68 |
|||
pfDivCnt : Uint32; //72 |
|||
end; |
end; |
||
tpPrimeFac = ^tprimeFac; |
|||
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac; |
|||
tPrimes = array[0..65535] of Uint32; |
|||
tItem = NativeUint; |
|||
tDivisors = array of tItem; |
|||
tpDivisor = pNativeUint; |
|||
var |
|||
SmallPrimes: tSmallPrimes; |
|||
primeDecomp: tprimeFac; |
|||
procedure InsertSort(pDiv:tpDivisor; Left, Right : NativeInt ); |
|||
var |
var |
||
SmallPrimes: tPrimes; |
|||
I, J: NativeInt; |
|||
PrimeDecompField :tPrimeDecompField; |
|||
Pivot : tItem; |
|||
pdfIDX,pdfOfs: NativeInt; |
|||
begin |
|||
for i:= 1 + Left to Right do |
|||
begin |
|||
Pivot:= pDiv[i]; |
|||
j:= i - 1; |
|||
while (j >= Left) and (pDiv[j] > Pivot) do |
|||
begin |
|||
pDiv[j+1]:=pDiv[j]; |
|||
Dec(j); |
|||
end; |
|||
pDiv[j+1]:= pivot; |
|||
end; |
|||
end; |
|||
procedure InitSmallPrimes; |
procedure InitSmallPrimes; |
||
//get primes Number 0..65535.Sieving only odd numbers |
|||
const |
|||
MAXLIMIT = (821641-1) shr 1; |
|||
var |
var |
||
pr : array[0..MAXLIMIT] of byte; |
|||
pr,testPr,j,maxprimidx: Uint32; |
|||
p,j,d,flipflop :NativeUInt; |
|||
isPrime : boolean; |
|||
Begin |
Begin |
||
maxprimidx := 0; |
|||
SmallPrimes[0] := 2; |
SmallPrimes[0] := 2; |
||
fillchar(pr[0],SizeOf(pr),#0); |
|||
pr := 3; |
|||
p := 0; |
|||
repeat |
repeat |
||
isprime := true; |
|||
j := 0; |
|||
repeat |
repeat |
||
p +=1 |
|||
until pr[p]= 0; |
|||
j := (p+1)*p*2; |
|||
if j>MAXLIMIT then |
|||
BREAK; |
|||
d := 2*p+1; |
|||
repeat |
|||
pr[j] := 1; |
|||
j += d; |
|||
until j>MAXLIMIT; |
|||
until false; |
|||
SmallPrimes[1] := 3; |
|||
SmallPrimes[2] := 5; |
|||
j := 3; |
|||
d := 7; |
|||
flipflop := (2+1)-1;//7+2*2,11+2*1,13,17,19,23 |
|||
p := 3; |
|||
repeat |
|||
if pr[p] = 0 then |
|||
begin |
|||
SmallPrimes[j] := d; |
|||
inc(j); |
|||
end; |
|||
d += 2*flipflop; |
|||
p+=flipflop; |
|||
flipflop := 3-flipflop; |
|||
until (p > MAXLIMIT) OR (j>High(SmallPrimes)); |
|||
end; |
|||
function OutPots(pD:tpPrimeFac;n:NativeInt):Ansistring; |
|||
var |
|||
s: String[31]; |
|||
chk,p,i: NativeInt; |
|||
Begin |
|||
str(n,s); |
|||
result := s+' :'; |
|||
with pd^ do |
|||
begin |
|||
str(pfDivCnt:3,s); |
|||
result += s+' : '; |
|||
chk := 1; |
|||
For n := 0 to pfMaxIdx-1 do |
|||
Begin |
|||
if n>0 then |
|||
result += '*'; |
|||
p := pFpotPrim[n]; |
|||
chk *= p; |
|||
str(p,s); |
|||
result += s; |
|||
i := pfpotMax[n]; |
|||
if i >1 then |
|||
Begin |
Begin |
||
str(pfpotMax[n],s); |
|||
result += '^'+s; |
|||
repeat |
|||
chk *= p; |
|||
dec(i); |
|||
until i <= 1; |
|||
end; |
end; |
||
inc(j); |
|||
until false; |
|||
end; |
|||
p := pfRemain; |
|||
If p >1 then |
|||
Begin |
Begin |
||
str(p,s); |
|||
chk *= p; |
|||
result += '*'+s; |
|||
end; |
end; |
||
str(chk,s); |
|||
result += '_chk_'+s+'<'; |
|||
until pr > 1 shl 16 -1; |
|||
str(pfSumOfDivs,s); |
|||
result += '_SoD_'+s+'<'; |
|||
end; |
|||
end; |
end; |
||
function smplPrimeDecomp(n:Uint64):tprimeFac; |
|||
procedure PrimeFacOut(proper:Boolean=true); |
|||
var |
var |
||
i, |
pr,i,pot,fac,q :NativeUInt; |
||
Begin |
|||
begin |
|||
with |
with result do |
||
Begin |
Begin |
||
pfDivCnt := 1; |
|||
pfSumOfDivs := 1; |
|||
pfRemain := n; |
|||
pfMaxIdx := 0; |
|||
with pfPrims[i] do |
|||
pfpotPrim[0] := 1; |
|||
pfpotMax[0] := 0; |
|||
write(potPrim,'*') |
|||
else |
|||
if Not(ODD(n)) then |
|||
write(potPrim,'^',potMax,'*'); |
|||
begin |
|||
with pfPrims[k] do |
|||
pot := BsfQWord(n); |
|||
pfMaxIdx := 1; |
|||
pfpotPrim[0] := 2; |
|||
pfpotMax[0] := pot; |
|||
write(potPrim,'^',potMax); |
|||
n := n shr pot; |
|||
pfRemain := n; |
|||
writeln(' got ',pfDivCnt-1,' proper divisors with sum : ',pfSumOfDivs-pfNum) |
|||
pfSumOfDivs := (1 shl (pot+1))-1; |
|||
else |
|||
pfDivCnt := pot+1; |
|||
end; |
|||
i := 1; |
|||
while i < High(SmallPrimes) do |
|||
begin |
|||
pr := SmallPrimes[i]; |
|||
q := n DIV pr; |
|||
if pr > q then |
|||
BREAK; |
|||
if n = pr*q then |
|||
Begin |
|||
pfpotPrim[pfMaxIdx] := pr; |
|||
pot := 0; |
|||
fac := pr; |
|||
repeat |
|||
n := q; |
|||
q := n div pr; |
|||
pot+=1; |
|||
fac *= pr; |
|||
until n <> pr*q; |
|||
pfpotMax[pfMaxIdx] := pot; |
|||
pfDivCnt *= pot+1; |
|||
pfSumOfDivs *= (fac-1)DIV(pr-1); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
inc(i); |
|||
end; |
|||
pfRemain := n; |
|||
if n > 1 then |
|||
Begin |
|||
pfDivCnt *= 2; |
|||
pfSumOfDivs *= n+1 |
|||
end; |
|||
end; |
end; |
||
end; |
end; |
||
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt; |
|||
function DivCount(const primeDecomp:tprimeFac):NativeUInt;inline; |
|||
//n must be multiple of base |
|||
begin |
|||
var |
|||
result := primeDecomp.pfDivCnt; |
|||
q,r: Uint64; |
|||
end; |
|||
i : NativeInt; |
|||
Begin |
|||
fillchar(dgt,SizeOf(dgt),#0); |
|||
i := 0; |
|||
// dgtNum:= n; |
|||
n := n div base; |
|||
result := 0; |
|||
repeat |
|||
r := n; |
|||
q := n div base; |
|||
r -= q*base; |
|||
n := q; |
|||
dgt[i] := r; |
|||
inc(i); |
|||
until (q = 0); |
|||
result := 0; |
|||
function SumOfDiv(const primeDecomp:tprimeFac):NativeUInt;inline; |
|||
while (result<i) AND (dgt[result] = 0) do |
|||
begin |
|||
inc(result); |
|||
result := primeDecomp.pfSumOfDivs; |
|||
inc(result); |
|||
end; |
end; |
||
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt; |
|||
procedure PrimeDecomposition(n:Uint32;var res:tprimeFac); |
|||
var |
var |
||
q :NativeInt; |
|||
i,pr,fac,cnt,DivCnt,quot{to minimize divisions} : NativeUint; |
|||
Begin |
Begin |
||
result := 0; |
|||
q := dgt[result]+1; |
|||
if q = base then |
|||
repeat |
|||
dgt[result] := 0; |
|||
inc(result); |
|||
Begin |
|||
q := dgt[result]+1; |
|||
until q <> base; |
|||
dgt[result] := q; |
|||
result +=1; |
|||
end; |
|||
procedure SieveOneSieve(var pdf:tPrimeDecompField); |
|||
var |
|||
dgt:tDigits; |
|||
i, j, k,pr,fac,n : NativeUInt; |
|||
begin |
|||
n := pdfOfs; |
|||
//init |
|||
for i := 0 to SizePrDeFe-1 do |
|||
begin |
|||
with pdf[i] do |
|||
Begin |
Begin |
||
pfDivCnt := 1; |
|||
pfSumOfDivs := 1; |
|||
pfRemain := n+i; |
|||
pfMaxIdx := 0; |
|||
pfpotPrim[0] := 1; |
|||
pfpotMax[0] := 0; |
|||
end; |
end; |
||
end; |
|||
//first factor 2. Make n+i even |
|||
end |
|||
else |
|||
i := (pdfIdx+n) AND 1; |
|||
IF (n = 0) AND (pdfIdx<2) then |
|||
i := 2; |
|||
repeat |
repeat |
||
with pdf[i] do |
|||
begin |
|||
j := BsfQWord(n+i); |
|||
pfMaxIdx := 1; |
|||
pfpotPrim[0] := 2; |
|||
pfpotMax[0] := j; |
|||
pfRemain := (n+i) shr j; |
|||
pfSumOfDivs := (1 shl (j+1))-1; |
|||
pfDivCnt := j+1; |
|||
end; |
|||
i += 2; |
|||
until i >=SizePrDeFe; |
|||
// i now index in SmallPrimes |
|||
i := 0; |
|||
repeat |
|||
//search next prime that is in bounds of sieve |
|||
repeat |
|||
inc(i); |
|||
if i >= High(SmallPrimes) then |
|||
BREAK; |
|||
pr := SmallPrimes[i]; |
|||
k := pr-n MOD pr; |
|||
if (k = pr) AND (n>0) then |
|||
k:= 0; |
|||
if k < SizePrDeFe then |
|||
break; |
|||
until false; |
|||
if i >= High(SmallPrimes) then |
|||
BREAK; |
|||
//no need to use higher primes |
|||
if pr*pr > n+SizePrDeFe then |
|||
BREAK; |
|||
// j is power of prime |
|||
j := CnvtoBASE(dgt,n+k,pr); |
|||
IF pr*quot = n then |
|||
repeat |
|||
with pdf[k] do |
|||
Begin |
Begin |
||
pfpotPrim[pfMaxIdx] := pr; |
|||
pfpotMax[pfMaxIdx] := j; |
|||
pfDivCnt *= j+1; |
|||
fac := pr; |
|||
repeat |
|||
pfRemain := pfRemain DIV pr; |
|||
dec(j); |
|||
fac *= pr; |
|||
until j<= 0; |
|||
pfSumOfDivs *= (fac-1)DIV(pr-1); |
|||
inc(pfMaxIdx); |
|||
k += pr; |
|||
j := IncByBaseInBase(dgt,pr); |
|||
inc(Cnt); |
|||
end; |
end; |
||
until k >= SizePrDeFe; |
|||
until false; |
until false; |
||
//a big prime left over? |
|||
//correct sum of & count of divisors |
|||
IF n > 1 then |
|||
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(Cnt); |
|||
DivCnt *= 2; |
|||
end; |
end; |
||
end; |
|||
res.pfCnt:= cnt; |
|||
res.pfDivCnt := DivCnt; |
|||
res.pfSumOfDivs := 0; |
|||
end; |
end; |
||
procedure NextSieve; |
|||
procedure GetDivs(var primeDecomp:tprimeFac;var Divs:tDivisors); |
|||
begin |
|||
dec(pdfIDX,SizePrDeFe); |
|||
inc(pdfOfs,SizePrDeFe); |
|||
SieveOneSieve(PrimeDecompField); |
|||
end; |
|||
function GetNextPrimeDecomp:tpPrimeFac; |
|||
begin |
|||
if pdfIDX >= SizePrDeFe then |
|||
NextSieve; |
|||
result := @PrimeDecompField[pdfIDX]; |
|||
inc(pdfIDX); |
|||
end; |
|||
procedure Init_Sieve(n:NativeUint); |
|||
//Init Sieve pdfIdx,pdfOfs are Global |
|||
begin |
|||
pdfIdx := n MOD SizePrDeFe; |
|||
pdfOfs := n-pdfIdx; |
|||
SieveOneSieve(PrimeDecompField); |
|||
end; |
|||
procedure InsertSort(pDiv:tpDivisor; Left, Right : NativeInt ); |
|||
var |
|||
I, J: NativeInt; |
|||
Pivot : tItem; |
|||
begin |
|||
for i:= 1 + Left to Right do |
|||
begin |
|||
Pivot:= pDiv[i]; |
|||
j:= i - 1; |
|||
while (j >= Left) and (pDiv[j] > Pivot) do |
|||
begin |
|||
pDiv[j+1]:=pDiv[j]; |
|||
Dec(j); |
|||
end; |
|||
pDiv[j+1]:= pivot; |
|||
end; |
|||
end; |
|||
procedure GetDivisors(pD:tpPrimeFac;var Divs:tDivisors); |
|||
var |
var |
||
pDivs : tpDivisor; |
pDivs : tpDivisor; |
||
pPot : UInt64; |
|||
i,len,j,l,p,k: Int32; |
|||
Begin |
Begin |
||
i := |
i := pD^.pfDivCnt; |
||
IF i > Length(Divs) then |
|||
setlength(Divs,i); |
|||
pDivs := @Divs[0]; |
pDivs := @Divs[0]; |
||
pDivs[0] := 1; |
pDivs[0] := 1; |
||
len := 1; |
len := 1; |
||
l := |
l := 1; |
||
with pD^ do |
|||
Begin |
|||
For i := 0 to primeDecomp.pfCnt-1 do |
|||
For i := 0 to pfMaxIdx-1 do |
|||
with primeDecomp.pfPrims[i] do |
|||
begin |
|||
//Multiply every divisor before with the new primefactors |
//Multiply every divisor before with the new primefactors |
||
//and append them to the list |
//and append them to the list |
||
k := |
k := pfpotMax[i]; |
||
p := |
p := pfpotPrim[i]; |
||
pPot :=1; |
pPot :=1; |
||
repeat |
repeat |
||
Line 4,471: | Line 4,664: | ||
Begin |
Begin |
||
pDivs[l]:= pPot*pDivs[j]; |
pDivs[l]:= pPot*pDivs[j]; |
||
sum += pDivs[l]; |
|||
inc(l); |
inc(l); |
||
end; |
end; |
||
dec(k); |
dec(k); |
||
until k<0; |
until k<=0; |
||
len := l; |
len := l; |
||
end; |
end; |
||
p := pfRemain; |
|||
primeDecomp.pfSumOfDivs := sum; |
|||
If p >1 then |
|||
begin |
|||
For j := 0 to len-1 do |
|||
Begin |
|||
pDivs[l]:= p*pDivs[j]; |
|||
inc(l); |
|||
end; |
|||
len := l; |
|||
end; |
|||
end; |
|||
//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; |
end; |
||
procedure AllFacsOut(pD:tpPrimeFac;var Divs:tdivisors;proper:boolean=true); |
|||
Function GetDivisors(n:Uint32;var Divs:tDivisors):Int32; |
|||
var |
|||
i:Int32; |
|||
Begin |
|||
PrimeDecomposition(n,primeDecomp); |
|||
i := DivCount(primeDecomp); |
|||
IF i > Length(Divs) then |
|||
setlength(Divs,i+1); |
|||
GetDivs(primeDecomp,Divs); |
|||
result := DivCount(primeDecomp); |
|||
end; |
|||
procedure AllFacsOut(n: Uint32;Divs:tDivisors;proper:boolean=true); |
|||
var |
var |
||
k,j: Int32; |
k,j: Int32; |
||
Begin |
Begin |
||
k := |
k := pd.pfDivCnt-1; |
||
PrimeFacOut(proper); |
|||
IF proper then |
IF proper then |
||
dec(k); |
dec(k); |
||
IF k > 0 then |
IF k > 0 then |
||
Begin |
|||
For j := 0 to k-1 do |
For j := 0 to k-1 do |
||
write(Divs[j],','); |
write(Divs[j],','); |
||
writeln(Divs[k]); |
|||
end; |
|||
end; |
end; |
||
procedure SpeedTest(Limit:Uint32); |
|||
var |
var |
||
pPrimeDecomp :tpPrimeFac; |
|||
Ticks,SumDivCnt : Int64; |
|||
Mypd : tPrimeFac; |
|||
Divs:tDivisors; |
|||
T0:Int64; |
|||
n : NativeUInt; |
|||
Begin |
Begin |
||
Ticks := GetTickCount64; |
|||
SumDivCnt := 0; |
|||
For number := 1 to Limit do |
|||
inc(SumDivCnt,GetDivisors(number,Divisors)); |
|||
writeln('SpeedTest ',(GetTickCount64-Ticks)/1000:0:3,' secs for 1..',Limit); |
|||
writeln('mean count of divisors ',SumDivCnt/limit:0:3); |
|||
writeln; |
|||
end; |
|||
var |
|||
Divisors : tDivisors; |
|||
number: Int32; |
|||
BEGIN |
|||
InitSmallPrimes; |
InitSmallPrimes; |
||
setlength(Divisors,1); |
|||
SpeedTest(1000*1000); |
|||
T0 := GetTickCount64; |
|||
writeln('Enter a number between 1 and 4294967295: '); |
|||
n := 0; |
|||
writeln('3491888400 is a nice choice :'); |
|||
readln(number); |
|||
IF number >= 0 then |
|||
Begin |
|||
writeln('Proper number version'); |
|||
AllFacsOut(number,Divisors); |
|||
Init_Sieve(0); |
|||
writeln('including n version'); |
|||
repeat |
|||
AllFacsOut(number,Divisors,false); |
|||
pPrimeDecomp:= GetNextPrimeDecomp; |
|||
end; |
|||
// GetDivisors(pPrimeDecomp,Divs); |
|||
//https://en.wikipedia.org/wiki/Highly_composite_number <= HCN |
|||
inc(n); |
|||
//http://wwwhomes.uni-bielefeld.de/achim/highly.txt the first 1200 HCN |
|||
until n > 10*1000*1000+1; |
|||
END.</lang> |
|||
T0 := GetTickCount64-T0; |
|||
{out}} |
|||
writeln('runtime ',T0/1000:0:3,' s'); |
|||
<pre> |
|||
GetDivisors(pPrimeDecomp,Divs); |
|||
SpeedTest 0.296 secs for 1..1000000 |
|||
AllFacsOut(pPrimeDecomp,Divs,true); |
|||
AllFacsOut(pPrimeDecomp,Divs,false); |
|||
writeln('simple version'); |
|||
T0 := GetTickCount64; |
|||
n := 0; |
|||
repeatSpeedTest 0.296 secs for 1..1000000 |
|||
mean count of divisors 13.970 |
mean count of divisors 13.970 |
||
**SpeedTest 5.707 secs for 1..10000000 |
**SpeedTest 5.707 secs for 1..10000000 |
||
Line 4,565: | Line 4,739: | ||
including n version |
including n version |
||
123456789 = 3^2*3607*3803 got 12 divisors with sum : 178422816 |
123456789 = 3^2*3607*3803 got 12 divisors with sum : 178422816 |
||
1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789 |
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; |
|||
T0 := GetTickCount64-T0; |
|||
writeln('runtime ',T0/1000:0:3,' s'); |
|||
GetDivisors(@Mypd,Divs); |
|||
AllFacsOut(@Mypd,Divs,true); |
|||
AllFacsOut(@Mypd,Divs,false); |
|||
end. |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
runtime 1.216 s |
|||
1,11,909091 |
|||
1,11,909091,10000001 |
|||
simple version |
|||
runtime 5.854 s |
|||
1,11,909091 |
|||
1,11,909091,10000001 |
|||
//out-commented GetDivisors, but still calculates sum of divisors and count of divisors |
|||
// 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 4.922 s |
|||
1,11,909091 |
|||
1,11,909091,10000001</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |