%I M3985 #75 Feb 27 2022 15:48:27
%S 1,5,37,353,4081,55205,854197,14876033,288018721,6138913925,
%T 142882295557,3606682364513,98158402127761,2865624738913445,
%U 89338394736560917,2962542872271918593,104128401379446177601
%N a(n) = n * (2*n - 1)!! - Sum_{k=0..n-1} a(k) * (2*n - 2*k - 1)!!.
%C a(n+1) is the moment of order n for the probability density function rho(x) = Pi^(-3/2)*sqrt(x/2)*exp(x/2)/(1-erf^2(i*sqrt(x/2))) on the interval 0..infinity, where erf is the error function and i=sqrt(-1). - _Groux Roland_, Nov 10 2009
%D E. W. Bowen, Letter to N. J. A. Sloane, Aug 27 1976.
%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).
%H Vincenzo Librandi, <a href="/A004208/b004208.txt">Table of n, a(n) for n = 1..300</a>
%H E. W. Bowen & N. J. A. Sloane, <a href="/A003319/a003319.pdf">Correspondence, 1976</a>
%H Colin K, <a href="http://mathematica.stackexchange.com/questions/6036/">Is it possible to make Mathematica reformulate an expression in a more numerically stable way?</a>. see second answer by whuber.
%F a(n) = (1/2) * A000698(n+1), n > 0.
%F x + (5/2)*x^2 + (37/3)*x^3 + (353/4)*x^4 + (4081/5)*x^5 + (55205/6)*x^6 + ... = log(1 + x + 3*x^2 + 15*x^3 + 105*x^4 + 945*x^5 + 10395*x^6 + ...) where [1, 1, 3, 15, 105, 945, 10395, ...] = A001147(double factorials). - _Philippe Deléham_, Jun 20 2006
%F G.f.: ( 1/Q(0) - 1)/x where Q(k) = 1 - x*(2*k+1)/(1 - x*(2*k+4)/Q(k+1)); (continued fraction). - _Sergei N. Gladkovskii_, Mar 19 2013
%F G.f.: (2/x)/G(0) - 1/x, where G(k) = 1 + 1/(1 - 2*x*(2*k+1)/(2*x*(2*k+1) - 1 + 2*x*(2*k+4)/G(k+1))); (continued fraction). - _Sergei N. Gladkovskii_, May 31 2013
%F G.f.: 1/(2*x^2) - 1/(2*x) - G(0)/(2*x^2), where G(k) = 1 - x*(k+1)/G(k+1); (continued fraction). - _Sergei N. Gladkovskii_, Aug 15 2013
%F L.g.f.: log(1/(1 - x/(1 - 2*x/(1 - 3*x/(1 - 4*x/(1 - 5*x/(1 - ...))))))) = Sum_{n>=1} a(n)*x^n/n. - _Ilya Gutkovskiy_, May 10 2017
%p df := proc(n) product(2*k-1,k=1..n) end: a[1] := 1: for n from 2 to 30 do a[n] := n*df(n)-sum(a[k]*df(n-k),k=1..n-1) od;
%t CoefficientList[Series[D[Log[Sum[(2n-1)!!x^n,{n,0,17}]],x],{x,0,16}],x] [From _Wouter Meeussen_, Mar 21 2009]
%t a[ n_] := If[ n < 1, 0, n Coefficient[ Normal[ Series[ Log @ Erfc @ Sqrt @ x, {x, Infinity, n}] + x + Log[ Sqrt [Pi x]]] /. x -> -1 / 2 / x, x, n]] (* _Michael Somos_, May 28 2012 *)
%o (PARI) {a(n) = if( n<1, 0, n++; polcoeff( 1 - 1 / (2 * sum( k=0, n, x^k * (2*k)! / (2^k * k!), x * O(x^n))), n))} /* _Michael Somos_, May 28 2012 */
%Y Cf. A000698.
%K nonn
%O 1,2
%A _N. J. A. Sloane_, following a suggestion from E. W. Bowen, Aug 27 1976
%E Description corrected by Jeremy Magland (magland(AT)math.byu.edu), Jan 07 2000
%E More terms from _Emeric Deutsch_, Dec 21 2003