A358451 Inverse Euler transform of the Riordan numbers, (A005043).

1, 0, 1, 1, 2, 5, 11, 28, 68, 174, 445, 1166, 3068, 8190, 21994, 59585, 162360, 445145, 1226376, 3394654, 9434260, 26317865, 73661588, 206809307, 582255448, 1643536725, 4650250254, 13186484316, 37468566744, 106666821221, 304200399505, 868977304140, 2486163857424

Offset: 0

Author

Peter Luschny, Nov 20 2022

Keywords

• nonn

Links

• OEIS Wiki, Euler transform

Crossrefs

Cf. A005043.

Programs

Maple

EulerInvTransform := proc(f) local c, b;
    c := proc(n) option remember;
    ifelse(n = 0, f(0), f(n) - b(n, n-1)) end:
    b := proc(n, k) option remember;
    if n = 0 then return 1 elif k < 1 then return 0 fi;
    add(binomial(c(k) + j - 1, j)*b(n-k*j, k-1), j=0..n/k) end:
c end:
a := EulerInvTransform(A005043): seq(a(n), n = 0..32);

Mathematica

EulerInvTransform[seq_List] := Module[{final = {}}, Do[AppendTo[final, i*seq[[i]] - Sum[final[[d]]*seq[[i-d]], {d, i-1}]], {i, Length[seq]}]; Table[Sum[MoebiusMu[i/d]*final[[d]], {d, Divisors[i]}]/i, {i, Length[seq] }]];
A005043[n_] := A005043[n] = If[n <= 1, 1-n, (n-1)*(2*A005043[n-1] + 3*A005043[n-2])/(n+1)];
Join[{1}, EulerInvTransform[Array[A005043, 32]]] (* Jean-François Alcover, Jun 15 2024 *)

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
a = EulerInvTransform(A005043)
print([a(n) for n in range(33)])

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)])

Formula

a(n) ~ 3^(n + 1/2) / (4*sqrt(Pi)*n^(3/2)). - Vaclav Kotesovec, Jun 15 2024