Double Factorial

Most mathematically-inclined folks know what a factorial is. The simplest, recursive definition is given by:

f(0) = 1
f(n) = n * f(n-1) = n!

Figure 1: Definition of Factorial

For all n in the set of Natural Numbers, ℕ.

In the set of real numbers, ℝ, this can be extended to the Gamma function, Γ, which has the form:

\Gamma(x) = \int_0^\infty t^{x-1} \mathrm{e}^{-t} \mathrm{d}t
\forall n \in \mathbb{N}: \Gamma(n) = (n-1)!

Figure 2: Definition of Gamma

Both Figure 1 and Figure 2 are typically well known definitions. But what’s less known is the double, triple, and even quadruple factorial:

n!! = n (n-2)!!
n!!! = n (n-3)!!!
n!!!! = n (n-4)!!!!
\vdots

Figure 3: Higher Level Factorial Functions

The question is, how do you extend this from Natural Numbers to the general case with Real Numbers, like with the Gamma function from Figure 2?

Fortunately, Python’s math library has math.gamma(n) to compute the Gamma function on a given number. Thus, rather than the recursive definition in Figure 1, we can directly define factorial in terms of Gamma:

def factorial(x):
    return math.gamma(x+1)

Figure 4: Definition of factorial(x)

The question is, can the Double and Higher Factorials be defined in terms of gamma?

The Pochhammer Function

The Pochhammer function is defined in therms of Gamma. Specifically:

(x)_n = \frac{\Gamma (x+n)}{\Gamma (x)} = x(x+1) \cdots (x+n-1)

Figure 5: The Pochhammer Function

The corresponding Python looks like this:

def Pochhammer(n, k):
    return math.gamma(k+m)/math.gamma(n)

Figure 6: Definition of Pochhammer

The Pochhammer brings us closer to our ideal continuous multi-factorial method. For instance, we could assume a version of Factorial over the Rational Numbers, ℚ:

def Rational_Multifactorial(n, numerator, denominator):
    return pow(numerator/denominator, n*denominator/numerator) * Pochhammer(1, n*denominator/numerator)

Figure 7: Definition of Rational_Multifactorial

Unfortunately, that doesn’t give us a Real solution, ℝ. Instead, if we’re to ask WolframAlpha to work it out for us, we end up with something completely different:

x!! = 2^{\frac{x}{2}} (c_2(-1)^x+c_1) \Gamma (1+\frac{x}{2})
where,
c_1 = \frac{1}{2} + \frac{1}{\sqrt{2 \pi}}
c_2 = \frac{1}{2} - \frac{1}{\sqrt{2 \pi}}

Figure 8: WolframAlpha definition of Double Factorial over the Real

The way that was calculated was the simple recursion in Figure 3 was entered and WolframAlpha was told to find an arithmetic solution. The answer was expressed in terms of the constants c1 and c2, which I then solved for knowing a few fixed values for the multiplication.

The result in Python, becomes:

def double_factorial(x):
    c1 = 0.5 + 1/math.sqrt(2*math.pi)
    c2 = 0.5 - 1/math.sqrt(2*math.pi)

    return pow(2, x/2) * (c2*pow(-1, x)+c1)*math.gamma(x/2+1)

Figure 9: Definition of double_factorial

What’s more, the rational solution used Pochhammer, but the real solution used only Gamma. I was able to use WolframAlpha to compute analytic functions for Triple and Quadruple factorials, with their three and four constants, respectively, but there didn’t seem to be any pattern to the calculations and the solutions move on to using Sine and Cosine. And none of them used Pochhammer. Yet the Rational solution in Figure 7 should work for all levels of Factorial.

I’ve writing up all my functions in Factorial.py, which you can find on my Subversion repository and I’m currently available for hire.