OFFSET
0,3
COMMENTS
Euler transform of the octagonal numbers (A000567).
LINKS
Vaclav Kotesovec, Table of n, a(n) for n = 0..2000
M. Bernstein and N. J. A. Sloane, Some canonical sequences of integers, Linear Alg. Applications, 226-228 (1995), 57-72; erratum 320 (2000), 210; arXiv:math/0205301 [math.CO], 2002.
M. Bernstein and N. J. A. Sloane, Some canonical sequences of integers, Linear Alg. Applications, 226-228 (1995), 57-72; erratum 320 (2000), 210. [Link to Lin. Alg. Applic. version together with omitted figures]
N. J. A. Sloane, Transforms
Eric Weisstein's World of Mathematics, Octagonal Number
FORMULA
G.f.: Product_{k>=1} 1/(1 - x^k)^(k*(3*k-2)).
a(n) ~ exp(4*Pi*n^(3/4) / (3*5^(1/4)) - 2*Zeta(3) * sqrt(5*n) / Pi^2 - 10*Zeta(3)^2 * (5*n)^(1/4) / Pi^5 - 200*Zeta(3)^3 / (3*Pi^8) - 3*Zeta(3) / (4*Pi^2) - 1/6) * A^2 / (2^(3/2) * 5^(1/12) * Pi^(1/6) * n^(7/12)), where A is the Glaisher-Kinkelin constant A074962. - Vaclav Kotesovec, Nov 08 2017
MAPLE
with(numtheory):
a:= proc(n) option remember; `if`(n=0, 1, add(add(
d^2*(3*d-2), d=divisors(j))*a(n-j), j=1..n)/n)
end:
seq(a(n), n=0..35); # Alois P. Heinz, Dec 02 2016
MATHEMATICA
nmax=32; CoefficientList[Series[Product[1/(1 - x^k)^(k (3 k - 2)), {k, 1, nmax}], {x, 0, nmax}], x]
PROG
(Python)
from sympy import divisors
from sympy.core.cache import cacheit
@cacheit
def a(n): return 1 if n==0 else sum(sum(d**2*(3*d - 2) for d in divisors(j))*a(n - j) for j in range(1, n + 1))//n
print([a(n) for n in range(51)]) # Indranil Ghosh, Aug 06 2017, after Maple code
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
Ilya Gutkovskiy, Nov 30 2016
STATUS
approved