|
%I #14 Sep 01 2023 03:57:58
%S 6,75,133,1005,1603,4258,5299,84292,89944,170568,192901,303003,695364,
%T 1633303
%N Consider the unitary aliquot parts, in ascending order, of a composite number. Take their sum and repeat the process deleting the minimum number and adding the previous sum. The sequence lists the numbers that after some iterations reach a sum equal to the reverse of themselves.
%H Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/UnitaryDivisor.html">Unitary Divisor</a>
%H Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/UnitaryDivisorFunction.html">Unitary Divisor Function</a>
%H Eric Weisstein's World of Mathematics, <a href="https://mathworld.wolfram.com/UnitaryPerfectNumber.html">Unitary Perfect Number</a>
%H Wikipedia, <a href="https://en.wikipedia.org/wiki/Unitary_divisor">Unitary divisor</a>
%e Unitary aliquot parts of 6 are 1, 2, 3. We have: 1 + 2 + 3 = 6 that is equal to its reverse.
%e Unitary aliquot parts of 75 are 1, 3, 25. We have: 1 + 3 + 25 = 29; 3 + 25 + 29 = 57 that is the reverse of 75.
%e Unitary aliquot parts of 84292 are 1, 4, 13, 52, 1621, 6484, 21073. We have: 1 + 4 + 13 + 52 + 1621 + 6484 + 21073 = 29248 that is the reverse of 84292.
%p with(numtheory): R:=proc(w) local x, y; x:=w; y:=0;while x>0 do
%p y:=10*y+(x mod 10); x:=trunc(x/10); od: y; end:
%p P:=proc(q, h) local a,b,c,k,n,t,v; v:=array(1..h);
%p for n from 1 to q do if not isprime(n) then a:=sort([op(divisors(n))]);
%p b:=[]; c:=ilog10(n)+1; for k from 1 to nops(a)-1 do if gcd(a[k],n/a[k])=1
%p then b:=[op(b),a[k]]; fi; od; if nops(b)>1 then
%p for k from 1 to nops(b) do v[k]:=b[k]; od; t:=nops(b)+1; v[t]:=add(v[k],k=1..nops(b)); if R(v[t])=n then print(n); else
%p while ilog10(v[t])+1<=c do t:=t+1; v[t]:=add(v[k], k=t-nops(b)..t-1);
%p if R(v[t])=n then print(n); break; fi; od; fi; fi; fi; od;
%p end: P(10^9,1000);
%Y Cf. A246544, A247012, A258142.
%K nonn,more,base
%O 1,1
%A _Paolo P. Lava_, May 25 2015
| 