|
%I #9 Sep 08 2022 08:46:26
%S 1,161,209,221,4265,12690,15941,22217,24041,25637,30377,38117,39077,
%T 48617,49097,55877,68441,73817,76457,80357,88457,95237,98117,99941,
%U 105641,110057,115397,122537,130217,131141,136517,143237,147941,148697,152357,154457,159077
%N Numbers k for which sigma(k) = k + k'', where k'' is the second derivative of k (A068346).
%C If p and q are different prime numbers and p + q is in A007850 (Giuga numbers) then m = p*q is a term because sigma(m) = sigma(p*q) = p*q + p + q + 1 and m + m'' = p*q + (p + q)' = p*q + p + q + 1 and sigma(m) = m + m''.
%e sigma(1) = 1 and 1 + 1'' = 1 so 1 is a term.
%e sigma(161) = 1 + 7 + 23 + 161 = 192 and 161 + 161'' = 161 + 30' = 161 + 31 = 192 so 161 is a term.
%e sigma(12690) = sigma(2*3^3*5*47) = 34560 and 12690 + 12690'' = 12690 + A068346(12690) = 12690 + 21870 = 34560 so 12690 is a term.
%t d[0] = d[1] = 0; d[n_] := n * Plus @@ ((Last[#]/First[#]) & /@ FactorInteger[n]); Select[Range[160000], DivisorSigma[1, #] == # + d[d[#]] &] (* _Amiram Eldar_, Oct 18 2021 *)
%o (Magma) f:=func<n|n le 1 select 0 else n*(&+[Factorisation(n)[i][2] / Factorisation(n)[i][1]: i in [1..#Factorisation(n)]])>; [n:n in [1..160000]| DivisorSigma(1,n) eq n+Floor(f(Floor(f(n))))];
%o (PARI) ad(n) = vecsum([n/f[1]*f[2]|f<-factor(n+!n)~]); \\ A003415
%o isok(k) = sigma(k) == k+ad(ad(k)); \\ _Michel Marcus_, Oct 18 2021
%Y Cf. A000203, A003415, A068346, A007850.
%K nonn
%O 1,2
%A _Marius A. Burtea_, Oct 18 2021
|
