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
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)$$