How to put a series in hypergeometric form

Posted on 18 November 2021 by John

I skipped a step in the previous post, not showing how a series was put into the form of a hypergeometric function.

Actually I skipped two steps. I first said that a series was not obviously hypergeometric, and yet at first glance it sorta is.

I’d like to make up for both of these omissions, but with a simpler example. It looks like a hypergeometric series, and it is, but it’s not the one you might think. The example will be the function below.

$f(x) = \sum_{n=0}^\infty \binom{5n}{n} \frac{x^n}{n!}$

Rising powers

Hypergeometric functions are defined by series whose coefficients are rising powers.

The nth rising power of a is

(a)_n = a (a+1) (a + 2) … (a + n – 1).

More on rising powers here. The coefficient of xⁿ in a series defining a hypergeometric function is a ratio of nth rising powers of constants.

Why the example isn’t obvious

You can write factorials as rising powers: n! = (1)_n. The binomial coefficient in our example is a ratio of rising powers:

(5n)! = (1)_5n

in the numerator and

n! (4n)! = (1)_n (1)_4n

in the denominator. So why aren’t we done? Because although these are all rising powers, they’re not all nth rising powers. The subscripts on the (a)_n are all different: 5n, 4n, and n.

So is our function not hypergeometric? It is, but we’ll have to introduce more parameters.

Rational polynomials

Another way to define hypergeometric series is to say that the ratio of consecutive coefficients is a rational function of the index n. This is very succinct, but not as explicit as the definition in terms of rising powers, though they’re equivalent.

In addition to brevity, this definition has another advantage: the hypergeometric parameters are the negatives of the zeros and poles of said rational function. The zeros are the upper parameters and the poles are the lower parameters.

This is how you put a series in hypergeometric form, if it has such a form. It’s also how you test whether a series is hypergeometric: a series is hypergeometric if and only if the ratio of consecutive terms is a single rational function of the index.

Back to our example

The ratio of the (n+1)st term to the nth term in our series is

$\begin{align*} \frac{t_{n+1}}{t_n} &= \left.\binom{5(n+1)}{n+1} \frac{z^{n+1}}{(n+1)!} \middle/{\binom{5n}{n}} \frac{z^n}{n!} \right. \\ &= \frac{(5n+5)!}{(n+1)! \,(4n+4)!} \,\frac{n!\,(4n)!}{(5n)!} \, \frac{n!}{(n+1)!}\frac{5^5}{4^4} z\\ &= \frac{(5n+5)(5n+4)(5n+3)(5n+2)(5n+1)}{(4n+4)(4n+3)(4n+2)(4n+1)(n+1)^2} \, \frac{5^5}{4^4} z \end{align*}$

and the final expression is a rational function of n. We can read the hypergeometric function parameters directly from this rational function. The upper parameters are 1, 4/5, 3/5, 2/5, and 1/5. The lower parameters are 1, 3/4, 2/4, 1/4, and 1. Identical factors in the upper and lower parameters cancel each other out, so we can remove the 1 from the list of upper parameters and remove one of the 1’s from the list of lower parameters [1].

So our example function is

${}_4 F_4\left(\frac{1}{5}, \frac{2}{5}, \frac{3}{5}, \frac{4}{5};\,\, \frac{1}{4}, \frac{1}{2}, \frac{3}{4}, 1; \,\,\frac{3125}{256}x \right)$

The order of the upper parameters doesn’t matter, and neither does the order of the lower parameters, though it matters very much which ones are upstairs and which are downstairs. The subscript to the left of the F gives the number of upper parameters and the subscript to the right gives the number of lower parameters. These subscripts are redundant since you could just count the number of parameters between the semicolons, but they make it easier to tell at a glance what class of hypergeometric function you’re working with.

Spot check

Let’s evaluate our original series and our hypergeometric series at a point to make sure they agree. I chose x = 0.01 to stay within the radius of convergence of the hypergeometric series.

Both

    Sum[Binomial[5 n, n] (0.01^n)/n!, {n, 0, Infinity}]

and

    HypergeometricPFQ[
        {1/5, 2/5, 3/5, 4/5}, 
        {1/4, 2/4, 3/4, 1},
        0.01*3125/256]

return the same value, namely 1.05233.

[1] The extra factor of n! in the denominator of hypergeometric functions complicates things. It seems like an unfortunately historical artifact to me, but maybe it’s a net benefit for reasons I’m not thinking of. When we look at zeros of the numerator and denominator, we have a single 1 on top and three ones on bottom. The 1 on top cancels out one of the 1s on bottom, but why don’t we have two 1s in the lower parameters? Because one of them corresponds to the extra n! in the denominator.

Quintic trinomial root

Posted on 18 November 2021 by John

This post looks at an exercise from Special Functions, exercise 6 in Appendix E.

Suppose that x^m+1 + ax – b = 0. Show that

$x = \frac{b}{a} - \frac{b^{m+1}}{a^{m+2}} + \frac{2m + 2}{2!} \frac{b^{2m+1}}{a^{2m+3}} - \frac{(3m+2)(3m+3)}{3!} \frac{b^{3m+1}}{a^{3m+4}} + \cdots$

Use this formula to find a solution to x⁵ + 4x + 2 = 0 to four decimal places of accuracy. When m = 0 this series reduces to the geometric series. Write this sum as a hypergeometric series.

There are several things I find interesting about this exercise.

It’s a possibly useful result.
It’s an application of Lagrange’s inversion formula; that’s where the series comes from.
The equation has m+1 roots. What’s special about the root given by the series?
The sum is not obviously hypergeometric.
What happens when the sum diverges? Is it a useful asymptotic series? Does the analytic continuation work?

I want to address the last two points in this post. Maybe I’ll go back and address the others in the future. To simplify things a little, I will assume m = 4.

The parameters of the hypergeometric series are not apparent from the expression above, nor is it even apparent how many parameters you need. It turns out you need m upper parameters and m-1 lower parameters. Since we’re interested in m = 4, that means we have a hypergeometric function of the form ₄F₃.

$x = \frac{b}{a}\,\, {}_4 F_3\left(\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5};\frac{1}{2}, \frac{3}{4}, \frac{5}{4} ; -\frac{3125b^4}{256a^5}\right)$

We can evaluate this expression in Mathematica as

    f[a_, b_] := (b/a) HypergeometricPFQ[
        {1/5, 2/5, 3/5, 4/5}, 
        {1/2, 3/4, 5/4}, 
        -3125 b^4 / (256 a^5)
    ]

When we evaluate f[4, -2] we get -0.492739, which is the only real root of

x⁵ + 4x + 2 = 0.

Recall that the sign on b is negative, so we call our function with b = -2.

Now lets, try another example, solving for a root of

x⁵ + 3x – 4 = 0.

If we plug a = 3 and b = 4 into the series at the top of the post, the series diverges immediately, not giving a useful asymptotic series before it diverges.

The series defining our hypergeometric function diverges for |z| > 1, and we’re evaluating it at z = -3125/243. However, the function can be extended beyond the unit disk by analytic continuation, and when we ask Mathematica to numerically evaluate the root with by running

    N{ f[3, 4] ]

we get 1, which is clearly a root of x⁵ + 3x – 4.

Numbering the branches of the Lambert W function

Posted on 15 November 2021 by John

The previous post used the Lambert W function to solve an equation that came up in counting partitions of an integer. The first solution found didn’t make sense in context, but another solution, one on a different branch, did. The default branch k = 0 wasn’t what we were after, but the branch k = -1 was.

Branches 0 and -1

The branch corresponding to k = 0 is the principal branch. In mathematical notation, the k is usually a subscript, i.e W₀(z) is the principal branch of the Lambert W function. For real z ≥ -1/e it returns a real value. It’s often what you want, and that’s why it’s the default in software implementations such as SciPy and Mathematica. More on that below.

The only other branch that returns real values for real inputs is W_-1 which returns real values for -1/e ≤ z < 0.

If you’re working strictly with real numbers, you only need to be concerned with branches 0 and -1. If branch 0 doesn’t give you what you expect, try branch -1, if your argument is negative.

SciPy and Mathematica

SciPy implements W_k with lambertw(z). The parameter k defaults to 0.

The Mathematica function ProductLog[z] implements W₀(z), and ProductLog[k, z] implements W_k.

Note that Mathematica and Python list their arguments in opposite orders.

Branch numbering

The branches are numbered in order of their imaginary parts. That is, The imaginary part of W_k (z) is an increasing function of k. For example, here is a list of the numerical values of the imaginary parts of W_k (1) for k running from -3 to 3.

    Table[N[Im[ProductLog[k, 1]]], {k, -3, 3}]
    {-17.1135, -10.7763, -4.37519, 0., 4.37519, 10.7763, 17.1135}

Note the zero in the middle because W₀(1) is real.

Recovering k

Given z and a value W_k (z) with k unknown, you can determine k with

$k = \frac{W_k(z) + \log W_k(z) - \log z}{2\pi i}$

with the exception that if the expression above is 0 and -1/e ≥ z < 0 then k = -1. See [1].

For example, let z = 1 + 2i and w = W₃ (z).

    z = 1 + 2 I
    w = ProductLog[3, z]

Now pretend you don’t know how w was computed. When you compute

    N[(w + Log[w] - Log[z])/(2 Pi I)]

the result is 3.

Partition example

The discussion of the W function here grew out of a problem with partitions. Specifically, we wanted to solve for n such that
$f(n) = \frac{1}{4n\sqrt{3}} \exp\left( \pi \sqrt{\frac{2n}{3}} \right)$

is approximately a trillion. We found in the previous post that solutions of the equation

$a w^b \exp(c w^d) = x$

are given by

$w = \left(\frac{b}{cd} \,\,W\left(\frac{cd}{b} \left(\frac{x}{a} \right )^{d/b} \right )\right )^{1/d}$

Our partition problem corresponds to a = 1/√48, b = -1, c = π √(2/3), and d = 1/2. We found that the solution we were after came from the k = -1 branch.

    N[( (b/(c d)) ProductLog[-1, (x/a)^(d/b) c d /b])^(1/d)]
    183.852

Even though our final value was real, the intermediate values were not. The invocation of W_-1 in the middle returns a complex value:

    N[ProductLog[-1, ((x/a)^(d/b) c d /b)^(1/d)]]
    -32.5568 - 3.24081 I

We weren’t careful about specifying ranges and branch cuts, but just sorta bulldozed our way to a solution. But it worked: it’s easy to verify that 183.852 is the solution we were after.

***

[1] Unwinding the branches of the Lambert W function by Jeffrey, Hare, and Corless. The Mathematical Scientist 21, 1–7 (1996)

Solving equations with Lambert’s W function

Posted on 14 November 2021 by John

In the previous post we wanted to find a value of n such that f(n) = 10¹² where

$f(n) = \frac{1}{4n\sqrt{3}} \exp\left( \pi \sqrt{\frac{2n}{3}} \right)$

and we took a rough guess n = 200. Turns out f(200) ≈ 4 × 10¹² and that was good enough for our purposes. But what if we wanted to solve f(n) = x accurately?

We will work our way up to solving f(n) = 10¹² exactly using the Lambert W function.

Lambert’s W function

If you can solve one equation involving exponential functions you can bootstrap your solution solve a lot more equations.

The Lambert W function is defined to be the function W(x) that maps each x to a solution of the equation

w exp(w) = x.

This function is implemented Python under scipy.special.lambertw. Let’s see how we could use it solve similar equations to the one that define it.

Basic bootstrapping

Given a constant c we can solve

cw exp(w) = x

simply by solving

w exp(w) = x/c,

which says w = W(x/c).

And we can solve

w exp(cw) = x

by setting y = cw and solving

y exp(y) = cx.

and so cw = W(cx) and w = W(cx)/c.

Combining the two approaches above shows that the solution to

aw exp(bw) = x

w = W(bx/a) / b.

We can solve

exp(w) / w = x

by taking the reciprocal of both sides and applying the result above with a = 1 and b = -1.

More generally, we can solve

w^c exp(w) = x

by raising both sides to the power 1/c; the reciprocal the special case c = 1.

Similarly, we can solve

w exp(w^c) = x

by setting y = w^c , changing the problem to

y^-1/c exp(y) = x.

Putting it all together

We’ve now found how to deal with constants and exponents on w, inside and outside the exponential function. We now have all the elements to solve our original problem.

We can solve the general equation

$a w^b \exp(c w^d) = x$

with

$w = \left(\frac{b}{cd} \,\,W\left(\frac{cd}{b} \left(\frac{x}{a} \right )^{d/b} \right )\right )^{1/d}$

The equation at the top of the post corresponds to a = 1/√48, b = -1, c = π √(2/3), and d = 1/2.

We can code this up in Python as follows.

    from scipy.special import lambertw

    def soln(a, b, c, d, x):
        t = (x/a)**(d/b) *c*d/b
        return (lambertw(t)*b/(c*d))**(1/d)

We can verify that our solution really is a solution by running it though

    from numpy import exp, pi

    def f(a, b, c, d, w):
        return a*w**b * exp(c*w**d)

to make sure we get our x back.

When we run

    a = 0.25*3**-0.5
    b = -1
    c = pi*(2/3)**0.5
    d = 0.5
    x = 10**12

    w = soln(a, b, c, d, x)
    print(f(a, b, c, d, w))

we do indeed get x back.

    a = 0.25*3**-0.5
    b = -1
    c = pi*(2/3)**0.5
    d = 0.5

    w = soln(a, b, c, d, 10)
    print(f(a, b, c, d, w))

When we take a look at the solution w, we get 1.443377079584187e-13. In other words, we get a number near zero. But our initial guess was w = 200. What went wrong?

Nothing went wrong. We got a solution to our equation. It just wasn’t the solution we expected.

The Lambert W function has multiple branches when viewed as a function on complex numbers. It turns out the solution we were expecting is on the branch SciPy denotes with k = -1. If we add this to the call to lambertw above we get the solution 183.85249773880142 which is more in line with what we expected.

Mathematical plot twist: q-binomial coefficients

Posted on 12 November 2021 by John

The simplest instance of a q-analog in math is to replace a positive integer n by a polynomial in q that converges to n as q → 1. Specifically, we associate with n the polynomial

[n]_q = 1 + q + q² + q³ + … + q^n-1.

which obviously converges to n when q → 1.

More generally we can associate with a real number a the rational polynomial

[a]_q = (1 – q^a) / (1 – q).

When a is a positive integer, the two definitions are the same. When a is not an integer, you can see by L’Hôpital’s rule that

lim_{q → 1} [a]_q = a.

Motivation

Why would you do this? I believe the motivation was to be able to define a limit of a series where the most direct approach fails or is difficult to work with. Move over to the land of q-analogs, do all your manipulations there where requiring |q| < 1 avoids problems with convergence, then take the limit as q goes to 1 when you’re done.

q-factorials and q-binomial coefficients

The factorial of a positive integer is

n! = 1 × 2 × 3 × … × n

The q-factorial of n replaces each term in the definition of factorial with its q-analog:

[n]_q! = [1]_q [2]_q [3]_q … [n]_q

(It helps to think of []_q! as a single operator defining the q-factorial as above. This is not the factorial of [n]_q.)

Similarly, you can start with the definition of a binomial coefficient

n! / (k! (n–k)!)

and define the q-binomial coefficient by replacing factorials with q-factorials.

[n]_q! / ([k]_q! [n–k]_q!)

You can apply the same process to rising and falling powers, and then replace the rising powers in the definition of a hypergeometric function to define q-analogs of these functions.

(The q-analogs of hypergeometric functions are confusingly called “basic hypergeometric functions.” They are not basic in the sense of simple examples; they’re not examples at all but new functions. And they’re not basic in the sense of a basis in linear algebra. Instead, the name comes from q being the base of an exponent.)

Plot twist

And now for the plot twist promised in the title.

I said above that the q-analogs were motivated by subtle problems with convergence. There’s the idea in the background that eventually someone is going to let q approach 1 in a limit. But what if they don’t?

The theorems proved about q-analogs were first developed with the implicit assumption that |q| < 1, but at some point someone realized that the theorems don’t actually require this, and are generally true for any q. (Maybe you have to avoid q = 1 to not divide by 0 or something like that.)

Someone observed that q-analog theorems are useful when q is a prime power. I wonder who first thought of this and why. It could have been inspired by the notation, since q often denotes a prime power in the context of number theory. It could have been as simple as a sort of mathematical pun. What if I take the symbol q in this context and interpret it the way q is interpreted in another context?

Now suppose we have a finite field F with q elements; a theorem says q has to be a prime power. Form an n-dimensional vector space over F. How many subspaces of dimension k are there in this vector space? The answer is the q-binomial coefficient of n and k. I explore this further in the next post.

Three error function series

Posted on 7 November 2021 by John

A common technique for numerically evaluating functions is to use a power series for small arguments and an asymptotic series for large arguments. This might be refined by using a third approach, such as rational approximation, in the middle.

The error function erf(x) has alternating series on both ends: its power series and asymptotic series both are alternating series, and so both are examples that go along with the previous post.

The error function is also a hypergeometric function, so it’s also an example of other things I’ve written about recently, though with a small twist.

The error function is a minor variation on the normal distribution CDF Φ. Analysts prefer to work with erf(x) while statisticians prefer to work with Φ(x). See these notes for translating between the two.

Power series

The power series for the error function is

$\text{erf}(x) = \frac{2}{\sqrt{\pi}} \sum_{n=0}^\infty (-1)^n \frac{x^{2n+1}}{n!(2n+1)}$

The series converges for all x, but it is impractical for large x. For every x the alternating series theorem eventually applies, but for large x you may have to go far out in the series before this happens.

Asymptotic series

The complementary error function erfc(x) is given by

erfc(x) = 1 – erf(x).

Even though the relation between the two functions is trivial, it’s theoretically and numerically useful to have names for both functions. See this post on math functions that don’t seem necessary for more explanation.

The function erfc(x) has an alternating asymptotic series.

$\text{erfc}(x) = \frac{e^{-x^2}}{\sqrt{\pi}x} \sum_{n=0}^\infty (-1)^n\frac{(2n-1)!!}{2^n} x^{-2n}$

Here a!!, a double factorial, is the product of the positive integers less than or equal to a and of the same parity as a. Since in our case a = 2n – 1, a is odd and so we have the product of the positive odd integers up to and including 2n – 1. The first term in the series involves (-1)!! which is defined to be 1: there are no positive integers less than or equal to -1, so the double factorial of -1 is an empty product, and empty products are defined to be 1.

Confluent hypergeometric series

The error function is a hypergeometric function, but of a special kind. (As is often the case, it’s not literally a hypergeometric function but trivially related to a hypergeometric function.)

The term “hypergeometric function” is overloaded to mean two different things. The strictest use of the term is for functions F(a, b; c; x) where the parameters a and b specify terms in the numerator of coefficients and the parameter c specifies terms in the denominator. You can find full details here.

The more general use of the term hypergeometric function refers to functions that can have any number of numerator and denominator terms. Again, full details are here.

The special case of one numerator parameter and one denominator parameter is known as a confluent hypergeometric function. Why “confluent”? The name comes from the use of these functions in solving differential equations. Confluent hypergeometric functions correspond to a case in which two singularities of a differential equation merge, like the confluence of two rivers. More on hypergeometric functions and differential equations here.

Without further ado, here are two ways to relate the error function to confluent hypergeometric functions.

$\text{erf}(x) = \frac{2x}{\sqrt{\pi}} F\left(\frac{1}{2}; \frac{3}{2}; -x^2\right) = \frac{2x}{\sqrt{\pi}} e^{-x^2} F\left(1; \frac{3}{2}; x^2\right)$

The middle expression is the same as power series above. The third expression is new.

Escaping the unit disk

Posted on 3 November 2021 by John

Hypergeometric functions are defined in terms of their power series representation. This power series converges inside the unit circle and diverges outside. The functions extend to most [1] of the complex plane, and do so in a computationally nice way.

Analytic continuation is often sort of existential: we prove that it can be done, but we may not be able to say much explicitly about the result. But for hypergeometric functions, there is an explicit formula.

This formula is so complicated that it’s possible to miss its basic structure. Here’s an outline version:

$F(\ldots z) = \ldots F(\ldots z^{-1}) + \ldots F(\ldots z^{-1})$

This says that F evaluated at a point z outside of the unit circle is equal to the sum of two terms involving F evaluated at 1/z, a point inside the unit circle. This means you can evaluate F everywhere [1] if you can evaluate it inside the unit circle.

The parameters in each of the F‘s above are different, so the formula doesn’t relate a single function’s value at z to the same function’s values at 1/z. Instead, it relates the value of one instance of a class of functions at z to values of other instances of the class at 1/z.

And now without further ado, here is the analytic continuation formula in all its glory:

$\begin{align*} \frac{\Gamma(a)\Gamma(b)}{\Gamma(c)}F(a, b;c; z) &= \frac{\Gamma(a)\Gamma(b-a)}{\Gamma(c-a)} \,\,(-z)^{-a}\,\,F(a, 1-c+a; 1-b+a; z^{-1}) \\ &+ \,\frac{\Gamma(b)\Gamma(a-b)}{\Gamma(c-b)} \,\,(-z)^{-b}\,\,F(b,\, 1-c+b; 1-a+b; z^{-1}) \end{align*}$

[1] There may be a singularity at 1, and there may be a branch cut along the real axis from 1 to infinity.

Connections

Posted on 31 October 2021 by John

One of my favorite things in math is connecting together things that do not seem related at first. This post will elaborate on a connection between two recent posts.

Gauss’s constant

A couple weeks ago I wrote about Gauss’s constant

$g \equiv \frac{1}{\text{AGM}(1, \sqrt{2})} = 0.834626\ldots$

and several things it was related to. We can calculate g with n iterations of the process defining the AGM (arithmetic-geometric mean) as follows.

    def f(n):
        a, b = 1, 2**0.5
        for _ in range(n):
            a, b = 0.5*(a+b), (a*b)**0.5
        return 1/a

This process converges very quickly, calculating g to full machine precision when n = 4.

Elliptic integrals

Gauss’s constant, and the AGM more generally, can be computed in terms of an elliptic integral, and that in turn can be computed in terms of a hypergeometric function.

Elliptic integrals are so named because historically they are related to the problem of finding the perimeter of an ellipse. The particular elliptic integral we’re interested here is Legendre’s K function

$K(k) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - k^2 \sin^2\theta}}$

(There are several conventions for parameterizing elliptic integrals, so this definition may not match what you see elsewhere. See last section.)

With this definition of K,

$\frac{1}{\text{AGM}(x,y)} = \frac{2}{\pi y} K(\sqrt{1 - y^2/x^2})$

and by the symmetry of the AGM, you can reverse the order of x and y. It follows that

$g = \frac{\sqrt{2}}{\pi}K\left(\frac{1}{\sqrt{2}}\right) = \frac{2}{\pi}K(i)$

Hypergeometric functions

The previous post briefly introduced hypergeometric functions. The function K is related to Gauss’s hypergeometric function by

$K(k) = \frac{\pi}{2} F\left(\frac{1}{2}, \frac{1}{2}; 1; k^2 \right)$

Expanding F as a series shows the connection to double factorials.

$\begin{align*} F\left(\frac{1}{2}, \frac{1}{2}; 1; x \right) &= \sum_{n=0}^\infty \frac{(1/2)_n (1/2)_n}{(1)_n n!} x^n \\ &= \sum_{n=0}^\infty \left( \frac{(2n-1)!!}{(2n)!!} \right)^2 x^n \end{align*}$

This means you could numerically evaluate a power series to compute Gauss’s constant, thought that would be inefficient. Because the AGM converges so quickly, you’d want to evaluate elliptic integrals in terms of the AGM, not evaluate the AGM in terms of elliptic integrals.

Conventions

There is some confusing terminology and variation in conventions when defining elliptic integrals. We’ve defined K here in terms of the “elliptic modulus” k. Some software, such as Mathematica and SciPy, define K in terms of “the parameter” m = k². So to evaluate K above using Mathematica, you’d need to square the argument first.

    g = EllipticK[0.5] Sqrt[2]/Pi

    g = 2 EllipticK[-1]/Pi

Similarly, in Python we’d write

    from scipy.special import ellipk
    from math import pi
    g = ellipk(0.5)*2**0.5/pi

    g = 2*ellipk(-1)/pi

Foreshadowing hypergeometric functions

Posted on 28 October 2021 by John

I’ve said in several blog posts that multi-factorials come up often in practice, and given examples. This post will give a glimpse of why this is.

Rising powers

The kth rising power of a is

(a)_k = a (a+1) (a + 2) … (a + k – 1).

So, for example, (5)₃ = 5*6*7 and (1)_k = n!.

Double factorials

Now let’s look at rising powers of fractions.

(1/2)_k = (1/2) (3/2) (5/2) … (2k-1)/2 = (2k-1)!! / 2^k.

This means double factorials of odd numbers can be written as rising powers.

(2k-1)!! = 2^k (1/2)_k .

Now for even numbers,

(2k)!! = 2^k k! = 2^k (1)_k

and so even double factorials can be written in terms of rising powers too. This means that the term

(2k – 1)!! / (2k)!!

from the series in the previous post could be written as

(1/2)_k / (1)_k .

Hmm. Ratios of rising powers sound useful.

Generalizing

From the discussion above you can imagine how to create triple factorials out of rising powers of 1/3 and 2/3. In fact it’s easy to see that all multifactorials can be made out of rising powers of rational numbers.

Ratios of multifactorials come up often in coefficients of power series, and these can be rewritten as ratios of rising powers. Maybe we should have a name for functions whose coefficients are ratios of rising powers. And we do: hypergeometric functions. You can find the exact definition here.

At first the definition of a hypergeometric function seems arbitrary. I hope my previous posts, and especially this post, has been a sort of warm up to motivate why you might be interested in such functions.

Applications

In an old post I quoted W. W. Sawyer saying

There must be many universities today where 95 percent, if not 100 percent, of the functions studied by physics, engineering, and even mathematics students, are covered by the single symbol F(a, b; c; x).

The notation F(a, b; c; x) here denotes the hypergeometric function whose kth power series coefficient is

(a)_k (b)_k / (c)_k k!

More general hypergeometric functions have two lists of parameters, separated by a semicolon. The first set corresponds to rising powers in the numerator, and the second set corresponds to rising powers in the denominator.

The exponential function is the simplest hypergeometric function: both parameter lists are empty, leaving only the k! term in the denominator.

The arcsine function mentioned above is an example of a function that is not hypergeometric per se, but can be written simply in terms of a hypergeometric function:

sin^-1(x) = x F(1/2, 1/2; 3/2, x²).

This is more what Sawyer had in mind. Elementary functions like sine, cosine, and log note pure hypergeometric functions but are simply related to hypergeometric functions, similar to arcsine. The same is true of more advanced functions, like elliptic integrals and the CDFs of common probability distributions.

Quadruple factorials and Legendre functions

Posted on 18 October 2021 by John

Last week I claimed that double, triple, and quadruple factorials came up in applications. The previous post showed how triple factorials come up in solutions to Airy’s equation. This post will show how quadruple factorials come up in solutions to Legendre’s equation.

Legendre’s differential equation is

$(1-x^2) y'' - 2x y' + \nu(\nu+1)y = 0$

The Legendre functions P_ν and Q_ν are independent solutions to the equation for given values of ν.

When ν is of the form n + ½ for an integer n, the values of P_ν and Q_ν at 0 involve quadruple factorials and also Gauss’s constant.

For example, if ν = 1/2, 5/2, 9/2, …, then P_ν(0) is given by

$(-1)^{(2\nu - 1)/4} \frac{\sqrt{2}(2\nu - 2)!!!!}{(2\nu)!!!!\pi g}$

Source: An Atlas of Functions

Rising powers

Why the example isn’t obvious

Rational polynomials

Back to our example

Spot check

Related posts

Related posts

Branches 0 and -1

SciPy and Mathematica

Branch numbering

Recovering k

Partition example

Lambert’s W function

Basic bootstrapping

Putting it all together

More Lambert W posts

Motivation

q-factorials and q-binomial coefficients

Plot twist

Related posts

Power series

Asymptotic series

Confluent hypergeometric series

Gauss’s constant

Elliptic integrals

Hypergeometric functions

Conventions

Rising powers

Double factorials

Generalizing

Applications

More on hypergeometric functions