# First read 110 terms of
A000081 into array b1
M:=100;
t1:=1/mul((1-x*y^i)^b1[i+1], i=2..M):
t2:=series(t1, y, M):
t3:=series(t2, x, M):
a:=(n, k)->coeff(coeff(t3, x, k), y, n);
g:=n->add(a(n-1+i, i), i=1..n-1);
# Alternative:
g:= proc(n, i) option remember; `if`(n=0, 1, `if`(i<1, 0,
add(binomial(T(i-1)+j-1, j) *g(n-i*j, i-1), j=0..n/i)))
end:
T:= n-> g(n, n):
b:= proc(n, i) option remember; `if`(n=0, 1, `if`(i<1, 0,
add(binomial(T(i)+j-1, j) *b(n-i*j, i-1), j=0..n/i)))
end:
a:= n-> b(n, n):
# Alternative:
g:= proc(n) option remember; `if`(n<=1, n, (add(add(d*
g(d), d=numtheory[divisors](j))*g(n-j), j=1..n-1))/(n-1))
end:
a:= proc(n) option remember; `if`(n=0, 1, add(add(d*
g(d+1), d=numtheory[divisors](j))*a(n-j), j=1..n)/n)
end: