|
%I #34 May 03 2020 13:42:54
%S 1,1,1,2,3,4,6,8,11,16,21,27,38,49,63,84,109,138,180,228,289,369,463,
%T 578,732,911,1128,1407,1741,2140,2646,3243,3968,4862,5925,7198,8770,
%U 10620,12833,15524,18718,22502,27075,32467,38873,46537,55565,66220,78946
%N Expansion of Product_{k>=1} (1+x^k)^A001055(k).
%H Seiichi Manyama, <a href="/A066806/b066806.txt">Table of n, a(n) for n = 0..10000</a> (terms 0..1000 from Alois P. Heinz)
%F 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).
%p with(numtheory):
%p g:= proc(n, k) option remember; `if`(n>k, 0, 1)+
%p `if`(isprime(n), 0, add(`if`(d>k, 0, g(n/d, d)),
%p d=divisors(n) minus {1, n}))
%p end:
%p b:= proc(n) b(n):= add((-1)^(n/d+1)*d*g(d$2), d=divisors(n)) end:
%p a:= proc(n) a(n):= `if`(n=0, 1, add(a(n-k)*b(k), k=1..n)/n) end:
%p seq(a(n), n=0..60); # _Alois P. Heinz_, May 16 2014
%t 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}}]];
%t b[n_] := b[n] = Sum[(-1)^(n/d + 1)*d*g[d, d], {d, Divisors[n]}];
%t a[n_] := a[n] = If[n == 0, 1, Sum[a[n - k]*b[k], {k, 1, n}]/n];
%t Table[a[n], {n, 0, 60}] (* _Jean-François Alcover_, Mar 23 2017, after _Alois P. Heinz_ *)
%o (Python)
%o from sympy.core.cache import cacheit
%o from sympy import divisors, isprime
%o @cacheit
%o 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]))
%o @cacheit
%o def b(n): return sum((-1)**(n//d + 1)*d*g(d, d) for d in divisors(n))
%o @cacheit
%o def a(n): return 1 if n==0 else sum(a(n - k)*b(k) for k in range(1, n + 1))//n
%o print([a(n) for n in range(61)]) # _Indranil Ghosh_, Aug 19 2017, after Maple code
%Y Cf. A000041, A001055, A066739, A057567, A321460.
%K nonn
%O 0,4
%A _Vladeta Jovovic_, Jan 19 2002
|
