Eulerian integrals and generalized factorials

Anybody with one-to-two years of basic calculus exposure can easily work out the following definite integral of $x^n$: $$\int_{0}^k x^n \,dx = \frac{k^{n+1}}{n+1} = {\rm L}(n + 1)$$ where ${\rm L}(n) = {k^n}/{n}$. This is a mid-17th-century result known to both Newton and Leibniz, but their results were not written in modern notationIn his 1669 paper “De Analysi per Aequationes Numero Terminorum Infinitas" (Of analysis by equations of an infinite number of terms), Newton wrote: Let the base $AB$ of any curve $AD$ have $BD$ for its perpendicular ordinate; and call $AB = x$, $BD = y$, and let $a, b, c,$ etc be given quantities, and $m$ and $n$ whole numbers, then following rule (Rule I) can be used to compute the quadrature of simple curves. If $ax^{m/n}=y$, then $$\frac{an}{m+n}x^{\frac{m+n}{n}} = {\rm area} \, ABD$$ Newton then proceeded to give a number of examples. He wrote: The thing will be evident by an example. If $x^2 = y$, that is $a = 1 = n$, and $m = 2$, it shall be $$\tfrac{1}{3}x^3={\rm area}\,ABD$$ Many years later, in his “A treatise of the method of fluxions and infinite series: With its application to the geometry of curve lines" (published in 1737), Newton used the following notation. If the order of curve is $dz^{n-1}=y$, then the value of the area is $$\tfrac{d}{n}z^n = t$$ (see the first example in: A table of some curves related to rectilinear figures, constructed by Problem VII, p. 137.) .

We then fast-forward to 18th century St. Petersburg to see how this result can be extended if you put $\frac{1}{e^{x}}$ inside the box. $$\int_{0}^k x^n \, \square \,dx = \underbrace{ \bigcirc = \Gamma(n+1)}_{\textrm{when $k \to \infty$ }}$$ and the thing that you should put in the circle is actually our factorial function $$\bigcirc = n!$$

To illustrate this, we suppose $f(n, k)$ is a function defined by $$f(n, k) = \int_0^k \frac{x^n}{ e^{x}}\,dx$$ Then using integration by parts ($v = -e^{-x}$ is used as the first seed) $$\int_0^k u v'\,dx = \left. uv \right|_{0}^k - \int_0^k u' v\,dx$$ it is easy to show that $$f(n, k) = n!- n! \cdot g(n, k)$$ for the first few integer values of $n$, where $$\begin{align} g(1, k) &= & \frac{1+k}{e^k}\\ g(2, k) &= & \frac{1+k+\frac{1}{2!}k^2}{e^k}\\ g(3, k) &= & \frac{1+k+\frac{1}{2!}k^2+\frac{1}{3!}k^3}{e^k}\\ g(4, k) &= & \frac{1+k+\frac{1}{2!}k^2+\frac{1}{3!}k^3+\frac{1}{4!}k^4}{e^k}\\ \end{align}$$ and so on. In general we can write: $$g(n,k) = \frac{1 + k + \frac{1}{2!}k^2 + \frac{1}{3!}k^3 + \ldots + \frac{1}{n!}k^n}{1 + k + \frac{1}{2!}k^2 + \frac{1}{3!}k^3 + \ldots + \frac{1}{n!}k^n + \textrm{ad infinitum}}$$ When $k \to \infty$, it is easy to see that $g(n, \infty) = 0$ and it is very tempting to write: $$f(n, \infty) = n!$$ In other words, it says that the value of $n!$ can be calculated from an integral, that is, $$n! = \int_0^\infty \frac{x^n}{ e^{x}}\,dx$$

But I suppose Euler didn't derive his initial results in factorial function through this route in the late 1720s Euler number $e$ was discovered by Jacob Bernoulli in 1683 via a gedankenexperiment with bank interest, and Euler's father is a friend to the Bernoulli family since he used the following integral $$\int_0^1 (-\ln u)^n\,du$$ in his letterThis is actually Euler's second letter to Goldbach. The first letter was dated October 13, 1729. In the first letter, the following formula: $$m! = \frac{1\cdot 2^m}{1+m}\frac{2^{1-m}\cdot3^m}{2+m}\frac{3^{1-m}\cdot4^m}{3+m}\frac{4^{1-m}\cdot5^m}{4+m}\cdots$$ was announced. Actually, Goldbach received another from letter (in French, dated October 6, 1729) from Daniel Bernoulli on the same topic. Bernoulli claimed that when $A$ approaches infinity, we must have: $$x! = \left(A + \frac{x}{2}\right)^{x-1} \frac{2}{1+x}\frac{3}{2+x}\frac{4}{3+x}\cdots\frac{A}{A-1+x}$$ Note that Bernoullis's infinite product looks very different from that of Euler's. Bernoulli gave two numerical examples in his letter. He showed by taking $x = \frac{3}{2}, A = 8$, he was able to get $\frac{3}{2}! \approx 1.3005$. (Actually Bernoulli miscalculated $\frac{3}{2}!$ since $\sqrt{\frac{35}{4}}\frac{4}{5}\frac{6}{7}\frac{8}{9}\frac{10}{11}\frac{12}{13}\frac{14}{15}\frac{16}{17}\approx 1.32907$, he must have mis-substituted $3$ for $\frac{3}{2}$ in the surd expression since he wrote $\sqrt{\frac{19}{2}}$ instead of $\sqrt{\frac{35}{4}}$) And when $x = 3, A = 16$, he obtained $3! \approx 6\frac{1}{204}$. to Christian Goldbach on January 8, 1730. Later, Euler must have realized that his logarithmic definition equivalent to the exponential definition if we let $x = -\ln u = \ln \frac{1}{u}$The following notation was used by Euler $$y = \int_0^1 du \left( \ln \tfrac{1}{u}\right)^p \quad {\rm and} \quad y = \int_0^\infty e^{-v}v^p\,dv$$ when he presented his paper to St. Petersburg Academy on December 19, 1765. See Euler (1765) De curva hypergeometrica hac aequatione expressa $y = 1\cdot 2\cdot 3 \cdots x$, Novi Commentarii academiae scientiarum Petropolitanae, Volume 13, pp. 3-66.. This logarithmic integral is often called Eulerian integral of the first kind.

The integral symbol ∫ introduced in 1675 was based on the ſ (long S) character, and was chosen because Leibniz thought of it as an ∞ sum pic.twitter.com/ZFr3WFKbDF

— Fermat's Library (@fermatslibrary) August 26, 2017

Since the logarithmic integral predates the exponential form. $$\int_0^1 (-\ln x)^n\,dx \quad {\rm predates} \quad \int_0^\infty \frac{x^n}{e^x}\,dx$$ The exponential form of the integral $$\int_0^\infty \frac{x^n}{e^x}\,dx$$ is often called the Eulerian integral of the second kind in the math literature.

Eulerian integral is important because it gives us a way to compute non-integer value or even negative values of $n$. For example, the factorial of $\frac{1}{2}$ is given by: $$\left(\tfrac{1}{2}\right)! = \int_0^\infty \frac{x^{1/2}}{ e^{x}}\,dx = \tfrac{1}{2}\sqrt{\pi} = 0.886227\ldots$$ To show that this is the case, we put $x = s^2$ in the integral, and write: $$\begin{align} \left(\tfrac{1}{2}\right)! &= \int_0^\infty s e^{-s^2}2s\,ds = \int_0^\infty s (2s e^{-s^2})ds\\ &= \left[s (-{e^{-s^2}})\right]_0^\infty -\int_0^\infty (-e^{-s^2})ds = \int_0^\infty e^{-s^2}ds = \tfrac{1}{2}\sqrt{\pi} \end{align}$$ The standard trick to evaluate the gaussian integral is to employ the transformation $s = r \cos \theta$, $t = r\sin \theta$ $$\int_0^\infty \int_0^\infty e^{-s^2} e^{-t^2} ds\,dt = \int_0^{\pi/2}\int_0^\infty e^{-r^2}r\,dr\,d\theta = \tfrac{1}{4}\pi$$ This double-integral trick was, according to Sturm in his is Cours d’Analyse de l’ecole polytechnique, first devised by Poisson.

The double-integral trick in Sturm's Cours d’Analyse, Volume II (1877 edition), the book was published posthumuously in 1859.

Around the same time in London, the Frenchman de Moivre was busying himself with probability distributions, and he needed to handle large factorials in his calculations, e.g. related to binomial coefficients like $\binom{2n}{n} = \frac{(2n)!}{(n!)^2}$. de Moivre solved his problems with a constant which he termed $B$, as in: $$\int_0^\infty \frac{x^n}{ e^{x}}\,dx \approx B\sqrt{n}\left(\frac{n}{e}\right)^n$$ where $B \approx 2.5074$In the second supplementary fascicle (titled: Approximatio ad Summam Terminorum Binomii $(a + b)^n$ in Seriem expansi, dated: November 12, 1733) to his Miscellanea Analytica (first published in 1730), de Moivre was able to obtain a infinite series for $\ln B$: $$\ln B = 1 - \tfrac{1}{12} + \tfrac{1}{360} - \tfrac{1}{1260} + \tfrac{1}{1680} +\ldots$$ The gap for the closed-form expression for $B$ was filled later by the Scottish mathematician James Stirling (1692 – 1770), who was a friend of Newton..

Common logarithms of factorials ($\log_{10} 10!$ to $\log_{10} 900!$) in the 1733 edition of Miscellanea Analytica (pp. 103-104). Note that the figures past the fifth decimal are incorrect (the numbers were eventually corrected when de Moivre published his Doctrine of Chances) The problem was highlighted to de Moivre by Stirling in a letter in 1730 (when the first edition of the book is released). Stirling told de Moivre in the letter that he knew the value of $B$ is exactly $\sqrt{2\pi}$.

We would like to show that the de Moivre's constant is actually $$B = \sqrt{8}\left(\tfrac{1}{2}\right)! = \sqrt{2\pi}$$ Now suppose $y = \frac{x^n}{ e^{x}}$, then we can easily show that $y$ is at its maximum when $x = n$, that is $y_{\rm max} = \left(\frac{n}{e}\right)^n$. Suppose that we define $y$ in term of its maximum value as $y = \left(\frac{n}{e}\right)^n e^{-t^2}$, that is, $$\int_0^\infty x^n e^{-x}\,dx = \int_0^\infty \left(\frac{n}{e}\right)^n e^{-t^2}\,dx = \left(\frac{n}{e}\right)^n\int_0^\infty e^{-t^2}\,\color{brown}{dx}$$ The key to the evaluation of this integral is that the element $\color{brown}{dx}$ must be re-expressed in term of $t$. To do this, we suppose $x$ can be expanded into a power series $x = n + \theta(t) =n + \sum_{k=1}^\infty a_k t^k$, from which the following relation must be satisified: $$(n+\theta)^n e^{-(n + \theta)} = \left(\frac{n}{e}\right)^n e^{-t^2}$$ from which $t$ can be readily expressed a function of $\theta$: $$t =\sqrt{\theta - n \ln \left(1 + \tfrac{\theta}{n}\right)} = \tfrac{1}{\sqrt{2n}}\left(\theta - \tfrac{1}{3n}\theta^2 + \tfrac{7}{36n^2}\theta^3 - \tfrac{73}{540n^3}\theta^4 +\ldots\right)$$ which can then be rewritten as: $$\begin{align} \tfrac{t}{1/\sqrt{2n}} &\equiv \left(\sum_{k=1}^\infty a_k t^k\right) - \tfrac{1}{3n}\left(\sum_{k=1}^\infty a_k t^k\right)^2\\ & \qquad + \tfrac{7}{36n^2}\left(\sum_{k=1}^\infty a_k t^k\right)^3 - \tfrac{73}{540n^3}\left(\sum_{k=1}^\infty a_k t^k\right)^4 +\ldots\\ &\equiv b_1 t + b_2 t^2 + b_3 t^3 + b_4 t^4 + \ldots \end{align}$$ where $b_k = 0, k > 1$. Comparing the value of coefficients $b_k$ gives: $$\begin{align} \sqrt{2n} &= a_1\\ 0 &= -\tfrac{1}{3n}a_1^2 + a_2\\ 0 &= \tfrac{7}{36n^2}a_1^3 - \tfrac{2}{3n}a_1 a_2 + a_3\\ 0 &= -\tfrac{73}{540n^3}a_1^4 + \tfrac{7}{12n^2}a_1^2 a_2 - \tfrac{1}{3n}(a_2^2 + 2a_1 a_3) + a_4\\ 0 &= -\tfrac{73}{135n^3}a_1^3 a_2 + \tfrac{7}{12n^2}a_1 (a_2^2 + a_1 a_3) - \tfrac{2}{3n}(a_2 a_3 + a_1 a_4)\\ \vdots &= \vdots \end{align}$$ Therefore, we can obtain an expression for $dx = d\theta$: $$\begin{align} \theta &= t\sqrt{2n} + \tfrac{2}{3}t^2 + \tfrac{1}{9\sqrt{2n}}t^3 - \tfrac{2}{135n}t^4 + \tfrac{1}{540n\sqrt{2n}}t^5 + \ldots\\ d\theta &= dt \left[\sqrt{2n} + \tfrac{4}{3}t + \tfrac{1}{3\sqrt{2n}}t^2 - \tfrac{1}{270n}t^3 + \tfrac{1}{108n\sqrt{2n}}t^4 + \ldots\right] \end{align}$$ where $t \in (-\infty, \infty)$ for $x \in [0, \infty)$.

Graphical representation of $3!$: The area under the blue curve is given by: $$3!=\int_0^\infty u^3 e^{-u}\,du$$

The area under the orange curve is given by:

$$3! = \left(\tfrac{3}{e}\right)^3\int_0^\infty e^{-u^2}dw = \left(\tfrac{3}{e}\right)^3\int_{-\infty}^\infty e^{-u^2}du\left(\sqrt{6}+\tfrac{4u}{3}+\tfrac{u^2}{3\sqrt{6}}-\tfrac{u^3}{810}+\tfrac{u^4}{324\sqrt{6}}+\ldots\right)$$

$n!$ can now be rewritten as: $$\begin{align} n! &= \left(\frac{n}{e}\right)^n\int_{-\infty}^\infty e^{-t^2}dt \left[\sqrt{2n} + \tfrac{4}{3}t + \tfrac{1}{3\sqrt{2n}}t^2 - \tfrac{1}{270n}t^3 + \tfrac{1}{108n\sqrt{2n}}t^4 +\ldots\right]\\ &= \left(\frac{n}{e}\right)^n \sqrt{2n\pi} \left(1 + \tfrac{1}{12n} + \tfrac{1}{288n^2} + \ldots \right) \end{align}$$ The series on the right hand side of the equation actually has a name. It is called Stirling series: $$1 + \tfrac{1}{12n} + \tfrac{1}{288n^2} -\tfrac{139}{51840n^3} - \tfrac{571}{2488320n^4} + \tfrac{163879}{209018880n^5} + \tfrac{5246819}{75246796800n^6}+ \ldots$$

In general, we can use the Laplace method to handle asymptotic expansion associated with function of the following form: $$I_n = \int_a^b g(x) e^{-n \,f(x)}\,dx \sim g(x_0) e^{-n \,f(x_0)} \sqrt{\frac{2\pi}{n f''(x_0)}}$$ where $f'(x_0) = 0$ and $f''(x_0) > 0$.

The Laplace transform method of solving linear differential equations in control engineering. For example, if we let $x = st$, then $$\int_0^\infty (st)^n e^{-st}\,d(st) = \int_0^\infty s^{n+1}\,t^n e^{-st}\,dt=s^{n+1}\int_0^\infty t^n e^{-st}\,dt$$ that is, $$\mathcal{L}[t^n] = \int_0^\infty t^n e^{-st}\,dt = \frac{\Gamma(n+1)}{s^{n+1}}$$ obviously is byproduct of Laplace's exploration of this type of integral.