%I #15 Jul 01 2021 03:35:28

%S 1,0,1,2,12,64,485,4038,39991,441992,5492322,75171700,1127989577,

%T 18381446004,323527186957,6114296752718,123513004310640,

%U 2655648779976640,60554669008300565,1459559515622280282,37079264125376670955,990226180225789628660,27733277682719819190246,812818183963966524137332,24880254143735238825011057

%N Expansion of the e.g.f. sqrt(1 / (2*exp(x) - 2*x*exp(x) - 1)).

%F E.g.f. y(x) satisfies y' = x*exp(x)*y^3.

%F a(0)=1, a(n) = Sum_{k=1..floor(n/2)} A014307(k)*A008306(n,k) for n >= 1.

%F For all p prime, a(p) == -1 (mod p).

%F a(n) ~ sqrt(2*c) * n^n / ((1-c)^(n+1) * exp(n)), where c = -LambertW(-exp(-1)/2). - _Vaclav Kotesovec_, Jun 25 2021

%e sqrt(1/(2*exp(x)-2*x*exp(x)-1)) = 1 + x^2/2! + 2*x^3/3! + 12*x^4/4! + 64*x^5/5! + 485*x^6/6! + 4038*x^7/7! + 39991*x^8/8! + 441992*x^9/9! + ...

%e a(13) = Sum_{k=1..6} A014307(k)*A008306(13,k) = 18381446004.

%e A014307(1)*A008306(13,1) == -1 (mod 13), because A014307(1) = 1 and A008306(13,1) = (13-1)!

%e For k>=2, A008306(13,k) == 0 (mod 13), result a(13) == -1 (mod 13).

%p A014307 := proc(n) option remember; `if`(n=0, 1 , 1+add((-1+binomial(n, k))*A014307(k), k=1..n-1)) end:

%p A008306 := proc(n, k): if k=1 then (n-1)! ; elif n<=2*k-1 then 0; else (n-1)*procname(n-1, k)+(n-1)*procname(n-2, k-1) ; end if; end proc:

%p a := n-> add((A014307(k)*A008306(n,k)), k=1..floor(n/2)):a(0):=1 ;

%p seq(a(n), n=0..24);

%p # second program:

%p a := series(sqrt((1/(2*exp(x)-2*x*exp(x)-1))), x=0, 25):

%p seq(n!*coeff(a, x, n), n=0..24);

%t CoefficientList[Series[Sqrt[1/(2*E^x-2*x*E^x-1)], {x, 0, 24}], x] * Range[0, 24]!

%o (PARI) my(x='x+O('x^25)); Vec(serlaplace(sqrt(1 / (2*exp(x) - 2*x*exp(x) -1)))) \\ _Michel Marcus_, Jun 24 2021

%Y Cf. A000166, A008306, A014307, A327006.

%K nonn

%O 0,4

%A _Mélika Tebni_, Jun 24 2021