So I was wondering how we approximate the inverse of the Gamma function, where I tried a few methods:

Lagrange inversion theorem:

$$\Gamma^{-1}(z)=a+\sum_{n=1}^{\infty}\lim_{w\to a}\frac{(z-\Gamma(a))^n}{\Gamma(n+1)}\frac{d^{n-1}}{dw^{n-1}}\left(\frac{w-a}{\Gamma(w)-\Gamma(a)}\right)^n$$

I also found that one could find the inverse of Stirling's approximation using the Lambert W function.

And of course, there are other ones, but which method gets the most accuracy for possibly easier calculations?

  • 1
    $\begingroup$ How exactly does one solve for $\bigg(\dfrac ne\bigg)^n\sqrt n=a$ in terms of Lambert's W function ? $\endgroup$
    – Lucian
    Jan 24, 2016 at 1:04
  • 2
    $\begingroup$ @Lucian math.stackexchange.com/questions/430167/… $\endgroup$ Jan 24, 2016 at 1:18
  • $\begingroup$ Is this question about the principal inverse of the gamma function as defined in this paper by Mitsuru Uchiyama? $\endgroup$
    – njuffa
    Jan 24, 2016 at 1:23
  • $\begingroup$ @njuffa Yes. I would think it rather difficult to find for all branches. I'll be looking through that paper. $\endgroup$ Jan 24, 2016 at 1:36

2 Answers 2


Factorial inverse approximation is given by this formula split into couple of steps for simplicity

$$m=\frac{\ln(x)}{\ln(\ln(x)+1)} $$



Integer version

$$m=\left \lceil \frac{\ln(x)}{\ln(\ln(x)+1)} \right \rceil$$

$$k= \frac{\ln(\frac{m!}{x})}{\ln(m+1)}$$

$$x¡=\left \lfloor m-k+ \frac{\ln(\frac{x}{(m-k)!})}{\ln(m-k+1)} \right \rceil$$

(Mind the ceiling/floor marking, the last is round off, first ceiling is for better convergence only and not to have to calculate non-integer factorial twice.)

The idea is to reach $y$ from $x=y!$ as close as possible using $\frac{\ln(x)}{\ln(\ln(x)+1)}$ then, since this is undershooting, to estimate how low we are and how many times we need to multiply by $(m+1)(m+2)...(m+k)$ more, rounding this by $(m+1)^k$ since the numbers are close, finally assuming we are close enough using the estimation $(n+\epsilon)! \sim n!(n+1)^{\epsilon}$ from Gamma limit we get to the final result.

(We marked inverse factorial with inverse exclamation mark.)

The formula works for $x>1$, however, simple $\frac{(x+1)!}{x+1}=x!$ can shift it upwards.

The main advantage is that the method is giving quite a precise evaluation even for non-integer values already, certainly sufficiently for integer factorial even beyond $10^7!$. However if a better precision is needed it is sufficient to repeat the last step


It could be some special task that would require repeating this more than once, likely for small values but the step can be repeated as many times as needed.

High precision version

$$m=\frac{\ln(x)}{\ln(\ln(x)+1)} $$




Integer binary algorithm

Assume that we have integer input $x=n!$

  1. If $x=1$ return $1$
  2. Count the number of trailing zeros, $t$, in the binary representation of $x$
  3. Calculate $w=\frac{x}{t!}$
  4. Calculate the number of times $t+1$ divides $w$, rounded down to the nearest integer giving this result, but simplified and optimized for whatever programming language is used $$r= \left \lfloor \log_{t+1}(w) \right \rfloor$$
  5. Return $t+r$

4.a. If 4. is to be implemented in form of a loop, instead of the above, we can start dividing $w$ with $t+1$, $w_1=\frac{w}{t+1}$. We keep on dividing $w_k=\frac{w_{k-1}}{t+k}$ until reach 1. The result is the highest $t+k$ we have reached.

  • $\begingroup$ In the first approximation, you have ln(m - k + 1) in the denominator. Then you say you can repeat the last step and you show the expression with 'r' substituting 'm - k', such that you have ln(r + 1) in the denominator. Makes sense. But then in the higher-precision version where you do that step twice, you have ln(m - k + 1/2) and ln(r + 1/2) in the denominators. Which is correct: +1 or +1/2? Or am I missing something? $\endgroup$ Feb 4, 2019 at 2:51
  • 1
    $\begingroup$ @BrianKennedy Those are two different methods. There is no ceiling floor in the version with 1/2. Integer method means: using integer evaluations. $\endgroup$
    – Alex Peter
    Feb 8, 2019 at 16:24

I was wondering about that for a year or so, but since I found this question, I decided to share my thoughts here.

Not sure if this is the best way or not, but my idea was to use the digamma function. Why? Because it's rather easy to compute. I use asymptotic expansion (the first $3$ terms only) and shift the variable to higher values:

$$\psi(z) \approx \log (z) -\frac{1}{2z}-\frac{1}{12 z^2}, \qquad z \geq z_m$$

$$\psi(z) = \psi(z+k)-\sum_{l=0}^{k-1} \frac{1}{z+l}, \qquad z < z_m, \quad k= \lceil z_m-z \rceil$$

As for the inverse gamma, I simply use the definition of digamma:


Inverting we have:

$$\frac{dz}{d \Gamma}=\frac{1}{\psi(z) \Gamma}$$

Now use any suitable numerical scheme with initial value we want, for example:


Here's an (improved) example in R, using the two step explicit scheme which appears to give the best agreement.

z0 <- 2;
zm <- 15;
g0 <- 1;
gm <- 5041;
psi <- function(z){if(z < zm){k <- floor(zm-z)}else{k <- 0};
h <- 1/20;
N <- floor((gm-g0)/h);
z <- 1:(N+1);
g <- g0+(0:N)*h;
z[1] <- z0;
z[2] <- z[1] + h/psi(z[1])/g[1];
z[2] <- z[1] + h/psi(z[2])/g[2];
n <- 2;
while (n <= N){n <- n+1;
                z[n] <- z[n-2] + 2*h/psi(z[n-1])/g[n-1];
       z[n] <- z[n-1] + h/psi(z[n])/g[n-1];
z[which(g == 6)]
z[which(g == 24)]
z[which(g == 120)]
z[which(g == 720)]
z[which(g == 5040)]

The results are:

enter image description here

> z[which(g == 6)]
[1] 3.999231
> z[which(g == 24)]
[1] 5.001148
> z[which(g == 120)]
[1] 6.001462
> z[which(g == 720)]
[1] 7.00142
> z[which(g == 5040)]
[1] 8.001333

So the error is about $10^{-3}$.

Remark: Of course, I intentionally started with the most simple interval where gamma is monotone, so the inverse will be single-valued.


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.