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?

  • $\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 '16 at 1:04
  • 1
    $\begingroup$ @Lucian math.stackexchange.com/questions/430167/… $\endgroup$ – Simply Beautiful Art Jan 24 '16 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 '16 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$ – Simply Beautiful Art Jan 24 '16 at 1:36

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 leapfrog method.

z0 <- 2;
zm <- 20;
g0 <- 1;
gm <- 24;
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];
n <- 2;
while (n <= N){n <- n+1;
                z[n] <- z[n-2] + 2*h/psi(z[n-1])/g[n-1];

The result is:

enter image description here

As can be seen, the leapfrog scheme causes oscillations, however the mean value stays true to the correct function, which is better than the direct Euler integration.

We can prove this by using the arithmetic means of even and odd elements of vector z:

z <- (shift(z, n=1L, fill=z0)+z)/2;

The 'shift' function belongs to the package data.table, in case anyone would want to check the code.

The result is (for the same value of step h as before):

enter image description here

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


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$ – Brian Kennedy Feb 4 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 at 16:24

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.