Issues related to factorial2 function

[yingxingcheng]

The issue with the factorial2 function persists, as the GitHub test workflow fails.

Additionally, the second argument of factorial2 is of bool type. The following usage of factorial2 is quite confusing:

iodata/iodata/overlap.py

Lines 225 to 243 in 6b9c6e0

	class GaussianOverlap:
	"""Gaussian Overlap Class."""
	def __init__(self, n_max):
	"""Initialize class.
	Parameters
	----------
	n_max : int
	Maximum angular momentum.
	"""
	self.binomials = [
	[scipy.special.binom(n, i) for i in range(n + 1)] for n in range(n_max + 1)
	]
	facts = [factorial2(m, 2) for m in range(2 * n_max)]
	facts.insert(0, 1)
	self.facts = np.array(facts)

[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