A381889 - OEIS

%I #16 Feb 09 2026 06:38:08

%S 1,4,18,86,428,2192,11468,60986,328532,1788368,9819128,54302712,

%T 302157424,1690193728,9497996152,53588976802,303434431108,

%U 1723578967056,9818195961512,56071829010968,320970950634288,1841213871449152,10582333064327824,60929582362628968,351385363433883472

%N Expansion of e.g.f.: (BesselI(0, 2*x) + BesselI(1, 2*x))^2*exp(2*x).

%C Binomial transform of A151093.

%C For p prime, a(p) - 2 == 0 (mod 2*p).

%F a(n) = Sum_{k=0..n} binomial(n, k)*2^(n-k)*A005566(k).

%F a(n) = Sum_{k=0..n} binomial(n, k)*A001405(n-k)*A001700(k).

%F a(n) = Sum_{k=0..n} binomial(n, k)*A005773(n-k+1)*A005773(k+1). - _Mélika Tebni_, Mar 19 2025

%F a(n) ~ 6^(n+1) / (Pi*n). - _Vaclav Kotesovec_, Feb 09 2026

%p a := n-> add(binomial(n, k)*binomial(n-k, iquo(n-k,2))*binomial(2*k+1,k+1), k = 0 .. n): seq(a(n), n = 0 .. 24);

%t len := 24; Table[n!,{n, 0, len}] CoefficientList[Series[(BesselI[0, 2x] + BesselI[1, 2x])^2 Exp[2x], {x, 0, len}], x]  (* _Peter Luschny_, Mar 19 2025 *)

%o (Python)

%o from math import comb as C

%o def a(n):

%o     return sum(C(n, k)*2**(n-k)*C(k, k//2)*C(k+1, (k+1)//2) for k in range(n+1))

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

%o (PARI) my(x='x+O('x^30)); Vec(serlaplace((besseli(0, 2*x) + x*besseli(1, 2*x))^2*exp(2*x))) \\ _Michel Marcus_, Mar 11 2025

%Y Cf. A001405, A001700, A005566, A151093.

%K nonn

%O 0,2

%A _Mélika Tebni_, Mar 09 2025