%I #13 Sep 20 2023 16:00:11
%S 1,6,8,22,12,48,16,66,41,72,24,176,28,96,96,178,36,246,40,264,128,144,
%T 48,528,97,168,176,352,60,576,64,450,192,216,192,902,76,240,224,792,
%U 84,768,88,528,492,288,96,1424,177,582,288,616,108,1056,288,1056,320
%N Dirichlet convolution of sigma with Dedekind psi function.
%F a(n) = Sum{d|n} A000203(d) * A001615(n/d).
%F a(p) = A365647(p) + 2 = A365648(p) + 2 where p is a term of A000040.
%F Multiplicative with a(p^e) = (2 + ((e+1)*p^2 - 2*p - e - 1)*p^e)/(p-1)^2. - _Amiram Eldar_, Sep 15 2023
%F From _Vaclav Kotesovec_, Sep 15 2023: (Start)
%F Dirichlet g.f.: zeta(s)^2 * zeta(s-1)^2 / zeta(2*s).
%F Sum_{k=1..n} a(k) ~ 5*n^2 * (log(n)/4 + gamma/2 - 1/8 + 3*zeta'(2)/Pi^2 - 45*zeta'(4)/Pi^4), where gamma is the Euler-Mascheroni constant A001620. (End)
%t f[p_, e_] := (2 + ((e + 1)*p^2 - 2*p - e - 1)*p^e)/(p - 1)^2; a[1] = 1; a[n_] := Times @@ f @@@ FactorInteger[n]; Array[a, 100] (* _Amiram Eldar_, Sep 15 2023 *)
%o (PARI) for(n=1, 100, print1(direuler(p=2, n, (1+X) / ((1-X) * (1 - p*X)^2))[n], ", ")) \\ _Vaclav Kotesovec_, Sep 15 2023
%o (Python)
%o from sympy import divisors, primefactors, prod, reduced_totient, divisor_sigma
%o def psi(n):
%o return n*prod(p+1 for p in primefactors(n))//prod(primefactors(n))
%o def a(n): return sum(divisor_sigma(d, 1) * psi(n//d) for d in divisors(n))
%Y Cf. A000040, A000203, A001615, A365647, A365648.
%K nonn,easy,mult
%O 1,2
%A _Torlach Rush_, Sep 14 2023