OFFSET
0,3
LINKS
Alois P. Heinz, Table of n, a(n) for n = 0..1000
FORMULA
G.f.: Product_{j>=1} 1/(1-x^j)^A005817(j).
a(n) ~ c * 4^n / n^3, where c = 19.002514794... - Vaclav Kotesovec, May 30 2019
MAPLE
g:= proc(n) option remember; `if`(n<2, 1, (4*(2*n+3)*
g(n-1)+16*(n-1)*n*g(n-2))/((n+3)*(n+4)))
end:
a:= proc(n) option remember; `if`(n=0, 1, add(add(g(d)
*d, d=numtheory[divisors](j))*a(n-j), j=1..n)/n)
end:
seq(a(n), n=0..35);
MATHEMATICA
g[n_] := g[n] = If[n<2, 1, (4*(2*n+3)*g[n-1] + 16*(n-1)*n*g[n-2])/((n+3)* (n+4))];
a[n_] := a[n] = If[n == 0, 1, Sum[Sum[g[d]*d, {d, Divisors[j]}]*a[n-j], {j, 1, n}]/n];
Table[a[n], {n, 0, 35}] (* Jean-François Alcover, May 30 2019, after Alois P. Heinz *)
PROG
(Python)
from sympy.core.cache import cacheit
from sympy import divisors
@cacheit
def g(n): return 1 if n<2 else (4*(2*n + 3)*g(n - 1) + 16*(n - 1)*n*g(n - 2))//((n + 3)*(n + 4))
@cacheit
def a(n): return 1 if n==0 else sum(sum(g(d)*d for d in divisors(j))*a(n - j) for j in range(1, n + 1))//n
print([a(n) for n in range(36)]) # Indranil Ghosh, Oct 15 2017
CROSSREFS
KEYWORD
nonn
AUTHOR
Alois P. Heinz, Oct 15 2017
STATUS
approved