OFFSET
0,3
COMMENTS
Euler transform of the absolute values of A008683. - Tilman Neumann, Dec 13 2008
Euler transform of A008966. - Vaclav Kotesovec, Mar 31 2018
LINKS
Alois P. Heinz, Table of n, a(n) for n = 0..10000
FORMULA
G.f.: 1/Product_{k>0} (1-x^A005117(k)).
a(n) = 1/n*Sum_{k=1..n} A048250(k)*a(n-k).
G.f.: 1 + Sum_{i>=1} mu(i)^2*x^i / Product_{j=1..i} (1 - mu(j)^2*x^j). - Ilya Gutkovskiy, Jun 05 2017
a(n) ~ exp(2*sqrt(n)) / (4*Pi^(3/2)*n^(1/4)). - Vaclav Kotesovec, Mar 24 2018
MAPLE
with(numtheory):
a:= proc(n) option remember; `if`(n=0, 1, add(add(d*
abs(mobius(d)), d=divisors(j)) *a(n-j), j=1..n)/n)
end:
seq(a(n), n=0..60); # Alois P. Heinz, Mar 05 2015
MATHEMATICA
Table[Length[Select[Boole /@ Thread /@ SquareFreeQ /@ IntegerPartitions[n], FreeQ[#, 0] &]], {n, 48}] (* Jayanta Basu, Jul 02 2013 *)
a[n_] := a[n] = If[n==0, 1, Sum[Sum[d*Abs[MoebiusMu[d]], {d, Divisors[j]}] * a[n-j], {j, 1, n}]/n]; Table[a[n], {n, 0, 60}] (* Jean-François Alcover, Oct 10 2015, after Alois P. Heinz *)
nmax = 60; CoefficientList[Series[Exp[Sum[Sum[Abs[MoebiusMu[k]] * x^(j*k) / j, {k, 1, Floor[nmax/j] + 1}], {j, 1, nmax}]], {x, 0, nmax}], x] (* Vaclav Kotesovec, Mar 31 2018 *)
PROG
(Haskell)
a073576 = p a005117_list where
p _ 0 = 1
p ks'@(k:ks) m = if m < k then 0 else p ks' (m - k) + p ks m
-- Reinhard Zumkeller, Jun 01 2015
(Python)
from functools import lru_cache
from sympy import mobius, divisors
@lru_cache(maxsize=None)
def A073576(n): return sum(sum(d*abs(mobius(d)) for d in divisors(i, generator=True))*A073576(n-i) for i in range(1, n+1))//n if n else 1 # Chai Wah Wu, Aug 23 2024
CROSSREFS
KEYWORD
nonn
AUTHOR
Vladeta Jovovic, Aug 27 2002
STATUS
approved