The mathematics underlying this baroque expression of the double-factorial

On the OEIS page for the double factorial, there are three ways of getting the sequence in PARI. One of them is this curious bit of PARI code:

a(n)=local(E); E=exp(x^2/2+x*O(x^n)); n!*polcoeff(1+E*x*(1+intformal(1/E)), n)


Compared to the others, the above may seem extraordinarily baroque:

a(n) = prod(i=0, floor((n-1)/2), n-2*i );
a(n)=if(n<2, 1, n*a(n-2))


In the example I'm asking about, exp stands for the exponential function, O for the zeroes of a p-adic or power series, intformal for formal integration.

Is this a well-known use of the exponential function?

Is it known who first came up with this?

All of that seems like an amazing amount of mathematics to calculate the double factorial, and I'd like to know more about it.

-
It's a generating function, a very powerful tool for understanding sequences. There is more on MathWorld mathworld.wolfram.com/DoubleFactorial.html and in this article jstor.org/stable/10.4169/math.mag.85.3.177 ("Double fun with double factorials") – user31373 Jun 7 '12 at 17:17

To settle this:

As Leonid says, that PARI code you saw is in fact a way to generate the double factorials from an appropriate (exponential) generating function.

$$1+x \exp\left(\frac{x^2}{2}\right)\left(1+\int_0^x \exp\left(-\frac{t^2}{2}\right) \mathrm dt\right)=\sum_{k=0}^\infty \frac{k!!}{k!}x^k$$
or, in terms of the error function $\mathrm{erf}(x)$,
$$1+x \exp\left(\frac{x^2}{2}\right)\left(1+\sqrt{\frac{\pi}{2}}\mathrm{erf}\left(\frac{x}{\sqrt 2}\right)\right)=\sum_{k=0}^\infty \frac{k!!}{k!}x^k$$
(In Mathematica, you can show this with the snippet FullSimplify[ExponentialGeneratingFunction[k!!, k, x]].)
Thus, to obtain the double factorials, one differentiates the exponential generating function an appropriate number of times (which is precisely what you must do to obtain a series coefficient), evaluate at $x=0$, and then multiply by the corresponding factorial.