96 lines
3 KiB
Text
96 lines
3 KiB
Text
% This program uses the 'bigint' cluster from PCLU's 'misc.lib'
|
|
|
|
% Remove leading and trailing whitespace (bigint$unparse adds a lot)
|
|
strip = proc (s: string) returns (string)
|
|
ac = array[char]
|
|
sc = sequence[char]
|
|
cs: ac := string$s2ac(s)
|
|
while ~ac$empty(cs) cand ac$bottom(cs)=' ' do ac$reml(cs) end
|
|
while ~ac$empty(cs) cand ac$top(cs)=' ' do ac$remh(cs) end
|
|
% There's a bug in ac2s that makes it not return all elements
|
|
% This is a workaround
|
|
return(string$sc2s(sc$a2s(cs)))
|
|
end strip
|
|
|
|
divisor_sum = proc (n: bigint) returns (bigint)
|
|
own zero: bigint := bigint$i2bi(0)
|
|
own one: bigint := bigint$i2bi(1)
|
|
own two: bigint := bigint$i2bi(2)
|
|
own three: bigint := bigint$i2bi(3)
|
|
|
|
total: bigint := one
|
|
power: bigint := two
|
|
while n//two=zero do
|
|
total := total + power
|
|
power := power * two
|
|
n := n / two
|
|
end
|
|
p: bigint := three
|
|
while p*p <= n do
|
|
sum: bigint := one
|
|
power := p
|
|
while n//p = zero do
|
|
sum := sum + power
|
|
power := power * p
|
|
n := n/p
|
|
end
|
|
total := total * sum
|
|
p := p + two
|
|
end
|
|
if n>one then total := total * (n+one) end
|
|
return(total)
|
|
end divisor_sum
|
|
|
|
classify_aliquot_sequence = proc (n: bigint)
|
|
LIMIT = 16
|
|
abi = array[bigint]
|
|
own zero: bigint := bigint$i2bi(0)
|
|
po: stream := stream$primary_output()
|
|
|
|
terms: array[bigint] := abi$predict(0,LIMIT)
|
|
abi$addh(terms, n)
|
|
|
|
classification: string := "non-terminating"
|
|
for i: int in int$from_to(1, limit-1) do
|
|
abi$addh(terms, divisor_sum(abi$top(terms)) - abi$top(terms))
|
|
if abi$top(terms) = n then
|
|
if i=1 then classification := "perfect"
|
|
elseif i=2 then classification := "amicable"
|
|
else classification := "sociable"
|
|
end
|
|
break
|
|
end
|
|
j: int := 1
|
|
while j<i cand terms[i] ~= terms[i-j] do j := j+1 end
|
|
if j<i then
|
|
if j=1 then classification := "aspiring"
|
|
else classification := "cyclic"
|
|
end
|
|
break
|
|
end
|
|
if abi$top(terms) = zero then
|
|
classification := "terminating"
|
|
break
|
|
end
|
|
end
|
|
|
|
stream$puts(po, strip(bigint$unparse(n)) || ": " || classification || ", sequence: "
|
|
|| strip(bigint$unparse(terms[0])))
|
|
for i: int in int$from_to(1, abi$high(terms)) do
|
|
if terms[i] = terms[i-1] then break end
|
|
stream$puts(po, " " || strip(bigint$unparse(terms[i])))
|
|
end
|
|
stream$putl(po, "")
|
|
end classify_aliquot_sequence
|
|
|
|
start_up = proc ()
|
|
for i: int in int$from_to(1, 10) do
|
|
classify_aliquot_sequence(bigint$i2bi(i))
|
|
end
|
|
for i: int in array[int]$elements(array[int]$
|
|
[11,12,28,496,220,1184,12496,1264460,790,909,562,1064,1488]) do
|
|
classify_aliquot_sequence(bigint$i2bi(i))
|
|
end
|
|
classify_aliquot_sequence(bigint$parse("15355717786080"))
|
|
classify_aliquot_sequence(bigint$parse("153557177860800"))
|
|
end start_up
|