1
$\begingroup$

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:

http://mathforum.org/kb/message.jspa?messageID=342551&tstart=0 (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
$\endgroup$
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Brian Kennedy Feb 4 at 8:34
  • $\begingroup$ @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. $\endgroup$ – Claude Leibovici Feb 4 at 9:45
  • $\begingroup$ 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. $\endgroup$ – Brian Kennedy Feb 4 at 16:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.