%I #11 Mar 31 2018 06:44:58
%S 1,1,4,8,20,38,88,162,336,624,1211,2195,4109,7295,13190,23072,40618,
%T 69838,120486,204006,345595,577387,962961,1588483,2613930,4262138,
%U 6928799,11179251,17976330,28720552,45729595,72401921,114239033
%N G.f.: A(x) = exp( Sum_{n>=1} sigma(n^2)*x^n/n ), a power series in x with integer coefficients.
%C Compare with g.f. for partition numbers: exp( Sum_{n>=1} sigma(n)*x^n/n ), where sigma(n) = A000203(n) is the sum of the divisors of n.
%H Vaclav Kotesovec, <a href="/A156303/b156303.txt">Table of n, a(n) for n = 0..1000</a>
%F a(n) = (1/n)*Sum_{k=1..n} sigma(k^2) * a(n-k) for n>0, with a(0)=1.
%F Euler transform of Dedekind psi function, cf. A001615. - _Vladeta Jovovic_, Feb 12 2009
%F a(n) ~ exp(3^(4/3) * (5*Zeta(3))^(1/3) * n^(2/3) / (2*Pi)^(2/3) - Pi^(2/3) * n^(1/3) / (2^(4/3) * (15*Zeta(3))^(1/3)) + 1/12 - Pi^2 / (720 * Zeta(3))) * (5*Zeta(3))^(7/36) / (A * 2^(2/9) * 3^(11/36) * Pi^(29/36) * n^(25/36)), where A is the Glaisher-Kinkelin constant A074962. - _Vaclav Kotesovec_, Mar 24 2018
%e G.f.: A(x) = 1 + x + 4*x^2 + 8*x^3 + 20*x^4 + 38*x^5 + 88*x^6 +...
%e log(A(x)) = x + 7*x^2/2 + 13*x^3/3 + 31*x^4/4 + 31*x^5/5 + 127*x^6/6 +...
%t nmax = 40; CoefficientList[Series[Exp[Sum[Sum[k*Sum[MoebiusMu[d]^2 / d, {d, Divisors @ k}] * x^(j*k) / j, {k, 1, Floor[nmax/j] + 1}], {j, 1, nmax}]], {x, 0, nmax}], x] (* _Vaclav Kotesovec_, Mar 31 2018 *)
%o (PARI) {a(n)=polcoeff(exp(sum(m=1,n,sigma(m^2)*x^m/m)+x*O(x^n)),n)}
%o (PARI) {a(n)=if(n==0,1,(1/n)*sum(k=1,n,sigma(k^2)*a(n-k)))}
%Y Cf. A000203 (sigma), A000041 (partitions), A001615, A301594.
%K nonn
%O 0,3
%A _Paul D. Hanna_, Feb 08 2009