3

Question in a nutshell:

Can anyone point me to an algorithm for computing to double-precision floating-point (roughly 16 digits) the inverse of either log gamma or log factorial?

In other words, if y = ln(factorial(x)), and I am given double-precision y, I want to compute double-precision x.

Details of what I have learned so far:

For some engineering work that I am doing, I have implemented both factorial(x) and lnfactorial(x) for non-negative real x using Lanczos Approximation. Now I need to implement the inverses of both. Inversing the Lanczos Approximation that I am using (a 13-element series expansion) is beyond my math-jitsu skills.

The best answer I have found is in a few posts by David W. Cantrell:

Link (Inverse gamma function)

http://mathforum.org/kb/message.jspa?messageID=342644 (23 kt. asymptotic expansion for gamma function)

https://groups.google.com/forum/#!topic/sci.math/EIAf7rQ6A_4 (A new convergent expansion for the gamma function)

I think those will get me the inverse factorial function that I need, once I implement a double-precision Lambert W function. I am working on that now.

My more detailed questions for you:

(1) Those posts are from 2001… what is the latest and greatest info on how to efficiently compute the inverse factorial/gamma?

(2) If that is the latest, how many coefficients are needed to give me 15 decimal digits of precision? And has anybody documented those n coefficients somewhere (given double-precision math is pretty common; I was able to find pre-computed coefficients for Lanczos that work well for double-precision)?

(3a) I am also trying to implement lnfactorial (natural log of the factorial) and its inverse… did Cantrell (or anybody else) develop a similar approximation for that? (3b) Or is anyone's math-jitsu skills strong enough to tell me how to use his inverse gamma approximation to compute inverse lngamma? (lnfactorial was computable with a very similar Lanczos approximation as factorial; so, I am hoping the same will be true of the inverses.)

Any guidance/leads you can give on any of that would be greatly appreciated!

Thanks,
    Brian
Glorfindel
  • 4,000
Brian Kennedy
  • 235

1 Answers1

3

Concerning the inverse of the gamma function, I would strongly suggest tha approximation proposed on this site by @robjohn in $2015$ (have a look to this post which corresponds to my latest answer on this topic and go backward to the linked pages).

For general $n$, the formula is $$n\sim e\exp\left(\operatorname{W}\left(\frac1{e}\log\left(\frac{\Gamma(n)}{\sqrt{2\pi}}\right)\right)\right)+\frac12$$ that you could write in a different way taking into account the fact that $z=W(z)\,e^{W(z)}$.

If you look for integer solution, then $$\color{blue}{n=\left\lceil e\exp\left(\operatorname{W}\left(\frac1{e}\log\left(\frac{\Gamma(n)}{\sqrt{2\pi}}\right)\right)\right)+\frac12 \right\rceil}$$ is stricly perfect.

Concerning the approximation of the factorial, I proposed here to use
$$n!\sim\sqrt{\pi}\,\left(\frac{n}{e}\right)^n\root\LARGE{6}\of{8n^3+4n^2+n+x(n)}$$ where $$x(n)=\frac{1}{30}+\frac{1455925-70798882\, n}{1544702880\, n^2+760646880\, n+981836548}$$ which I could still improve if this is of interest for you.

Edit

As I wrote above, it is possible to improve the calculation of $n!$ making $x(n)$ as a $[k,k+1]$ Padé approximant. For example, using $k=4$, $$\left( \begin{array}{ccc} 1 & 1 & 0.999998383049203632780 \\ 2 & 2 & 1.999999991752156826810 \\ 3 & 6 & 5.999999999668965248865 \\ 4 & 24 & 23.99999999995414936778 \\ 5 & 120 & 119.9999999999852711102 \\ 6 & 720 & 719.9999999999912556462 \\ 7 & 5040 & 5039.999999999991680317 \\ 8 & 40320 & 40319.99999999998847062 \\ 9 & 362880 & 362879.9999999999782436 \\ 10 & 3628800 & 3628799.999999999946786 \\ 11 & 39916800 & 39916799.99999999983748 \\ 12 & 479001600 & 479001599.9999999993984 \\ 13 & 6227020800 & 6227020799.999999997356 \\ 14 & 87178291200 & 87178291199.99999998649 \\ 15 & 1307674368000 & 1307674367999.999999921 \\ 16 & 20922789888000 & 20922789887999.99999948 \\ 17 & 355687428096000 & 355687428095999.9999961 \\ 18 & 6402373705728000 & 6402373705727999.999968 \\ 19 & 121645100408832000 & 121645100408831999.99971 \\ 20 & 2432902008176640000 & 2432902008176639999.9971 \end{array} \right)$$

Claude Leibovici
  • 289,558
  • From that link, it appears that formulation for n! is only accurate to about 6 digits when n is small and 12 digits when n is big. Good enough for integer solution, as you show; but I am not looking for an integer solution... within the range of doubles, an integer solution is easy enough with a 171 element lookup table. For non-integer (real) solution, my implementation of Lanczos approximation is giving me nearly full-precision double results. What I need is similarly near-double-precision accuracy of the inverse factorial and inverse lnfactorial functions. Can you improve it that much? – Brian Kennedy Feb 04 '19 at 08:34
  • @BrianKennedy. As I wrote, for the factorial, this could be improved using a few more terms and I am ready to do it (no guarantee for the fifteen digits but we could see). Concerning the inverse of the gamma function, waht I should do is one or two Newton iterations for the zero of $f(n)=\log (\Gamma (n))-k$ using for the derivative $f'(n)=\psi ^{(0)}(n)$ the same formula as used for the gamma function (whatever it is). I think it could be worth to try. – Claude Leibovici Feb 04 '19 at 09:45
  • Claude, if you can get 15 digits of precision out of a simpler formulation than the 13-series Lanczos approximation, I would certainly be interested... I think most of the world of people using the Lanczos approximation would. More importantly, there have been a lot of posts of people looking for inverse factorial / gamma and inverse log factorial / log gamma... and there is far less existing documentation on and experience with that. So, a good documented solution for that will be valuable to a lot more than me. – Brian Kennedy Feb 04 '19 at 16:10