Gamma and Factorial in the Monthly

TL;DR

This paper surveys fifty publications from The Monthly on the Gamma function and Stirling's formula, contextualizing them within the broader mathematical literature to highlight key developments and insights.

Contribution

It provides a comprehensive overview of the literature on the Gamma function and Stirling's formula published in The Monthly, emphasizing notable works and their significance.

Findings

Identifies key papers on Gamma and Stirling's formula in The Monthly.

Places these papers within the broader mathematical context.

Highlights influential contributions and trends in the literature.

Abstract

The Monthly has published roughly fifty papers on the $\Gamma$ function or Stirling's formula. We survey those papers (discussing only our favourites in any detail) and place them in the context of the larger mathematical literature on $\Gamma$.

Equations

$$\Gamma (z)=(z-1)!={\int }_0^{\infty }{t}^{z-1}{\mathrm{e}}^{-t}dt\mathrm{for}\mathrm{Re}z>0,$$

$$(-z)!z!=\frac{\pi z}{\sin \pi z}\mathrm{or}\Gamma (1-z)\Gamma (z)=\frac{\pi }{\sin \pi z}$$

$$z!$$

$$\frac{1}{\Gamma(z)}=z\mathrm{e}^{\gamma z}\prod_{k\geq 1}\bigg{(}1+\frac{z}{k}\bigg{)}\mathrm{e}^{-z/k}\>,$$

$$\psi (z)=\frac{d}{dz}\ln \Gamma (z)=\frac{{\Gamma }^{\prime }(z)}{\Gamma (z)}.$$

$$\Gamma (2z)$$

$$\zeta (z)$$

$$\ln \Gamma (z)$$

$$B(s,t)={\int }_{x=0}^1{x}^{s-1}(1-x{)}^{t-1}dx=\frac{\Gamma (s)\Gamma (t)}{\Gamma (s+t)}.$$

$$x(1+y)y^{\prime }-y=0.$$

$$H(z)=\frac{1}{\Gamma(1-z)}\frac{d}{dz}\ln\Bigg{[}\Gamma\bigg{(}\frac{1-z}{2}\bigg{)}\bigg{/}\Gamma\bigg{(}1-\frac{z}{2}\bigg{)}\Bigg{]}\>.$$

$$H(z)=\frac{\psi\bigg{(}1-\frac{z}{2}\bigg{)}-\psi\bigg{(}\frac{1-z}{2}\bigg{)}}{2\Gamma(1-z)}\>,$$

$$L(x)=-\gamma x+\sum_{k\geq 1}\Bigg{(}\frac{x}{k}-\ln\bigg{(}1+\frac{x}{k}\bigg{)}\Bigg{)}\>,$$

$$\frac{1}{\Gamma (s)}=\frac{1}{2\pi }{\int }_{\mu =-\infty }^{\infty }\frac{e^{a+i\mu }}{(a+i\mu {)}^s}d\mu$$

$$\frac{d^2}{d{z}^2}\ln f(z)=\sum _{n\ge 0}\frac{1}{(n+z{)}^2}.$$

$$\ln \Gamma (z)={\zeta }^{\prime }(0,z)-{\zeta }^{\prime }(0)$$

$$\zeta (s,z)=\sum _{k\ge 0}(k+z{)}^{-s}$$

$$\bigg{(}1-\frac{a}{x+a}\bigg{)}^{1-a}\leq\frac{\Gamma(x+a)}{x^{a}\Gamma(x)}\leq 1$$

$$\Gamma\bigg{(}\frac{1}{62}\bigg{)}\Gamma\bigg{(}\frac{33}{62}\bigg{)}\Gamma\bigg{(}\frac{35}{62}\bigg{)}\Gamma\bigg{(}\frac{39}{62}\bigg{)}\Gamma\bigg{(}\frac{47}{62}\bigg{)}=2^{4}\pi^{5/2}\>.$$

$$\frac{\Gamma\bigg{(}n+\dfrac{1+u}{2}\bigg{)}}{\Gamma\bigg{(}n+\dfrac{1-u}{2}\bigg{)}}\sim\bigg{(}n^{2}+\dfrac{1-u^{2}}{12}\bigg{)}^{u/2}$$

$${ϵ}_n(u)=\frac{u(1-u^2)(4-u^2)}{6!n^4}{F}_n(u)$$

$$\frac{d^2}{d{z}^2}\ln \Gamma (z)$$

$$=\sum _{\ell \ge 0}\frac{\ell !}{(\ell +1)}\cdot \frac{1}{z^{\overline{\ell }}}\mathrm{Re}z>0$$

$$F\left(\begin{array}{c}a,b\text{hfill}\\ c\end{array}|1\right)=\sum _{j\ge 0}\frac{a^j{b}^j}{c^j}\cdot \frac{1}{j!}=\frac{\Gamma (c)\Gamma (c-a-b)}{\Gamma (c-a)\Gamma (c-b)}$$

$$\frac{x}{\ln (1+x)}=1+L_{1}x+L_2{x}^2+\cdots$$

$$L(t)={\int }_{s=0}^1\left(\begin{array}{c}s\text{hfill}\\ t\end{array}\right)ds,$$

$$\frac{1}{\Gamma (z)}=\frac{1}{2\pi i}{\int }_C{s}^{-z}{\mathrm{e}}^{s}ds.$$

$$y=f(x)$$

$$q=1/f(p)$$

$$\ln \Gamma (z)$$

$$r$$

Full text

Gamma and Factorial in the Monthly

Jonathan M. Borwein and Robert M. Corless

1 Introduction

The Monthly has published roughly fifty papers on the $\Gamma$ function or Stirling’s formula. We survey those papers (discussing only our favourites in any detail) and place them in the context of the larger mathematical literature on $\Gamma$. We will adopt a convention of $z\in\mathbb{C}$, $n\in\mathbb{N}$, $x\in\mathbb{R}$, and often $x>0$. We begin with notation:

[TABLE]

which extends the positive integer factorial $n!=n(n-1)\cdots 2\cdot 1$ to the right half of the complex plane. The reflection formulas (discovered by Euler)

[TABLE]

explicitly show the analytic continuation and identify the poles at $z=-n$, nonpositive integers (for $\Gamma$), and negative integers (for factorial) with residues $(-1)^{n}/n!$ for $\Gamma$ and an equivalent shifted result for factorial.

We reluctantly bypass the amusing “notation war” where authors argue about the “minor but continual nuisance” of the shift by 1 in passing between factorial notation and $\Gamma$ notation, $z!=\Gamma(1+z)$. See for example the footnotes in **[26]**, **[46]**, **[51]**, the introduction in **[57]**, and the history in **[44]**, for some fullisades.

An important fact is that $\Gamma(x)$ is logarithmically convex, i.e. that $\ln\Gamma(x)$ is convex (see figure 1(b)) and that this distinguishes $\Gamma$ amongst all interpolants of the factorial. This is the Bohr-Mollerup-Artin theorem. A function that satisfies $f(z+1)=zf(z)$ and $f(1)=1$—but is not logarithmically convex—is called a pseudogamma function.

Euler’s definition of $z!$ is, for $z$ not a negative integer,

[TABLE]

We have $z!=z(z-1)!$ and $\Gamma(z+1)=z\Gamma(z)$. Weierstrass apparently preferred Newman’s formula:

[TABLE]

where $\gamma=\lim\limits_{n\rightarrow\infty}(\sum_{k=1}^{n}\nicefrac{{1}}{{k}}-\ln n)$ is the Euler-Mascheroni constant, $\gamma=0.577\ldots\>.$

The derivative of $\Gamma$ is less convenient than is the combination

[TABLE]

The function $\psi(z)$ is also called the digamma function.

Legendre derived the duplication formula:

[TABLE]

The reflection formula

[TABLE]

was known to Euler, but was first proved by Riemann. Here $\zeta(z)=\sum_{n\geq 1}1/n^{z}$ is the famous Riemann zeta function. Kummer (1847) derived a Fourier series:

[TABLE]

which is illustrated in figure (1(a)). Also, the Beta function is

[TABLE]

2 The $\Gamma$ and Factorial Functions: A cross–section of mathematics

The following Monthly papers (in alphabetical order by surname of the first author) are the ones that stand out, our “eleven favourites”111In [13] the authors highlighted what they called “$\pi*$” articles, i.e. those that had been cited at least $30$ times. In contrast, inclusion in this list here is purely subjective.:

Patrick Ahern and Walter Rudin, “Geometric Properties of the $\Gamma$ function” **[2]** 2. 2.

Wladimir de Azevedo Pribitkin, “Laplace’s Integral, the Gamma Function and Beyond” **[6]** 3. 3.

Richard Askey, “Ramanujan’s extensions to the $\Gamma$ and Beta functions” **[4]** 4. 4.

Bruce Berndt, “The $\Gamma$ function and the Hurwitz $\zeta$-function”  **[10]** 5. 5.

Manjul Bhargava, “The factorial function and generalizations” **[11]** 6. 6.

Philip J. Davis, “Leonhard Euler’s Integral” **[26]** 7. 7.

Louis Gordon, “A stochastic approach to the $\Gamma$ function” **[43]** 8. 8.

Detlef Laugwitz and Bernd Rodewald, “A simple characterization of the $\Gamma$ function” **[59]** 9. 9.

Reinhold Remmert, “Wielandt’s theorem about the $\Gamma$ function” **[81]** 10. 10.

Lee Rubel, “A survey of transcendentally transcendental functions”  **[85]** 11. 11.

Gopala Krishna Srinivasan, “The Gamma function: An eclectic tour” **[92]**

These eleven disparate papers provide a cross–section of mathematics: analysis, geometry, statistics, combinatorics, logic and number theory. We begin with the earliest, the Chauvenet–winning **[26]** by Philip J. Davis. We then tour through the others and end with Fields medallist Manjul Bhargava’s **[11]**.

Fundamental Properties

Davis tells the story of $\Gamma$ as an instance of how mathematics evolves. This required informed historical commentary and an appreciation of mathematical aesthetics. The paper was dedicated to the memory of Milton Abramowitz, who died (too young) in 1958. Davis was then chief of Numerical Analysis at The National Bureau of Standards of the United States, and a contributor to the great project that became (after six years of further hard work by Irene A. Stegun) the monumental Handbook of Mathematical Functions **[1]**. Indeed the graph of the real $\Gamma$ function in **[26]** is apparently the same as in the Handbook. This is not a surprise: Davis wrote that section.

This monumental work has recently been updated as The Digital Library of Mathematical Functions (http://dlmf.nist.gov) **[76]**. The DLMF and the similar-but-different INRIA project, Bruno Salvy’s Dynamic Dictionary (http://ddmf.msr-inria.inria.fr/1.9.1/ddmf) **[87]**, fill a need for good online mathematical material.

To return to $\Gamma$, the paper **[26]** gives a story that enriches the plain facts in the Handbook, and illustrates the evolution of mathematical thought. Davis spent a lot of space on the interpolation of $n!$ by $\Gamma(z+1)$, treating it with the philosophical attention it deserves. He did, however, note that $\Gamma$ arises in applications mostly because the integral defining it does.

The impressive paper **[43]** by L. Gordon proves many of these same facts by considering the $\Gamma$ distribution and taking the statistical point of view. Of course, $\Gamma$ is widely used in probability and statistics.

Now we come to the survey **[85]** by L. A. Rubel, which is concerned with $\Gamma$ only insofar as $\Gamma$ is the most famous of the “transcendentally transcendental” (TT) functions; that is, functions that do not satisfy any polynomial differential equation $P(x;y,y^{\prime},\cdots,y^{n})=0$. Functions that do satisfy such equations are called Differential Algebraic functions; for instance, the Lambert $W$ function **[24]** satisfies $y\mathrm{e}^{y}=x$ and thus is DA because

[TABLE]

It’s mentioned in **[26]** that the first proof that $\Gamma$ is TT was due to Hölder in 1887.

Rubel economically and intelligibly presents Ostrowski’s proof that $\Gamma$ is TT. The proof has many of the characteristics one sees in computer algebraic proofs: ordering terms in reverse lexicographic order (Rubel doesn’t call it that), and using divisibility arguments to arrive at a minimalist polynomial $P$ for $\Gamma$, which can then be shown self-contradictory. The only property of $\Gamma$ aside from differentiability that is used is the functional equation $\Gamma(x+1)=x\Gamma(x)$. Therefore, the proof runs through unchanged for (sufficiently differentiable) pseudogamma functions. In particular, Hadamard’s entire pseudogamma $H(z)$ (see equation (12) below), is TT. This was not pointed out in **[26]** and it’s not obvious from its definition (after all, the combination $\Gamma(x+1)/\Gamma(x)$ is algebraic; it’s not automatic that combinations of $\Gamma(\nicefrac{{(1-x)}}{{2}}),\Gamma(1-\nicefrac{{x}}{{2}})$ and $\Gamma(1-x)$ will be TT). Altogether this highly-cited paper provides a very readable introduction to the area, and fills a gap in the education of people who use $\Gamma$.

Characterizing $\Gamma$

The paper **[59]** by Laugwitz and Rodewald was a suprise to RMC, as was the paper **[81]** by Remmert. RMC had not been aware of satisfactory alternatives to the Bohr-Mollerup-Artin characterization of $\Gamma(x)$. Such a condition is necessary, because there are pseudogamma functions; for instance, Davis gives Hadamard’s example (“Hadamard’s pseudogamma function”):

[TABLE]

Maple simplifies this quite dramatically to

[TABLE]

where $\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$. In retrospect, it’s not so dramatic, merely replacing a derivative with $\psi$; nonetheless it looks simpler. For a more nuanced discussion of the satisfaction given by such formulae, see **[104]** by Wilf.

Curiously enough $H(z)$ is an entire function. Why, then, is $\Gamma(z)$ to be preferred to $H(z)$, given that $\Gamma$ has poles at non-positive integers? The “real” answer given by Davis is that the integral for $\Gamma$ appears frequently. But the logarithmic convexity characterization of Bohr-Mollerup-Artin is compelling. Davis eloquently explains “The desired condition was found in notions of convexity. A curve is convex if the following is true of it: take any two points on the curve and join them by a straight line; then the portion of the curve between the points lies below the line. A convex curve does not wiggle; it cannot look like a camel’s back. At the turn of the century, convexity was in the mathematical air$\ldots$” **[26, p. 867]**.

We point out that logarithmic convexity is stronger than mere convexity: the convexity of $f(x)$ does not imply the convexity of $\ln f(x)$, but the convexity of $\ln f(x)$ implies that of $f(x)$.

But there are alternative characterizations. The characterization discussed by Laugwitz and Rodewald in **[59]** is actually due to Euler himself, and is very elegant: “For a fixed infinitely large natural number $n$ and all finite natural numbers $m$, the expression $(n+m)!$ behaves approximately like a geometric sequence.” **[59, p.534]**. Laugwitz and Rodewald go on to put this in a more usual ‘limit’ language: they prove that if $L(0)=0,L(x+1)=\ln(x+1)+L(x)$, and $L(n+x)=L(n)+x\ln(n+1)+r_{n}(x)$ where $r_{n}(x)\to 0$ as $n\to\infty$ (here $x\geq 0$ and $n\in\mathbb{N}$) then $L(x)$ exists and is unique and is thus $\ln(\Gamma(x+1))$. They do this by showing that the conditions imply a convergent series for $L(x)$; they then verify that the resulting series,

[TABLE]

satisfies the conditions.

Another alternative characterization is apparently due to Wielandt. Until Remmert’s paper **[81]**, this seems to have been “hardly known”. Wielandt’s theorem states that if $F(z)$ is holomorphic in the right half plane, $F(z+1)=zF(z)$ for all such $z$, and is bounded on the strip $1\leq\operatorname{Re}z\leq 2$, then $F(z)=a\Gamma(z)$ and $a=F(1)$. This characterization seems to be quite easy to use, more so than logarithmic convexity.

The intriguing and concise paper **[2]** by Ahern and Rudin cites **[81]** but does not appear to actually use Wielandt’s theorem; instead it extends logarithmic convexity into the complex plane. The final theorem in the paper states that $\ln\Gamma(z)$ is univalent in $\operatorname{Re}\ z>x_{0}$, where $x_{0}$ is the only positive zero of $\psi(x)$. They note $1<x_{0}<2$; by a computer algebra system, $x_{0}=1.4616321\ldots$. There seems to be no standard name for this number. This result on univalence seems to be as close as it gets in the Monthly to a discussion of the functional inverse of $\Gamma$ (see later in this section).

The very interesting (and fun to read) paper **[6]** by W. de Azevedo Pribitkin does use Wielandt’s theorem. The paper is concerned with the Laplace integral for the reciprocal of $\Gamma$:

[TABLE]

where $a>0$ is any positive real number. This paper comes as close as any in the Monthly to the computation scheme of  **[89]**.

Srinivasan’s “Eclectic Tour” **[92]** also uses Wielandt’s theorem. This wide-ranging survey achieves a remarkable complementarity to Davis’ paper; indeed even the historical details are different (not contradictory, just complementary). New elegant proofs of many results are attained using a function-theoretic mindset and a systematic approach, apparently termed the “additive approach” by Remmert. The paper starts with the differential equation (which is also in **[26]**)

[TABLE]

With initial conditions $f(1)=1$ and $f^{\prime}(1)=-\gamma$, this gives $f(z)=\Gamma(z)$. Wielandt’s theorem is used in the proof. Altogether, this is a remarkable paper, which presents (with proofs!) a very significant body of knowledge of $\Gamma$, and of related functions.

Relationship to $\zeta$

Bruce Berndt’s short paper **[10]** illuminates old connections of $\Gamma$ to the Hurwitz $\zeta$ function. He gives several proofs using that connection. For instance, he proves equation (8). He first proves

[TABLE]

where ′ means differentiation with respect to $s$ (the first variable) and

[TABLE]

is the Hurwitz zeta function; $\zeta(s)=\zeta(s,1)$ is the Riemann zeta function.

Generalizations

We begin with Askey’s beautiful paper **[4]**. Askey concurred with Davis on the reason $\Gamma$ is useful: “Euler’s integral (1.3) occurs regularly and is the real reason for studying the gamma function.” **[4, p. 347]**. Other than that, Askey’s paper is quite different from Davis’s, although it, too, tells a story, namely the story of Ramanujan’s $q$-extensions of the $\Gamma$ and Beta functions. Askey provides much detail, many results, and concise proofs. He places the results in context, and explains why they are important. His starting point was that Hardy thought the results strange; Askey demonstrates that the results are profoundly connected to modern developments, and thus not so ‘strange’ today. Of course, part of the direction of modern developments is due to Ramanujan, so this conclusion might have been expected. Both Davis and Askey show relationships of $\Gamma$ to other functions and areas.

Now we come to the paper **[11]** by Manjul Bhargava. The first time RMC read the paper, RMC didn’t notice the author’s name222In his defence, it was in a BIG pile of $\Gamma$ papers!; RMC noted down that the paper was concerned with discrete generalizations of the factorial, was focused on number theoretic ideas and applications, and seemed quite substantial, much more so than some other “generalization” papers in the Monthly, namely **[73]** and **[35]** (although the last isn’t so lightweight, either).

On RMC’s next pass, he paid a bit more attention. Bhargava was awarded the Fields medal in 2014 for new methods in geometry of numbers. It seems a safe bet that this paper, although written 14 years previous to his prize, might perhaps contain something a bit more substantial than the average paper in the Monthly on the $\Gamma$ function333It’s interesting that the notion “an average paper in the Monthly on $\Gamma$” actually makes sense.. Of course, it does, even though it is written in a lively, almost joyful, style.

Instead of making a small change in the definition of $\Gamma$ (in contrast, the paper **[73]** by T. A. Newton changes integrals to sums, and derivatives to differences; the paper **[35]** by Tomlinson Fort changes the “$+1$” in the functional equation to a variable “$+h_{n}(x)$”) Bhargava isolates an invariant that characterizes the factorial: he writes $k!$ as the product over all primes $p$ of certain valuations. Then he is able to generalize this to work over subsets of the integers or even over more abstract objects. To do this he used something reminiscent of the Leja ordering for good-quality use of Newton interpolation, and indeed the bases he constructs are very like Newton interpolational bases. This generalization encompasses the $q$-factorials from enumerative combinatorics which we met in **[4]** and it seems the ones in **[73]** too, as well as others in the non-Monthly literature.

Bharghava ends his paper with a sequence of then-open questions, including questions on generalizations of this factorial to a generalized $\Gamma$ function, and on generalizations of Stirling’s formula.

Shorter papers and notes

There are quite a few short pieces on $\Gamma$ in the Monthly. One of the oldest papers in the Monthly is **[7]**, which uses explicit expressions for orthogonal polynomials and their orthogonality conditions to arrive at certain identities for binomial double sums.

The Note **[102]** uses Hölder’s inequality to establish

[TABLE]

for $0<a<1$ and $x>0$; letting $x\to\infty$ gives the classical $\Gamma(x+a)\sim x^{a}\Gamma(x)$, which is an asymptotic generalization of the functional equation $\Gamma(x+1)=x\Gamma(x)$. Nanjundiah’s paper **[69]** contains an interesting trick, using the asymptotic relation $B(\alpha,t)\sim\Gamma(\alpha)/t^{\alpha}$ as $t\to\infty$ between the Beta function and the $\Gamma$ function to get the exact relationship, $B(x,y)=\Gamma(x)\Gamma(y)/\Gamma(x+y)$ and from there through $\Gamma(x+a)\sim x^{a}\Gamma(x)$ (again!) to Euler’s product definition.

The very pretty paper **[74]** cites **[14]** and a Monthly problem **[38]** as motivation, and produces a recipe for finding some short products of Gammas that give simple values, such as

[TABLE]

The paper **[17]** merely finds $\Gamma$ functions in a particular solution to the strange (nonlinear) ODE $\frac{d^{n}y}{dx^{n}}\frac{d^{m}x}{dy^{m}}=1$.

The paper **[54]** by Knuth probably shouldn’t be included in this review, although the ratio $n!/n^{n}$ plays a starring role: this is Egorychev’s theorem, that the permanent (like the determinant, but with all plus signs instead of $(-1)^{\sigma}$ in the permutation expansion) of a doubly stochastic matrix cannot be less than this, and therefore the doubly stochastic matrix $\mathbf{A}$ with $a_{ij}=1/n$ has this minimal permanent. Egorychev had settled a long-standing conjecture of Pólya’s with this theorem, and the purpose of **[54]** was pure exposition. The paper also demonstrates that this review cannot be expected to be complete—the paper was not found by computer search, but rather by browsing, i.e. by serendipity.

Speaking of serendipity, the paper **[83]** uses an exploration of something related to fractional integration to discuss just that, in mathematical practice. The notion of fractional integration, i.e. interpolation in the dimension of a multiple integral, is an idea of the same type that led to the creation of $\Gamma$, of course, as discussed in **[26]**.

The paper **[36]** is concerned with the ratio and its approximation

[TABLE]

for large $n$. This formula is exact when $n\geq 1$ is an integer and $u=0,\pm 1,\pm 2$. The rest of the paper is concerned with the error, which the author shows is exp$(-\epsilon_{n}(u))$* with*

[TABLE]

where $0<F_{n}(u)\leq 1$ for $|u|\leq 1$. The paper **[50]** is also concerned with $\Gamma$ function ratios.

The expression (15) has a similar analogue in rising powers, namely

[TABLE]

where $z^{\overline{\ell}}=z(z+1)\cdots(z+\ell-1)$, $z$ to the $\ell$ rising. This equation can be derived from a formula of Gauss for a hypergeometric function with unit argument,

[TABLE]

if $\operatorname{Re}(c-a-b)>0$ and $c\not\in\mathbb{N}$. The paper **[86]** gives some details of a proof.

The paper **[47]** uses the $\Gamma$ function to find volumes of spheres in higher dimension.

The paper **[25]** by H. T. Davis considers the generating function for so–called logarithmic numbers $L_{n}$:

[TABLE]

and establishes the curious interpolation formula $L(t)=L_{n}$ when $t=n$, where

[TABLE]

among other results. H. T. Davis does not use the binomial notation, but rather a ratio of $\Gamma$ functions, which makes clear the meaning of the neater binomial coefficient used here.

Finally, a referee points out the delightful paper **[58]** in which Lange uses properties of $\Gamma$ to give a new continued fraction for $\pi$. Since $\Gamma$ didn’t appear in the title, we missed that one!

Gaps: computation, visualization, inverse $\Gamma$

One purpose of a survey such as this is to identify gaps. We noticed three. There’s not much visualization for $\Gamma$ in the Monthly, apart from in **[26]**. There’s almost no computational work. The most surprising omission, though, is the lack of exploration of the functional inverse of $\Gamma$, which we denote here by $\mathrm{inv}\Gamma$ or $\check{\Gamma}$. To be fair, we have found only four references in total on $\mathrm{inv}\Gamma$, and none in the Monthly. Still, to have studied a function since 1730, or 1894 from the start of the Monthly, without thinking of the inverse? It seems quite the blind spot.

Visualization

There has been significant progress in comprehensible visualization of complex functions, by the use of color phase plots **[100, 101]**. In the inset figure in figure 3(a) the $z$-plane is colored by arg$(z)=\mathrm{invtan}(y,x)$; the same colors are used in figure 3(a) for coloring by arg$(\Gamma(z))$. Of course, $z=x+iy$. The poles at [math], $-1$, $-2$, and so on are clearly visible. In figure 3(b) we’ve added the same phase coloring to the plot of $|\Gamma(z)|$, on the same scale as the famous hand-drafted plot of **[49]**. If one looks carefully, this figure illustrates Wielandt’s boundedness characterization: because the contours in the right half bend away from $y$ axis, the boundedness on (say) $1\leq x\leq 2$ is entirely believable.

A similar pair of plots are given for equation (13), in figure 4. We will see more visualizations shortly, including one for a portion of the Riemann surface for $\mathrm{inv}\Gamma$.

Computation

We did not find in the Monthly any substantive discussion (or even pointers) as to how to compute the $\Gamma$ function, apart from using the fundamental recurrence to push $z$ far enough into the right half-plane that Stirling’s asymptotic formula can be used. For a solid discussion of good methods, see **[57]**, **[89]**, and **[91]**. Schmelzer and Trefethen’s **[89]** uses the beautiful contour integral

[TABLE]

Here the contour C, called a deformed Bromwich contour, winds around the negative real axis counterclockwise. It turns out that the trapezoidal rule (or the midpoint rule) is spectrally accurate for this quadrature, and the contour integral is (almost) entirely practical! See **[89*]** for details. *

Inverse $\Gamma$

Plotting the real inverse Gamma function

As stated earlier, it’s surprising in hindsight that there’s so few (and only so recent) studies on the inverse Gamma function, $\mathrm{inv}\Gamma$, and none in the Monthly. We begin with the simplest things. Plotting inverse functions with computer software is very simple, even if there are no routines for the inverse function itself. The basic idea is explained in **[90]**: just exchange $x$ and $y$. Figure 5(b) was produced using Maple. The details of the code are available at http://www.apmaths.uwo.ca/~rcorless/frames/Gamma/GammaBorweinCorless.mw.

In figure 5(a) we see several locations where $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ is zero: between $x=1$ and $x=2$, (at about $x_{0}=1.4616\ldots$), between $x=-1$ and $x=0$ (at about $x_{-1}=-0.5048\ldots$), between $x=-2$ and $x=-1$ (at about $x_{-2}=-1.5734\ldots$), and so on. Denote $x_{-k}$ the root between $-k$ and $-k+1$, except $x_{0}\in(1,2)$. Each of these locations separates a pair of branches of $\mathrm{inv}\Gamma$ in figure 5(b). The branch with $\mathrm{inv}\Gamma(x)\geq x_{0}$ is of special interest, because it contains the inverse factorials: $\Gamma(n)=(n-1)!$ means that $\mathrm{inv}\Gamma(24)=5$, for instance, on this branch. We take up the asymptotics of this branch of $\mathrm{inv}\Gamma$ in section 3.

Plotting the Riemann Surface for the Inverse $\Gamma$ Function

Given a reliable routine to compute a complex function, it is possible to visualize the Riemann surface for inv$f$* automatically [23, 95, 94]. The basic idea is the same as above. If $z=x+iy$ and $f(z)=u+iv$, then a Riemann surface for inv*$f$, the functional inverse of $f$, can be plotted by plotting $[u,v,y]$ as a parametric plot with $x$ and $y$ as parameters. Note that a routine to compute inv$f$* is not required.*

There is one theoretical consideration, namely how to ensure that each point on the representation corresponds uniquely to a point on the four-dimensional map $(x,y)\leftrightarrow(u,v)$, which isn’t always possible because not all Riemann surfaces are embeddable in 3D, but this can be mitigated by the use of color as the fourth dimension. There is a practical consideration, as well: namely how to choose the portion of the surface to fit into the viewing frame. For the $\Gamma$ function, there is a further difficulty, namely its extreme dynamic range, even in its imaginary part.

Because the functional inverse of the reciprocal of $f$ is related simply to the functional inverse of $f$, as follows,

[TABLE]

we decided to plot a portion of the Riemann surface for $\mathrm{inv}\Gamma(1/z)$. See figure 6. The multilayer convolutions show a rich structure for investigation. Our code is archived at http://www.apmaths.uwo.ca/~rcorless/frames/Gamma/GammaBorweinCorless.mw.

As a practical matter, constructing branches of $\mathrm{inv}\Gamma$ that people can use seems to be at least as important. A colleague and RMC have completed some preliminary investigations, but we’ll leave reporting on this to future work. We note the foundations provided by **[78]**.

3 Stirling’s Formula

Judging by the number of articles in the Monthly on the subject, Stirling’s formula approximating $n!$ for large $n$ is by far the most popular aspect of the $\Gamma$ function. There are “some remarks”, “notes”, more “remarks”; there are “simple proofs”, “direct proofs”, “new proofs”, “corrections”, “short proofs”, “very short proofs”, “elementary” proofs, “probabilistic” proofs, “new derivations”, and (our favourite title) “The (n+1)th proof”. The earliest one we can find is **[48]**; the paper **[84]** is about a different formula444Which could be applied to $z!$ and $\log\Gamma(1+z)$. See [26]..

The name “Stirling’s formula” turns out not to be historically correct, as we will see. However, there will be no changing the convention. Here is what we mean conventionally by Stirling’s formula:

[TABLE]

where $\Gamma$ is asymptotic to a divergent series

[TABLE]

and the $B_{2n}$ are Bernoulli numbers, with generating function

[TABLE]

Stirling’s formula also occurs in an exponentiated form, $\Gamma(z)\sim\sqrt{2\pi}z^{z-\nicefrac{{1}}{{2}}}\mathrm{e}^{-z}\mathrm{e}^{r}$ and one can easily compute a series for $\mathrm{e}^{r}$ from (30). It turns out that

[TABLE]

with $a_{1}=1$ and for $n>1$

[TABLE]

Note that only the odd $n$ occur in (31) but that even $n$ are needed in the recurrence.

Most people don’t bother to distinguish equation (28) from its exponential: both are useful, and both are called Stirling’s formula or Stirling’s series. Note that equation (28) contains only odd powers and thus seems more efficient than equation (31): $\mathrm{e}^{r}$ contains all powers. Exponentiation of a series is simple. See e.g. **[22, ch. 2]**.

Trapezoidal Rule analyses

Hummel motivates the paper **[48]** by claiming

“The average college undergraduate, with nothing beyond an elementary calculus course, knows little about infinite series and still less about infinite products, hence gets very little out of one of the usual developments of Stirling’s formula.”

Students were different, then, as was the curriculum. Nevertheless the conclusions may still be true. Hummel then goes on in **[48]** to use the concavity of logarithm and its consequential inequality from the trapezoidal rule to establish

[TABLE]

and claims $\mathrm{e}^{11/12}\doteq\sqrt{2\pi};$ since $\mathrm{e}^{11/12}\doteq 2.5009$ while $\sqrt{2\pi}\doteq 2.5066$ this claim is entirely reasonable. The upper bound is not very tight.

Several Monthly papers such as Nanjundiah’s **[68]** and the highly-cited **[82]** by Robbins use a similar starting point, though few refer to **[48]**. The paper **[107]** uses it backwards to show Stirling’s formula implies the integral for $\ln x$. Using the trapezoidal rule is a natural idea: the trapezoidal rule gives a lower bound on an integral of a concave function. The paper **[52]** uses the trapezoidal rule to establish the existence of the limit—this is not a petty point. The paper **[20]** uses an average of trapezoids to get both a lower bound and an upper bound for $n!$.

Here is the basic element of the trapezoidal analysis: because $\ln t$ is concave, the trapezoidal rule gives a lower bound for the integral:

[TABLE]

We now do our own version. Adding up panels of width $1$ from $t=m$ to $t=n$, where $m<n$,

[TABLE]

Adding $\sum^{m}_{j=1}\ln j-\ln(m)/2+\ln(n)/2$ to both sides,

[TABLE]

or

[TABLE]

where

[TABLE]

Exponentiating,

[TABLE]

We’ve thus replaced the constant $c=\sqrt{2\pi}$ in Stirling’s formula ** by a computable upper bound. The larger we take $m$, the tighter the bound. When $m=1$ this is the upper bound quoted in **[48]**. Taking $m=n-1$ gives the tightest possible estimate, but is silly, because to compute $d_{n-1}$ requires $(n-1)!$. We will consider $m\geq n$ later. Taking $m$ small or of moderate size is reasonable in order to estimate $n!$ for large $n$. Taking $m=10$ gives $d_{10}\doteq 2.5276$, whereas $\sqrt{2\pi}\doteq 2.5066$. Taking $m=20$ gives $d_{20}\doteq 2.5171$. More work shows $d_{m}$ monotonically decreases.

Binet’s Formula

We pass so easily from $n!$ to $z!$ nowadays that it’s hard to notice the steps in passing from the asymptotics of $n!$ to the asymptotics of $z!$, and of course it’s obvious that the formula must be the same555Some proofs do appear in the Monthly, e.g. Blyth and Pathak [12] use $\Gamma(n+\alpha)\sim n^{\alpha}\Gamma(n)$ to do so.. Nonetheless , it’s satisfying to read in **[103]** a proof using Binet’s formula, namely (for $\operatorname{Re}z>0$)

[TABLE]

that they really are the same. In the Monthly, we find **[88]** which gives a similar proof of (40) and an intelligible recap of its use to derive the logarithmic variant (28) of Stirling’s formula, plus bounds guaranteeing the errors in truncation.

Other papers

One of the prettiest papers is **[3]**. The author, M. I. Aissen, gives a concise but lucid exposition that uses the simple inequality $n!<n^{n}$ and the properties of the ratios $u_{n}=n^{n}/n!$ to make the appearance of $\mathrm{e}$ seem inevitable and completely natural. From there the author gives two proofs of Stirling’s formula, up to $\sqrt{2\pi}$.

The shortest paper that we found is **[80]**: in just 17 lines, including title, references, and address, Pippenger establishes the following “striking companion to Wallis’ product”

[TABLE]

by appeal to Stirling’s formula.

The papers **[69]** and **[62]** give simple refinements of the result from the previously mentioned **[82]** by Robbins.

The author of **[32]**, a famous statistician from Princeton, was “deeply apologetic” in **[33]** for various errors. We point this out not in a spirit of schadenfreude but rather in sympathy; we hope there are no errors in this present work, but.

In **[77]** we find a “very short proof” that eschews the Central Limit Theorem as claimed to be used in **[12]** because that “cannot reasonably be considered elementary.” Instead, the author (Patin) uses the Lebesgue dominated convergence theorem! In fact the authors of **[12]** criticize earlier papers, namely **[53]** and **[106]**, for inappropriately using the Central Limit Theorem which they say gives results valid only for the asymptotics of $n!$ for integers $n$. They use moment generating functions instead.

The paper **[79]** uses the Poisson distribution to get Stirling’s formula (including the $\sqrt{2\pi}$). The paper **[60]** uses change of variables. The pair of papers by Reinhard Michel **[64]** and **[65]** use artful estimates of $\ln$ and skilful substitution to get several terms of Stirling’s formula in a really elementary way. The second paper has our favourite title, “The $(n+1)$th proof of Stirling’s formula”, and takes a swipe at Patin for calling the Lebesgue dominated convergence theorem “elementary”. The letter to the Editor **[70]** takes a different swipe and claims Patin’s proof is contained in Newman’s Problem Seminar.

Six years after the $(n+1)$th proof, we have another “new proof” in **[72]**.

The very impressive paper **[28]** gives a concise but complete version of Laplace’s argument, and a nice historical discussion, including mention of de Moivre.

The paper **[66]** gives an elementary proof that

[TABLE]

and thereby extends the power of the root test.

The paper **[67]** did something quite different, namely use the Legendre duplication formula (6) and generalization to get some new recurrence relations for the logarithmic variant, equation (28). This generated the follow-up paper **[27]** because this gave new formulae for Bernoulli numbers.

The paper **[31]** extends Stirling’s formula for use with a certain generalization of central binomial coefficients. Its focus seems to be statistical, and it uses discrete random variables and the Central Limit Theorem.

The paper **[63]** is not without flaws, and indeed was criticized in **[45]** for not being as “new” as its title claimed666We found [45] by browsing the Monthly, not by citation search. Its criticism of [63] is largely sound: Copson’s treatment is better in many respects. , having been anticipated by the textbook treatment in **[21]**. Indeed, **[21]** points further back, to Watson and Ramanujan for the main idea. Nonetheless **[63]** is an interesting and readable paper, and the recurrence relation developed there, which we have already presented in equation (32), is apparently not in **[21]** and seems to be new to the paper. This paper was discussed in some detail in **[15]**.

One can say more, nowadays: in **[16]** and **[71]** we find “explicit” formulas in terms of 2-associated Stirling numbers, and David Jeffrey reports (personal communication) the following expression:

[TABLE]

Here $\begin{bmatrix}m\\ n\end{bmatrix}_{\geq 2}$ means a $2$-associated Stirling cycle number. It is gratifying to see Stirling numbers in Stirling’s asymptotic series. [Even if it’s really due to de Moivre].

Gaps: midpoint rule, Stirling’s original formula, and asymptotics of $\mathrm{inv}\Gamma$

**

We saw in a previous section that several papers used the trapezoidal rule to estimate $\ln n!$ by $\int^{n}_{1}\ln t\ dt$; we used this to give a parameterized upper bound, namely equation (39).

Nowhere in the Monthly that we have seen does the midpoint rule appear to be used to give a lower bound for $n!$. The following theorem rectifies this omission.

Theorem 1**.**

If $1\leq k\leq n-1$, and

[TABLE]

then

[TABLE]

Proof.

For a concave function, the midpoint rule gives an upper bound777This is in several texts, including [22]— but one place to learn it is by teaching from [93]., so if $j\geq 1$

[TABLE]

Using panels of width 1, if $k<n$,

[TABLE]

Adding $\sum_{j=1}^{k}\ln j$ to both sides,

[TABLE]

or

[TABLE]

Exponentiating gives the theorem. ∎

Remark**.**

We now consider what happens when $k=n$ or $m=n$: the inequalities become equality, and we learn $n!=n!$ which is not a surprise. But what if $k>n$, or $m>n$? In those cases the sense of the inequalities reverse because the integrals change sign:

[TABLE]

if both $m>n$ and $k>n$. One can establish also that $d_{m}$ decreases and is bounded below by $\sqrt{2\pi}$, while $c_{k}$ increases and is bounded above by $\sqrt{2\pi}$; we have

[TABLE]

These bounds and the finitary bounds of the theorem, i.e.

[TABLE]

if $k,m<n,$ are similarly tight to those produced by examining the error term in the trapezoidal rule (e.g. [48]).

Remark**.**

The first draft of this paper contained the phrase “we believe this is new, although it’s difficult to be sure.” Jon read that and was skeptical. It turns out that he was correct in part, but the story is interesting. Our first real clue was the paper [30] which gives, at the end, an approximation containing the characteristic factor $(n-1/2)^{n-1/2}$ (for $\Gamma(n)$, not $n!$). They attribute this result to Burnside [18], in $1917$ (a hundred years ago!). But the form is much older than that, even: Stirling himself gives (using “$\ell$”, to denote the base–10 logarithm)

[TABLE]

“added to half the logarithm of the circumference of the circle whose radius is one” in example 2 of Proposition 28 **[97, p. 151]**. This is an expression for $\log_{10}(n!)$ but here $\boldsymbol{z=n+1/2}$; the constant $a$ is $1/\ln(10)$ and this means, in modern notation,

[TABLE]

This is not equation (28)! This is, as Tweddle points out in **[96]**, a mid-point formula! Stirling did not explicitly use the midpoint rule or the concavity of logarithm to establish an inequality; instead he used telescoping series, induction, and a (clear and algorithmic) recipe for computing more terms to give his series. Apparently, de Moivre was the first to give the simpler version that we call “Stirling’s formula” today, namely equation (28).

The following theorem can be proved using very similar steps as equation (40) is proved in **[103]** and **[88]** (and is rather fun to do). We leave this for the reader.

Theorem 2**.**

For $\operatorname{Re}z>1/2$,

[TABLE]

This should be contrasted to equation (40); the known series expansion for $\operatorname{csch}(z)$ (see e.g. **[1]**) can then be used to fully develop Stirling’s original formula, after an appeal to Watson’s lemma **[9]**:

[TABLE]

Burnside **[18]** gave, by entirely different means, the formula

[TABLE]

which, being convergent, gives therefore a double sum for this Binet-like integral when $-1+z=N$ is an integer. Burnside did not comment on the similarity to Stirling’s original formula. He did, however, compute a constant

[TABLE]

as “to six figures” as $0.418938$, and then identified it as $\ln\sqrt{2\pi/\mathrm{e}}=0.418938533\ldots$ using “Stirling’s Theorem”, but using $N^{N+1/2}$ not $(N+1/2)^{N+1/2}$. This suggests that he did not realize that his formula was actually closer to Stirling’s original.

Asymptotics of inv$\Gamma$

The first few terms of Stirling’s original formula give

[TABLE]

We wish to solve $x=\Gamma(y)$ for $y$ when $x$ is large.

Put $v=x/\sqrt{2\pi}$ and $u=y-1/2$. Then the formula above gives

[TABLE]

Ignoring the $\mathcal{O}(1/u)$ term, we have

[TABLE]

and, somewhat remarkably, we can solve this for $u$, as follows:

[TABLE]

so (using the Lambert $W$ function)

[TABLE]

Exponentiating,

[TABLE]

This formula approximately inverts $\Gamma$, so $x\doteq\Gamma(y_{0})$. Using James Stirling’s own example, $1000!$ or $\Gamma(1001)\doteq 4.0238\cdot 10^{23567}$, the formula above gives $1000.99999397$, a relative error of about $6\cdot 10^{-9}$, from the computed value of $1000!$.

*The branch of $W$ used above is the principal branch, $W_{0}(x)$, defined for $-1/\mathrm{e}\leq x$ and satisfying $W_{0}(x)\geq-1$. Thus $\ln(x/\sqrt{2\pi})/\mathrm{e}\geq-1/\mathrm{e}$, or $x\geq\sqrt{2\pi}/\mathrm{e}$. *

At $x=\sqrt{2\pi}/\mathrm{e}$, $y_{0}=1/2-1/W_{0}(-1/\mathrm{e})=1/2+1=3/2$. This states that $x=\Gamma(y_{0})$ should be approximately true when $x=\sqrt{2\pi}/\mathrm{e}$ and $y_{0}=3/2$. We can check this: $\Gamma(3/2)=\Gamma(1/2)/2=\sqrt{\pi}/2$. Now $\sqrt{2\pi}/\mathrm{e}\doteq 0.922137$, $\sqrt{\pi}/2\doteq 0.88622$. These differ only by about 4%. Note that this accuracy is attained at the lowest possible value of $x$ for which the formula is valid at all, and that (as we shall see) the error decays reasonably rapidly as $x\to\infty$. This is very nearly a universal formula. Of course this is because Stirling’s formula is so accurate888We remind the reader that divergent series are often remarkably accurate. There are two relevant limits, only one of which might be bad. However, Stirling’s formula is famously accurate, more so than most. And Stirling’s original formula is even better, because the midpoint rule is twice as accurate..

But it gets better. For $-1/\mathrm{e}\leq\ln(x/\sqrt{2\pi})/\mathrm{e}<0$, we may use $W_{-1}(x)$ as well. That is, on

[TABLE]

the formula

[TABLE]

also inverts Stirling’s original formula. Plotting, on the same graph as $\mathrm{inv}\Gamma$, we get the picture in figure 7. The agreement is astonishing.

We claim but do not prove here (this paper is already too long) that the relative error in $y$ is $\mathcal{O}(1/y^{2})$ as $x\to\infty$; more precisely, the next two terms in the series are

[TABLE]

where $u_{0}=\ln(x/\sqrt{2\pi})/w$ with $w=W_{0}(\ln(x/\sqrt{2\pi})/\mathrm{e})$. There is a surprising connection to one of the series in **[98]**, but this will be discussed another day.

We remark that de Moivre’s simpler formula

[TABLE]

is not as simple to invert:

[TABLE]

doesn’t have a solution in terms of Lambert $W$, and thus we’re forced to approximate:

[TABLE]

and ignoring the $-1/2\ln y$, we get the same equation as before, but now

[TABLE]

which always differs by $1/2$ from the previous.

Drawing the line somewhere

We could have included here a look at **[5]**, **[34]**, **[39*]**, *

* [40], [41], [56], [75], [105], or for that matter [55], but no. Read them yourselves!*

4 Concluding Remarks

This has been an amazing journey for us, through the Monthly and through a selection of the vast $\Gamma$ literature. The $\Gamma$ bibliography at milanmerkle.com had $986$ entries when we accessed it. This present survey has focussed on Monthly papers, with some others for context. Diego Dominici’s survey of Stirling’s formula **[29]** attempted some of this, and we found his paper (in the Nota di Matematica) helpful; but indeed this whole project was a wonderful excuse to read great mathematical works. The ostensibly novel contributions of this paper (the pictures of the Riemann surface for inv$\Gamma$, the asymptotic formula for inv$\Gamma$, and the elementary bracketing inequality (33) for $n!$ that has approximations to Stirling’s original formula on the left and de Moivre’s simpler formula that we call “Stirling’s formula” on the right), are small in comparison.

Stirling’s original series for $n!$ is largely unknown, although the historians have it right; and Wikipedia correctly attributes the series popularly known as “Stirling’s formula” to de Moivre. For more information on de Moivre’s contributions, see **[8]**. Recently, Peter Luschny **[61]** and Weiping Wang **[98]** have started calling the exponentiated version of Stirling’s original formula by the name “the de Moivre formula” **[98, pg. 584]**! In view of similar historical misnamings (Euler invented Newton’s method; Newton invented (symplectic) Euler’s method, for instance **[99]**), we’re ok with this.

After writing this paper, we discovered that David W. Cantrell in 2001 **[19]** found the asymptotic formula for $\mathrm{inv}\Gamma$, before us. However, we then found that it’s also in **[42]**, beating us all by decades!

Dedicated to the memory of Jonathan M. Borwein.

About the Paper

I proposed the idea of this paper to Jon in late May 2016, and he was immediately enthusiastic. Jon suggested that we start by scouring JSTOR for Monthly titles containing “Gamma” or “Stirling’s formula”. Later, we added the terms “factorial” and “psi” to the search. I wrote the first draft of the paper, containing the elements of every section except section 2, meaning to leave that to Jon. In our many lunchtime discussions of the evolving draft over the months of June and July, Jon made many contributions, including pointing out connections to some of his early work. He also suggested that one of my results in section 3 was not as ‘new’ as I had thought; as I now know, Jon was right.

Jon died early in the morning August 2, 2016, the week I was to pass the draft over to him. Jon’s influence on this paper is therefore much less than it could have been. Now we’ll just have to imagine what Jon would have done with section 2.

On a final note, the task that this paper attempts is impossible in reasonable time (and we both knew it). A sample of the Problems and Solutions section turns up $\Gamma(z)$ or $n!$ or $\psi$ more frequently than not, as the simplest rebuttal of a claim to completeness. But I like to think that this paper achieves 80$\%$* coverage. Had Jon lived, we might have achieved 90*$\%$.

Judith Borwein has graciously permitted Jonathan’s name to be listed as author; certainly I feel it’s justified. It wasn’t just details: part of this paper was his vision.

• ACKNOWLEDGMENT.

This work was supported by Western University through a Distinguished Visiting Research Fellowship for the first author, and by NSERC, the Fields Institute for Research in the Mathematical sciences, the Rotman Institute for Philosophy, and the Ontario Research Centre for Computer Algebra for the second. The manuscript was typed by Julia E. Jankowski and Chris Brimacombe.