(SageMath)
z = PowerSeriesRing(ZZ, 'z').gen().O(33)
g = 1 + z + sqrt(1 - 2*z - 3*z**2)
f = -z * g.derivative() / g
print([1] + [sum(moebius(n // d) * f[d]
for d in divisors(n)) // n for n in range(1, 33)])
(Python)
from typing import Callable
from functools import cache
from math import comb
# Define 'binomial' for compatibility with Maple.
def binomial(n: int, k: int) -> int:
if 0 <= k <= n: return comb(n, k)
if k <= n < 0: return comb(-k-1, n-k)*(-1)**(n-k)
if n < 0 <= k: return comb(-n+k-1, k)*(-1)**k
return 0
def EulerInvTransform(f: Callable) -> Callable:
@cache
def h(n: int, k: int) -> int:
if n == 0: return 1
if k < 1: return 0
return sum(binomial(b(k)+j-1, j) * h(n-k*j, k-1)
for j in range(1 + n // k))
@cache
def b(n: int) -> int:
if n == 0: return f(0)
return f(n) - h(n, n - 1)
return b
print([a(n) for n in range(33)])
