Practical numbers: Difference between revisions

Content added Content deleted
m (→‎{{header|alternative}}: only memorizing heigher divisors)
m (Fix pascal version to run in delphi)
Line 138: Line 138:
{$ELSE}
{$ELSE}
{$APPTYPE CONSOLE}
{$APPTYPE CONSOLE}

{$ENDIF}
{$ENDIF}

uses
uses
sysutils;
sysutils
{$IFNDEF FPC}
,Windows
{$ENDIF}
;

const
LOW_DIVS = 0;
HIGH_DIVS = 2048 - 1;

type
type
tdivs = record
tdivs = record
DivsVal : array[0..2048-1] of Uint32;
DivsVal: array[LOW_DIVS..HIGH_DIVS] of Uint32;
⚫
DivsMaxIdx, DivsNum, DivsSumProp: NativeUInt;
DivsMaxIdx,
⚫
end;
DivsNum,
⚫
DivsSumProp :NativeUInt;
⚫
end;


var
var
Divs: tDivs;
Divs: tDivs;
HasSum :array of byte;
HasSum: array of byte;

procedure GetDivisors(var Divs:tdivs;n:Uint32);
procedure GetDivisors(var Divs: tdivs; n: Uint32);
//calc all divisors,keep sorted
//calc all divisors,keep sorted
var
var
i,quot,ug,og : UInt32;
i, quot, ug, og: UInt32;
sum: UInt64;
sum: UInt64;
begin
Begin
with Divs do
with Divs do
Begin
begin
DivsNum := n;
DivsNum := n;
sum := 0;
sum := 0;
ug := Low(tdivs.DivsVal);
ug := 0;
og := High(tdivs.DivsVal);
og := HIGH_DIVS;
i := 1;
i := 1;

while i*i < n do
while i * i < n do
begin
begin
quot := n div i;
quot := n div i;
IF n - quot*i = 0 then
if n - quot * i = 0 then
Begin
begin
DivsVal[og] := quot;
DivsVal[og] := quot;
Divs.DivsVal[ug] := i;
Divs.DivsVal[ug] := i;
inc(sum,quot+i);
inc(sum, quot + i);
dec(og);
dec(og);
inc(ug);
inc(ug);
Line 178: Line 189:
inc(i);
inc(i);
end;
end;
If i*i = n then
if i * i = n then
Begin
begin
DivsVal[og] := i;
DivsVal[og] := i;
inc(sum,i);
inc(sum, i);
dec(og);
dec(og);
end;
end;
//move higher divisors down
//move higher divisors down
while og < High(tDivs.DivsVal) do
while og < high_DIVS do
begin
begin
inc(og);
inc(og);
Line 191: Line 202:
inc(ug);
inc(ug);
end;
end;
DivsMaxIdx := ug-2;
DivsMaxIdx := ug - 2;
DivsSumProp := sum-n;
DivsSumProp := sum - n;
end;//with
end; //with
end;
end;


function SumAllSetsForPractical(Limit:Uint32):boolean;
function SumAllSetsForPractical(Limit: Uint32): boolean;
//mark sum and than shift by next divisor == add
//mark sum and than shift by next divisor == add
//for practical numbers every sum must be marked
//for practical numbers every sum must be marked
var
var
hs0,hs1 : pByte;
hs0, hs1: pByte;
idx,j,loLimit,maxlimit,delta : NativeUint;
idx, j, loLimit, maxlimit, delta: NativeUint;
begin
Begin
Limit := trunc(Limit*(Limit/Divs.DivsSumProp));
Limit := trunc(Limit * (Limit / Divs.DivsSumProp));
LoLimit:=0;
loLimit := 0;
maxlimit := 0;
maxlimit := 0;
hs0 := @HasSum[0];
hs0 := @HasSum[0];
hs0[0] := 1;//empty set
hs0[0] := 1; //empty set
for idx := 0 to Divs.DivsMaxIdx do
for idx := 0 to Divs.DivsMaxIdx do
Begin
begin
delta := Divs.DivsVal[idx];
delta := Divs.DivsVal[idx];
hs1 := @hs0[delta];
hs1 := @hs0[delta];
For j := maxlimit downto 0 do
for j := maxlimit downto 0 do
hs1[j] := hs1[j] or hs0[j];
hs1[j] := hs1[j] or hs0[j];
maxlimit += delta;
maxlimit := maxlimit + delta;
while hs0[LoLimit]<> 0 do
while hs0[loLimit] <> 0 do
inc(LoLimit);
inc(loLimit);
//IF there is a 0 < delta, it will never be set
//IF there is a 0 < delta, it will never be set
//IF there are more than the Limit set,
//IF there are more than the Limit set,
//it will be copied by the following Delta's
//it will be copied by the following Delta's
if (LoLimit < delta) OR (LoLimit > Limit) then
if (loLimit < delta) or (loLimit > Limit) then
Break;
Break;
end;
end;
result := (LoLimit > Limit);
result := (loLimit > Limit);
end;
end;


function isPractical(n:Uint32):boolean;
function isPractical(n: Uint32): boolean;
var
var
i:NativeInt;
i: NativeInt;
sum:NativeUInt;
sum: NativeUInt;
begin
Begin
if n= 1 then
if n = 1 then
EXIT(True);
EXIT(True);
if ODD(n) then
if ODD(n) then
EXIT(false);
EXIT(false);
if (n > 2) AND NOT((n MOD 4 = 0) or (n Mod 6 = 0)) then
if (n > 2) and not ((n mod 4 = 0) or (n mod 6 = 0)) then
EXIT(false);
EXIT(false);


GetDivisors(Divs,n);
GetDivisors(Divs, n);
i := n-1;
i := n - 1;
sum := Divs.DivsSumProp;
sum := Divs.DivsSumProp;
if sum < i then
if sum < i then
result := false
result := false
else
else
Begin
begin
IF length(HasSum) > sum+1+1 then
if length(HasSum) > sum + 1 + 1 then
FillChar(HasSum[0],sum+1,#0)
FillChar(HasSum[0], sum + 1, #0)
else
else
Begin
begin
setlength(HasSum,0);
setlength(HasSum, 0);
setlength(HasSum,sum+8+1);
setlength(HasSum, sum + 8 + 1);
end;
end;
result:=SumAllSetsForPractical(i);
result := SumAllSetsForPractical(i);
end;
end;
end;
end;


procedure OutIsPractical(n:nativeInt);
procedure OutIsPractical(n: nativeInt);
begin
begin
If isPractical(n) then
if isPractical(n) then
writeln(n,' is practical')
writeln(n, ' is practical')
else
else
writeln(n,' is not practical');
writeln(n, ' is not practical');
end;
end;


const
const
ColCnt = 10;
ColCnt = 10;
MAX = 333;
MAX = 333;

var
var
T0 : Int64;
T0: Int64;
n,col,count : NativeInt;
n, col, count: NativeInt;

Begin
begin
col:=ColCnt;
col := ColCnt;
count := 0;
count := 0;
For n := 1 to MAX do
for n := 1 to MAX do
if isPractical(n) then
if isPractical(n) then
begin
begin
write(n:5);
write(n: 5);
inc(count);
inc(count);
dec(col);
dec(col);
if col = 0 then
if col = 0 then
Begin
begin
writeln;
writeln;
col :=ColCnt;
col := ColCnt;
end;
end;
end;
end;
writeln;
writeln;
writeln('There are ',count,' pratical numbers from 1 to ',MAX);
writeln('There are ', count, ' pratical numbers from 1 to ', MAX);
writeln;
writeln;


Line 297: Line 310:
OutIsPractical(5384);
OutIsPractical(5384);
OutIsPractical(1441440);
OutIsPractical(1441440);
writeln(Divs.DivsNum,' has ',Divs.DivsMaxIdx+1,' proper divisors');
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors');
writeln((GetTickCount64-T0)/1000:10:3,' s');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
T0 := GetTickCount64;
T0 := GetTickCount64;
OutIsPractical(99998640);
OutIsPractical(99998640);
writeln(Divs.DivsNum,' has ',Divs.DivsMaxIdx+1,' proper divisors ');
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors ');
writeln((GetTickCount64-T0)/1000:10:3,' s');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
T0 := GetTickCount64;
T0 := GetTickCount64;
OutIsPractical(99998640);
OutIsPractical(99998640);
writeln(Divs.DivsNum,' has ',Divs.DivsMaxIdx+1,' proper divisors ');
writeln(Divs.DivsNum, ' has ', Divs.DivsMaxIdx + 1, ' proper divisors ');
writeln((GetTickCount64-T0)/1000:10:3,' s');
writeln((GetTickCount64 - T0) / 1000: 10: 3, ' s');
setlength(HasSum,0);
setlength(HasSum, 0);
{$IFNDEF UNIX} readln; {$ENDIF}
⚫
end.</lang>
end.
⚫
</lang>
{{out}}
{{out}}
<pre> TIO.RUN.
<pre> TIO.RUN.
Retrieved from "https://rosettacode.org/wiki/Practical_numbers"