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

%I #9 Jun 15 2024 05:55:57

%S 1,0,1,1,2,5,11,28,68,174,445,1166,3068,8190,21994,59585,162360,

%T 445145,1226376,3394654,9434260,26317865,73661588,206809307,582255448,

%U 1643536725,4650250254,13186484316,37468566744,106666821221,304200399505,868977304140,2486163857424

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

%H OEIS Wiki, <a href="https://oeis.org/wiki/Euler_transform">Euler transform</a>

%F a(n) ~ 3^(n + 1/2) / (4*sqrt(Pi)*n^(3/2)). - _Vaclav Kotesovec_, Jun 15 2024

%p EulerInvTransform := proc(f) local c, b;

%p c := proc(n) option remember;

%p ifelse(n = 0, f(0), f(n) - b(n, n-1)) end:

%p b := proc(n, k) option remember;

%p if n = 0 then return 1 elif k < 1 then return 0 fi;

%p add(binomial(c(k) + j - 1, j)*b(n-k*j, k-1), j=0..n/k) end:

%p c end:

%p a := EulerInvTransform(A005043): seq(a(n), n = 0..32);

%t 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] }]];

%t A005043[n_] := A005043[n] = If[n <= 1, 1-n, (n-1)*(2*A005043[n-1] + 3*A005043[n-2])/(n+1)];

%t Join[{1}, EulerInvTransform[Array[A005043, 32]]] (* _Jean-François Alcover_, Jun 15 2024 *)

%o (SageMath)

%o z = PowerSeriesRing(ZZ, 'z').gen().O(33)

%o g = 1 + z + sqrt(1 - 2*z - 3*z**2)

%o f = -z * g.derivative() / g

%o print([1] + [sum(moebius(n // d) * f[d]

%o for d in divisors(n)) // n for n in range(1, 33)])

%o (Python)

%o from typing import Callable

%o from functools import cache

%o from math import comb

%o # Define 'binomial' for compatibility with Maple.

%o def binomial(n: int, k: int) -> int:

%o if 0 <= k <= n: return comb(n, k)

%o if k <= n < 0: return comb(-k-1, n-k)*(-1)**(n-k)

%o if n < 0 <= k: return comb(-n+k-1, k)*(-1)**k

%o return 0

%o def EulerInvTransform(f: Callable) -> Callable:

%o @cache

%o def h(n: int, k: int) -> int:

%o if n == 0: return 1

%o if k < 1: return 0

%o return sum(binomial(b(k)+j-1, j) * h(n-k*j, k-1)

%o for j in range(1 + n // k))

%o @cache

%o def b(n: int) -> int:

%o if n == 0: return f(0)

%o return f(n) - h(n, n - 1)

%o return b

%o a = EulerInvTransform(A005043)

%o print([a(n) for n in range(33)])

%Y Cf. A005043.

%K nonn

%O 0,5

%A _Peter Luschny_, Nov 20 2022