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 '20 at 13:35
• Purely by observation and not as accurate as the approximation in the question (but maybe a useful starting point of sorts?), $\frac{1}{\Gamma(x+1)}-1$ on $[0,1]$ is approximately $\frac{\sin(\pi(x+2))}{\pi(x+2)}$. I would think I am not the first to observe this. Sep 2 '20 at 5:30
• @njuffa. This is a nice and interesting obervation. I shall look around. Thanks. Sep 2 '20 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 '20 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 '20 at 12:43

$$\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 '20 at 7:51

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.

This answer is curve fitting, but it involves only 3 parameters. The accuracy is better than 0.13% over all 0<x<1; 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 '20 at 4:41

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 '20 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 '20 at 6:29
• @gary. No idea ? Cheers :-) Sep 12 '20 at 9:48

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 :

• Thanks, Erik ! Did you read my answer ? Sep 6 at 9:50
• @ClaudeLeibovici Yes and I can finish with the closed interval method ! Sep 6 at 9:53