OFFSET
1,2
COMMENTS
Euler transform of A001001.
LINKS
Alois P. Heinz, Table of n, a(n) for n = 1..9000
Tad White, Counting Free Abelian Actions, arXiv preprint arXiv:1304.2830, 2013
FORMULA
a(n) ~ exp(2^(7/4) * Pi^(3/2) * Zeta(3)^(1/4) * n^(3/4) / (3^(3/2) * 5^(1/4)) - sqrt(5*Zeta(3)*n) / (2^(3/2)*Pi) + (sqrt(Pi) * 5^(1/4) / (2^(15/4) * 3^(3/2) * Zeta(3)^(1/4)) - sqrt(3) * 5^(5/4) * Zeta(3)^(3/4) / (2^(15/4) * Pi^(7/2))) * n^(1/4) - 25*Zeta(3) / (16*Pi^6) + (5 - 2*Zeta(3)) / (192*Pi^2)) * Pi^(1/4) * Zeta(3)^(1/8) / (2^(13/8) * 3^(1/4) * 5^(1/8) * n^(5/8)). - Vaclav Kotesovec, Mar 26 2018
MAPLE
with(numtheory):
b:= proc(n) option remember; add(d*sigma(d), d=divisors(n)) end:
a:= proc(n) option remember; `if`(n=0, 1, add(add(d*
b(d), d=divisors(j)) *a(n-j), j=1..n)/n)
end:
seq(a(n), n=1..40); # Alois P. Heinz, Mar 06 2015
MATHEMATICA
b[n_] := b[n] = DivisorSum[n, #*DivisorSigma[1, #]&];
a[n_] := a[n] = If[n == 0, 1, Sum[DivisorSum[j, #*b[#]&]*a[n-j], {j, 1, n}] /n];
Array[a, 40] (* Jean-François Alcover, Mar 27 2017, after Alois P. Heinz *)
nmax = 40; Rest[CoefficientList[Series[Exp[Sum[Sum[Sum[d*DivisorSigma[1, 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 *)
PROG
(PARI) EulerT(v)={Vec(exp(x*Ser(dirmul(v, vector(#v, n, 1/n))))-1, -#v)}
seq(n)={my(v=vector(n, i, 1)); for(k=1, 2, v=dirmul(v, vector(n, i, i^k))); EulerT(v)} \\ Andrew Howroyd, May 09 2023
CROSSREFS
KEYWORD
nonn
AUTHOR
N. J. A. Sloane, Jun 08 2013
STATUS
approved