OFFSET
0,4
LINKS
Seiichi Manyama, Table of n, a(n) for n = 0..10000 (terms 0..1000 from Alois P. Heinz)
FORMULA
a(n) = (1/n)*Sum_{k=1..n} a(n-k)*b(k), n>0, a(0)=1, b(k)=Sum_{d|k} (-1)^(n/d+1)*d*A001055(d).
MAPLE
with(numtheory):
g:= proc(n, k) option remember; `if`(n>k, 0, 1)+
`if`(isprime(n), 0, add(`if`(d>k, 0, g(n/d, d)),
d=divisors(n) minus {1, n}))
end:
b:= proc(n) b(n):= add((-1)^(n/d+1)*d*g(d$2), d=divisors(n)) end:
a:= proc(n) a(n):= `if`(n=0, 1, add(a(n-k)*b(k), k=1..n)/n) end:
seq(a(n), n=0..60); # Alois P. Heinz, May 16 2014
MATHEMATICA
g[n_, k_] := g[n, k] = If[n > k, 0, 1] + If[PrimeQ[n], 0, Sum[If[d > k, 0, g[n/d, d]], {d, Divisors[n] ~Complement~ {1, n}}]];
b[n_] := b[n] = Sum[(-1)^(n/d + 1)*d*g[d, d], {d, Divisors[n]}];
a[n_] := a[n] = If[n == 0, 1, Sum[a[n - k]*b[k], {k, 1, n}]/n];
Table[a[n], {n, 0, 60}] (* Jean-François Alcover, Mar 23 2017, after Alois P. Heinz *)
PROG
(Python)
from sympy.core.cache import cacheit
from sympy import divisors, isprime
@cacheit
def g(n, k): return (0 if n>k else 1) + (0 if isprime(n) else sum(0 if d>k else g(n//d, d) for d in divisors(n)[1:-1]))
@cacheit
def b(n): return sum((-1)**(n//d + 1)*d*g(d, d) for d in divisors(n))
@cacheit
def a(n): return 1 if n==0 else sum(a(n - k)*b(k) for k in range(1, n + 1))//n
print([a(n) for n in range(61)]) # Indranil Ghosh, Aug 19 2017, after Maple code
CROSSREFS
KEYWORD
nonn
AUTHOR
Vladeta Jovovic, Jan 19 2002
STATUS
approved