OFFSET
1,1
COMMENTS
EXAMPLE
6’=5; 5’=1; 6=5+1 (k=2);
38’=21; 21’=10; 10’=7; 38=21+10+7 (k=3);
42’=41; 41’=1; 42=41+1 (k=2);
62’=33; 33’=14; 14’=9; 9'=6; 62=33+14+9+6 (k=4);
146’=75; 75’=55; 55’=16; 146=75+55+16 (k=3);
1145’=234; 234’=291; 291’=100; 100'=140; 140'=188; 188'=192; 1145=234+291+100+140+188+192 (k=6).
MAPLE
with(numtheory);
A216384:= proc(i)
local a, b, n, p, pfs;
for n from 1 to i do
pfs:=ifactors(n)[2]; a:=n*add(op(2, p)/op(1, p), p=pfs);
if a<n then b:=a;
while b<n do
pfs:=ifactors(a)[2]; a:=a*add(op(2, p)/op(1, p), p=pfs);
if a=0 then break; else b:=b+a; fi;
od;
if b=n then print(n); fi; fi; od;
end:
A216384 (10000000);
MATHEMATICA
d[1]=0; d[n_] := n*Total[#2/#1 & @@@ FactorInteger[n]]; seqQ[n_] := Module[{s = 0, k = n}, While[s < n && k > 0, k = d[k]; s += k]; k < n && s == n]; Select[ Range[16000], seqQ] (* Amiram Eldar, Mar 30 2019 *)
PROG
(Python)
from sympy import factorint
A216384 = []
for n in range(1, 10**5):
....ndsum = nd = sum([int(n*e/p) for p, e in factorint(n).items()])
....while ndsum <= n and nd > 1:
........nd = sum([int(nd*e/p) for p, e in factorint(nd).items()])
........ndsum += nd
........if ndsum == n:
............A216384.append(n)
# Chai Wah Wu, Aug 21 2014
(PARI) der(n) = sum(i=1, #f=factor(n)~, n/f[1, i]*f[2, i]);
isok(n) = {my(s = 0, kn = n, nb = 0, d); until (s == kn, d = der(n); if (d==0, return(0)); s += d; if (s > kn, return (0)); n = d; nb++; ); nb > 1; } \\ Michel Marcus, Mar 30 2019
CROSSREFS
KEYWORD
nonn,more
AUTHOR
Paolo P. Lava, Sep 06 2012
EXTENSIONS
a(20)-a(23) from Amiram Eldar, Mar 30 2019
STATUS
approved