# Approximation of $\Big[\Gamma(1+x)\Big]^{-1}$ for $0 \leq x \leq 1$ (for the art for art's sake).

Recently was asked on the site the question : how to compute $$f(t)=\int_0^t\frac{\mu^x}{\Gamma(x+1)}\ dx$$ for which I assumed $$0 \leq t \leq 1$$. There is no antiderivative available from any CAS I tried.

My first idea was to approximate $$\Gamma(x+1)$$ by a more or less constrained polynomial in $$x$$, then to use partial fraction decomposition. This leads to a linear combination of elliptic integrals.

Later, I found that $$\Big[\Gamma(1+x)\Big]^{-1}$$ looks to be very similar to a Gibbs excess energy diagram for a binary system. Some models of this physical properties are very simple (Van Laar, Margules), but they only have two parameters which is not enough for an "accurate" representation of the function. Much more accurate are Scatchard-Hildebrand, Wilson, NRTL or Uniquac models but their complexity would not allow the required integration.

So, with a limited choice, I decided to try with a model looking like the one proposed by Redlich-Kister which contains a pure polynomial component. So my idea was to write $$\frac{1}{\Gamma(x+1)}\sim 1+x(x-1) \sum_{k=0}^p d_k\, x^k$$ $$f(t)=\frac{\mu ^t-1}{L}+$$ $$\sum_{k=0}^p (-1)^k\,d_k L^{-(k+3)} (L\, \Gamma (k+2,-L t)+\Gamma (k+3,-L t)-(k+L+2) \Gamma (k+2))$$ where $$L=\log(\mu)$$ (provided that $$\Re(\log (\mu ))<0\land \Re(k)>-2$$).

Willing to stay non-empirical, in my answer, I used $$p=3$$ and the $$d_k$$'s were computed in order to match the function and first derivative values at $$x=0$$, $$x=\frac 12$$ and $$x=1$$. The result were no too bad.

For the art for art's sake, I decided to use $$p=6$$, the $$d_k$$'s were computed in order to match the function, first and second derivative values at $$x=0$$, $$x=\frac 12$$ and $$x=1$$. This leads to $$d_0=-\gamma \qquad d_1=-\gamma -\frac{\gamma ^2}{2}+\frac{\pi ^2}{12}$$

$$4 \sqrt \pi\,d_2=6 (178+3 \gamma (8+\gamma )) \sqrt{\pi }+32 \pi ^2-3 \pi ^{5/2}-64 \left(P^2+4 P+36\right)$$ $$3 \sqrt \pi\,d_3=-3 (1356+\gamma (111+16 \gamma )) \sqrt{\pi }-144 \pi ^2+8 \pi ^{5/2}+96 \left(3 P^2+8 P+92\right)$$ $$3 \sqrt \pi\,d_4=2 \left(3 (1376+\gamma (61+14 \gamma )) \sqrt{\pi }+156 \pi ^2-7 \pi ^{5/2}-24 \left(13 P^2+20 P+372\right)\right)$$ $$\frac{\sqrt{\pi }}{4}\,d_5=-\left(628+9 \gamma +6 \gamma ^2\right) \sqrt{\pi }-24 \pi ^2+\pi ^{5/2}+16 \left(3 P^2+2 P+84\right)$$ $$3 \sqrt \pi\,d_6=-4 \left(-6 (106+(\gamma -1) \gamma ) \sqrt{\pi }-24 \pi ^2+\pi ^{5/2}+48 \left(P^2+28\right)\right)$$ where $$P=\psi \left(\frac{3}{2}\right)=2-\gamma -2\log (2)$$

The maximum absolute error is $$1.5 \times 10^{-8}$$ which seems to be decent.

Making the numbers rational, this would give $$\frac{1}{\Gamma(x+1)}\sim 1+x(1-x) P_6(x)$$ $$P_6(x)=\frac{2807}{4863}-\frac{247}{3140}x-\frac{461 }{3820}x^2+\frac{66 }{1435}x^3+\frac{11 }{3303}x^4-\frac{15 }{2726}x^5+\frac{3 }{2750}x^6$$

A few values of the definite integral (just for comparison) $$\left( \begin{array}{cccc} \mu & t & \text{approximation} & \text{exact} \\ 0.25 & 0.2 & 0.182857067200 & 0.182857068268 \\ 0.25 & 0.4 & 0.329965795201 & 0.329965797571 \\ 0.25 & 0.6 & 0.443007596935 & 0.443007599310 \\ 0.25 & 0.8 & 0.526656210875 & 0.526656212492 \\ 0.25 & 1.0 & 0.586607844050 & 0.586607845209 \\ & & & \\ 0.50 & 0.2 & 0.195704099746 & 0.195704100926 \\ 0.50 & 0.4 & 0.376453536397 & 0.376453539149 \\ 0.50 & 0.6 & 0.535922898037 & 0.535922900791 \\ 0.50 & 0.8 & 0.671420159082 & 0.671420160584 \\ 0.50 & 1.0 & 0.782934567076 & 0.782934567750 \\ & & & \\ 0.75 & 0.2 & 0.203784473259 & 0.203784474510 \\ 0.75 & 0.4 & 0.407825165336 & 0.407825168343 \\ 0.75 & 0.6 & 0.602995831714 & 0.602995834721 \\ 0.75 & 0.8 & 0.782793915893 & 0.782793917219 \\ 0.75 & 1.0 & 0.943235581611 & 0.943235581767 \end{array} \right)$$

I am more than happy with these results but,again, just for the art for art's sake, I would like to be even better. For sure, I could use $$p=9$$ and obtain the parameters in order to match the function, first, second and third derivative values at $$x=0$$, $$x=\frac 12$$ and $$x=1$$; this would lead to monstreous expressions.

So, my question is : without any curve fit and avoiding so many terms, is there a way to obtain a better (mathematically based) approximation of $$\Big[\Gamma(1+x)\Big]^{-1}$$ for $$0 \leq x \leq 1$$

• See my own question math.stackexchange.com/questions/3709724/… Sep 1, 2020 at 13:35
• @njuffa. This is a nice and interesting obervation. I shall look around. Thanks. Sep 2, 2020 at 7:30
• Claude have you considered the function $f(x)=ax\ln(x)+1$ with $0<a<1$ to make an approximation of the $\Gamma(x+1)$ with $0\leq x\leq 1$ ?Good day ! Sep 5, 2020 at 12:29
• @c-love-garlic. Good idea but bad result ! This would give $a\sim -0.34$. Any idea or suggestion is welcomed. Cheers :-) Sep 5, 2020 at 12:43
• @ClaudeLeibovici have you considered $f(x)=\Gamma(x^x+1)$ the abscissa of the minimum is near from $e^{-e^{-1}}$ ? Oct 26, 2020 at 17:30

$$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$$ \begin{align} &\bbox[5px,#ffd]{1 \over \Gamma\pars{x + 1}} = {1 \over x!} = {x + 1 \over \pars{x + 1}!} = {\pars{x + 2}\pars{x + 1} \over \pars{x + 2}!} \\[5mm] = &\ \cdots = {\prod_{k = 1}^{n}\pars{x + k} \over \pars{x + n}!} \\[5mm] \approx &\ {\prod_{k = 1}^{n}\pars{x + k} \over \root{2\pi}\pars{x + n}^{x + n + 1/2}\,\,\,\expo{-n - x}} = {\tt mySol}\pars{x,n} \end{align}

• Very interesting approach ! Thanks for it. By the way, when I see some answers, I know they are from you since you have a very personal style. I suppose that you don't remember but we met in Caracas some 30+ years ago. Cheers :-) Nov 20, 2020 at 7:51

My try : we have on $$x\in[1,2]$$ :

$$1+4(x-1)(x-2)\left(1-\Gamma\left(\left(\frac{\left(1+\sqrt{2}\right)}{2}\right)^{2}\right)\right)\simeq \Gamma(x)$$

It's not really tight but it's simple in the form.

Hope it helps you Claude.

Edit :

Around $$x=2$$ we have :

$$\Gamma(x)\simeq ((1-x)\frac{1}{2}(\gamma-1)+1+\frac{1}{2}(\gamma-1))^x$$

Edit 2:

Using the inverse function of the Lambert's function and denotes by :

$$g(x)=\left(1+\frac{1}{2}xe^x\right)^{x+1}$$

We have $$1\leq x \leq 2$$:

$$\Gamma(x)\simeq (g(1-x))^x$$

We have on $$(-1,0)$$ the approximation :

$$(\operatorname{W}(-ex))^{(\operatorname{W}(-x))^{1.45}}\simeq (\Gamma(1-x))^{\frac{1}{1-x}}$$

For $$1\leq x\leq 2$$ it seems we have the inequality :

## $$f(x)=\left(\Gamma(x)\tanh(x)\right)^{\frac{2}{x(1+\tanh(e^{-1}x))}}$$

$$f(x)\leq \frac{f(1)-f(2)}{-1}(x-1)+f(1)$$

And the difference is less than $$1*10^{-3}$$

Last edit :

## $$f(x)=\left(\Gamma(x)\tanh(x)\right)^{\frac{2}{x(1+\tanh(0.25x))}}$$

For $$1\leq x\leq 2$$ it seems we have the inequality :

$$(f(x))''>0$$

So we can get using convexity an arbitrary accuracy following the lenght of the interval .

• Thanks, Erik ! Did you read my answer ? Sep 6, 2021 at 9:50
• @ClaudeLeibovici What do you think about my edit ? Sep 20, 2021 at 8:57
• This is interesting ! Thanks, Erik ! Sep 20, 2021 at 10:12
• @ClaudeLeibovici Last edit what do you think about ? Oct 3, 2021 at 17:01
• Interesting again ! Thanks, Erik. Oct 4, 2021 at 1:08

Letting $$x!=\Gamma(x+1)$$, the Gamma function may be defined through the following limit for an arbitrary $$\alpha$$ where $$x\approx0$$ and $$x\approx\alpha^2$$ for $$\alpha\approx1$$ are most accurate:

$$x!=\lim_{n\to\infty}\left(n+\alpha\right)^x\prod_{k=1}^n\left(1+\frac xk\right)^{-1}$$

which has simple inverse given by

$$x!^{-1}=\lim_{n\to\infty}\left(n+\alpha\right)^{-x}\prod_{k=1}^n\left(1+\frac xk\right)$$

In the given problem, this leads to an integrand of the form $$P(x)\phi^x$$ for polynomial $$P$$ and constant $$\phi$$, which can be evaluated elementarily.

It's also worth pointing out that refining the estimate by taking larger $$n$$ does not involve a lot of recalculation since it is just one additional term on the product. However, the convergence is not very fast.

Choosing something such as $$\alpha=0.9$$ should lead to a good balance of over- and under-estimates of the Gamma function.

What is below is exactly what I do not want to do. It is just posted to show that high order polynomials are highly significant.

What I did was to generate $$10001$$ equally spaced data points of $$f(x)=\frac{1}{(x-1) x}\left(\frac{1}{\Gamma (x+1)}-1\right)\qquad \qquad -\gamma \leq f(x) \leq -1+\gamma$$ using high precision.

Then, polynomial regressions; below is an example of the results for a polynomial of degree $$9$$ (it leads to maximum absolute errors $$\sim 10^{-11}$$).

$$\begin{array} \text{ } & \text{Estimate of } d_k & \sigma\text{ for }d_k \\ 1 & -0.5772156649 & 2.449\times 10^{-13} \\ x & +0.0786624078 & 1.399\times 10^{-11} \\ x^2 & +0.1206650090 & 2.600\times 10^{-10} \\ x^3 & -0.0458731951 & 2.221\times 10^{-9} \\ x^4 & -0.0036781095 & 1.028\times 10^{-8} \\ x^5 & +0.0059542170 & 2.790\times 10^{-8} \\ x^6 & -0.0012903855 & 4.563\times 10^{-8} \\ x^7 & -0.0000841706 & 4.421\times 10^{-8} \\ x^8 & +0.0000893923 & 2.336\times 10^{-8} \\ x^9 & -0.0000138358 & 5.183\times 10^{-9} \end{array}$$

Edit

In view of these results, I made the calculations for $$p=9$$ matching the function, first, second and third derivative values at $$x=0$$, $$x=\frac 12$$ and $$x=1$$; this effectively leads to monstrous expressions (I shall not type them but they are available); as one could expect, the same constants plus $$\zeta(3)$$ appear in their expressions. The maximum absolute error is $$1.1 \times 10^{-11}$$.

Now, the same table as before $$\left( \begin{array}{cccc} \mu & t & \text{approximation} & \text{exact} \\ 0.25 & 0.2 & 0.182857068266992 & 0.182857068267690 \\ 0.25 & 0.4 & 0.329965797569747 & 0.329965797571307 \\ 0.25 & 0.6 & 0.443007599308124 & 0.443007599309699 \\ 0.25 & 0.8 & 0.526656212489610 & 0.526656212491646 \\ 0.25 & 1.0 & 0.586607845206624 & 0.586607845208931 \\ & & & \\ 0.50 & 0.2 & 0.195704100925464 & 0.195704100926239 \\ 0.50 & 0.4 & 0.376453539146954 & 0.376453539148764 \\ 0.50 & 0.6 & 0.535922900789419 & 0.535922900791250 \\ 0.50 & 0.8 & 0.671420160581682 & 0.671420160584277 \\ 0.50 & 1.0 & 0.782934567746625 & 0.782934567749710 \\ & & & \\ 0.75 & 0.2 & 0.203784474509130 & 0.203784474509953 \\ 0.75 & 0.4 & 0.407825168340619 & 0.407825168342596 \\ 0.75 & 0.6 & 0.602995834719402 & 0.602995834721404 \\ 0.75 & 0.8 & 0.782793917216033 & 0.782793917219063 \\ 0.75 & 1.0 & 0.943235581763057 & 0.943235581766778 \end{array} \right)$$

• For comparison, coefficients of minimax $P_{9}$ (matches at $x=0$, $x=1$) with maximum absolute error $< 10^{-12}$: 5.7721566490547482e-1, -7.8662407402331200e-2, -1.2066501595952826e-1, 4.5873243335408945e-2, 3.6779360284088752e-3, -5.9538708486641407e-3, 1.2900033631901313e-3, 8.4380677149800901e-5, -8.9427608944409487e-5, 1.3828612244534447e-5 Sep 2, 2020 at 17:54
• @njuffa. Thanks ! This is nice but the problem with minimax is that you do not have any idea about the significance of resulting parameters. Sep 3, 2020 at 6:29
• @gary. No idea ? Cheers :-) Sep 12, 2020 at 9:48

Inspired by @FelixMarin's answer I found :

Let $$2/3 and define :

## $$f\left(x\right)=\sqrt{2\pi x}\left(\frac{x}{e}\right)^{x}e^{\frac{1}{12x+\left(e-2\right)x^{-x}\cdot2^{-x}-\left(e-2\right)x^{x^{x}}+1-C}}$$

Where $$C$$ is the Catalan's constant .

Then it seems we have :

$$0<\frac{f\left(x\right)}{f\left(1\right)}-x!<8*10^{-5}$$

This answer is curve fitting, but it involves only 3 parameters. The accuracy is better than 0.13% over all $$0; good, but nothing like Claude's 8 digits or so. The idea is that the local max is not exactly at 1/2, as implied the factor $$x(1-x)$$ in Claude's analysis.

$$\frac{1}{\Gamma(1+x)} \approx 1+ A \ x^b (1-x)^c$$ $$A=0.538 , \ b=0.971 ,\ c=1.095$$

The error is close to minimax; that is, the amplitudes of the difference curve are approximately equal.

• Thanks for answering ! But, as I wrote, I do not want any empirical curve fit. Next, how would you compute the integral $f(t)$ ? "My" maximum is not at $x=1/2$; for the rigorous function, the coordinates of the extremum are $(0.461632145,1.12917389)$ and for the approximation $(0.461632154,1.12917389)$. Cheers and thanks again for the effort. Sep 2, 2020 at 4:41

Hi Claude I find the following nice:

Let $$x\in(0,1)$$ and

$$f(x)=x!^{\frac{1}{x-1}},\quad g(x)=\frac{ f(0)-f(1) }{-1}\left(x-1\right)+f(1)+\sum_{n=1}^{\infty}\ln\left(1+\frac{x^{2n}\left(1-x\right)^{2n}}{\left(2\pi\right)^{n}}\right).$$

Then $$k=10$$

$$\left(g(x)\right)^{x-1}+k-\sum_{n=1}^{k}\exp\left(\frac{x^{\frac{4}{\pi}}\left(1-x\right)^{8+n}}{8^{2}}\right)>x!$$

• @Gary Thanks for the edit ! Jan 23, 2023 at 9:29
• We can also use : $$f(x)=x!^{\frac{b}{x-1}},b=\left(1+\frac{1}{x+25000}\right)^{x}$$ Jan 23, 2023 at 9:59
• @ClaudeLeibovici see the derivative of $$\frac{\left(x^{\ln\left(x-e^{-\left(x^{\ln\left(3\right)}\left(1-x^{\ln\left(3\right)}\right)\right)}+1-x^{\ln\left(3\right)}\left(1-x^{\ln\left(3\right)}\right)\right)}\right)}{\left(\frac{1}{x}\right)!}-x!$$ around $x=1$ Apr 27, 2023 at 10:50