%I #33 May 09 2021 22:13:13
%S 1,2,6,42,104,120,165,245,272,561,1806,47058,765625,1137501,3874128,
%T 9131793,2214502422,52495396602
%N Numbers m such that (m'+1)' = m-1, where m' is the arithmetic derivative of m.
%C The differential equation whose solutions are the primary pseudoperfect numbers is m' = k*m-1, with k a positive integer. Let us rewrite the equation as m'+1 = k*m and then take the derivative: (m'+1)' = (k*m)' = k'*m + k*m' = k'*m + k*(k*m-1) = (k'+k^2)*m-k. Let k=1: (m'+1)' = m-1. The solutions of this equation are the primary pseudoperfect numbers plus couples of numbers (x,y) for which x' = y-1 and y' = x-1.
%C A054377 is a subsequence of this sequence.
%C a(17) > 10^9. - _Michel Marcus_, Nov 05 2014
%C a(19) > 10^11. - _Giovanni Resta_, Jun 04 2016
%e 765625' = 1137500; (1137500 + 1)' = 1137501' = 765624 = 765625 - 1, so 765625 is a term.
%e 1137501' = 765624; (765624 + 1)' = 765625' = 1137500 = 1137501 - 1, so 1137501 is a term.
%p with(numtheory);
%p P:=proc(i)
%p local a,n,p,pfs;
%p for n from 1 to i do
%p pfs:=ifactors(n)[2]; a:=n*add(op(2,p)/op(1,p),p=pfs);
%p pfs:=ifactors(a+1)[2]; a:=(a+1)*add(op(2,p)/op(1,p),p=pfs);
%p if a=n-1 then print(n); fi;
%p od;
%p end:
%p P(10000000);
%t A003415[n_]:=If[Abs[n]<2,0,n*Total[#2/#1&@@@FactorInteger[Abs[n]]]];
%t Select[Range[1,100000],A003415[A003415[#]+1]==#-1&] (* _Julien Kluge_, Jul 08 2016 *)
%o (PARI) ad(n) = sum(i=1, #f=factor(n)~, n/f[1, i]*f[2, i]);
%o isok(n) = ad(ad(n)+1) == n-1; \\ _Michel Marcus_, Nov 05 2014
%Y Cf. A054377, A203617.
%K nonn,more
%O 1,2
%A _Paolo P. Lava_, Jan 20 2012
%E a(17)-a(18) from _Giovanni Resta_, Jun 04 2016