What approximations for the Gamma function's inverse appear to work 'best'?

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

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

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

$$x¡=m-k+\frac{\ln(\frac{x}{(m-k)!})}{\ln(m-k+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

$$r+\frac{\ln(\frac{x}{r!})}{\ln(r+1)}$$

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

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

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

$$x¡=r+\frac{\ln(\frac{x}{r!})}{\ln(r+\frac{1}{2})}$$

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.

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:

$$\frac{\Gamma'(z)}{\Gamma(z)}=\psi(z)$$

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:

$$z(1)=2$$

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};
                    return(log(z+k+1)-1/2/(z+k+1)-1/12/(z+k+1)^2-sum(1/(z+0:k)))};
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];
};
plot(g,z,type="l",col="red",xlab="Gamma(z)",ylab="z")
ge<-c(1,2,6,24,120,720,5040);
ze<-c(2,3,4,5,6,7,8);
points(ge,ze,pch=16)
z[which(g == 6)]
z[which(g == 24)]
z[which(g == 120)]
z[which(g == 720)]
z[which(g == 5040)]

The results are:

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