Issues related to factorial2 function

[yingxingcheng]

Issues related to factorial2 function

[PaulWAyers]

Yes, the documentation is embarrassingly bad.

It took me a bit of time to remember how the integral works. It would have helped if we would have documented it. The integral is

$$ \int_{-\infty}^{\infty} (x-x_1)^{n_1} (x-x_2)^{n_2} e^{-\alpha_1 (x-x_1)^2} e^{-\alpha_2 (x-x_2)^2} dx $$

but before this funtion is invoked the Gaussian product rule has been invoked, and the utility variable

$$ \mathtt{two}\textunderscore\mathtt{at} = 2\frac{\alpha_1 x_1 + \alpha_2 x_2}{\alpha_1 + \alpha_2} $$

was defined. Then the constant prefactor

$$ \sqrt{\frac{2\pi}{\mathtt{two}\textunderscore\mathtt{at}}} \exp\left(-\frac{\alpha_1 \alpha_2}{\alpha_1 + \alpha_2}(x_1-x_2)^2 \right) $$

is factored out. (At least these utility variables are relatively clearly). So this function actually computes

$$ \sqrt{\frac{\mathtt{two}\textunderscore\mathtt{at}}{2\pi}}\int_{-\infty}^{\infty} (x-x_1)^{n_1} (x-x_2)^{n_2} e^{-\tfrac{1}{2}\mathtt{two}\textunderscore\mathtt{at} \cdot x^2} dx $$

Then, using the standard Gaussian integral formula,

$$ \int_{-\infty}^{\infty} x^{i+j} e^{-\tfrac{1}{2}\mathtt{two}\textunderscore\mathtt{at} \cdot x^2} dx = \begin{cases} \frac{(i+j-1)!!}{\mathtt{two}\textunderscore\mathtt{at}^{i+j}}\sqrt{\frac{2\pi}{\mathtt{two}\textunderscore\mathtt{at}}}, & (i+j)\mod{2} = 0\\ 0 & (i+j)\mod{2} = 1 \end{cases} $$

We have:

$$ \sqrt{\frac{\mathtt{two}\textunderscore\mathtt{at}}{2\pi}}\int_{-\infty}^{\infty} (x-x_1)^{n_1} (x-x_2)^{n_2} e^{-\tfrac{1}{2}\mathtt{two}\textunderscore\mathtt{at} \cdot x^2} dx = \sum_{i=0}^{n_1} \sum_{{j=0} \atop {(i+j) \mod 2 = 0}}^{n_2} \binom{n_1}{i} \binom{n_2}{j} x_1^{n_1-i} x_2^{n_2 - j} \frac{(i+j-1)!!}{\mathtt{two}\textunderscore\mathtt{at}^{i+j}} $$

Realistically, this is overcomplicated: IMO some of the things we did to "save evaluations in the inner loop" were probably not worth the loss in readability, and at least demanded more documentation. It also seems to me that the double-factorial we are using is off by one, but I can't imagine how we managed to pass the old tests if that is the case, so I must be missing something.

[tovrstra]

Sorry, forgot to mention that this is fixed in #319