Inverse gamma function? This is an analysis question I remember thinking about in high school. Reading some of the other topics here reminded me of this, and I'd like to hear other people's solutions to this.
We have the gamma function, which has a fairly elementary form as we all know,
$\Gamma(z) = \int_0^\infty e^{-t} t^{z-1} dt = \int_0^1 \left[ \ln(t^{-1}) \right]^{z-1}$
Which satisfies of course, $\Gamma(n) = (n-1)!$, $n\in \mathbb{N}$, and the various recurrence relations and other identities that we can all look up on wikipedia or mathworld or wherever.  We note that the gamma function is increasing on the interval $[a,\infty]$ where $a\approx 1.46163$.
The question is--can we come up with an explicit inverse function to the gamma function on this interval which looks similarly simple?
My techniques at the time were to write down a differential equation that the inverse would satisfy, and solve it, which I could do in terms of a power series expansion (being in high school, ignoring the issues of convergence) to get an approximate solution. But I was never able to get a very nice looking or exact solution.  I have a few more sophisticated tricks now to do this, but I would be interested to see how people with more experience with these kinds of questions would go about answering this.
The gamma function also satisfies a reasonable number of somewhat interesting looking functional relations like $\Gamma(z)\Gamma(1-z)=\pi/\sin(\pi z)$.  Does the inverse function satisfy any similar relations?
 A: This article addresses the question:
https://www.ams.org/journals/proc/2012-140-04/S0002-9939-2011-11023-2/S0002-9939-2011-11023-2.pdf
Uchiyama, Mitsuru, The principal inverse of the gamma function, Proc. Am. Math. Soc. 140, No. 4, 1343-1348 (2012). ZBL1245.33005, MR2869117.
A: A better numerical value of the positive zero of the digamma function can be obtained by using the built-in Mathematica function FindInstance:
$k \approx 1.46163214496836234126265954233$
Then a better numerical value of $c$ is found and used:
$c=\sqrt{2\pi}/e - \Gamma(k)$
$c \approx 0.0365338144849004166003358471688$
Then the relation mentioned in other answers can be used:
$\Gamma_{approx}^{-1}(n)=L(x)/W(\frac{L(x)}{e}) + \frac{1}{2}$
However a more accurate real numerical value of the inverse gamma function is found by directly using FindInstance:
invgamnum[n_] = 
  FindInstance[Gamma[x] == n, x, Reals, WorkingPrecision -> 30]

The numerical value obtained was set to 30 decimal digits; this can be changed to an arbitrary numerical precision.
If Reals is not added in the code above, FindInstance usually yields a complex valued solution.
FindInstance can be used to find more than one solution (if they exist).
For example, for $\Gamma(x) = e$ , with Mathematica:
FindInstance[Gamma[x] == E, x, Reals, 2, WorkingPrecision -> 30]

The two numerical values obtained are:
$x=\Gamma^{-1}(e) \approx 0.328713113325854368520566502893$
$x=\Gamma^{-1}(e) \approx 3.31244088253916037027609875238$
Note that if the relation involving the Lambert W function were used, a less accurate single numerical value would be found:
$\Gamma^{-1}(e) \approx 3.3111260003863245271459515483$
A: Mathematica has an inverse gamma function. It is on the web page on special functions. This would suggest that the problem is at least simple enough for computer implementation.
I have just found more material on the inverse of the regularized incomplete gamma function from Mathematica. There are downloads on the site with information as well. Some of them are Mathematica notebooks and need the player which is free to be opened.
The information includes differential equations, representations through equivalent functions and series representations among other things.
A: David Cantrell gives a good approximation of $\Gamma^{-1}(n)$ on this page.
I'll copy the result here in case that page ever goes down:
$k$ = the positive zero of the digamma function, approximately $1.461632$
$c$ = $\sqrt{2\pi}/e - \Gamma(k)$, approximately $0.036534$
$L(x)$ = $\ln(\frac{x+c}{\sqrt{2\pi}})$
$W(x)$ = Lambert W function
$ApproxInvGamma(x)$ = $L(x)/W(\frac{L(x)}{e}) + \frac{1}{2}$
A: This should work for Mathematica:
c = 0.036534
l[x_] = Log[(x + c)/Sqrt[2*Pi]]
aig[x_] = l[x]/(ProductLog[l[x]/E]) + 1/2

A: 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!$


*

*If $x=1$ return $1$

*Count the number of trailing zeros, $t$, in the binary representation of $x$

*Calculate $w=\frac{x}{t!}$

*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 \frac{\ln(w)}{\ln(t+1)} \right \rfloor$$

*Return $t+r$

A: Here is some code I wrote for a similar problem a few years back. I used Stirling's approximation and Newton's method. It may not be the most efficient solution but it worked fine for my applications. Maybe this could help someone in the future. I would also like to emphasize that this is only an approximation and that it may be possible to get the incorrect answer if your input is sufficiently large (However I have not observed this)
Stirlings approximation: https://en.wikipedia.org/wiki/Stirling%27s_approximation
Newton's Method: https://en.wikipedia.org/wiki/Newton%27s_method
import math

q = math.log10(int(input())
def function(x):
    y = .5*math.log10(2*math.pi*x) + x*math.log10(x)-x*math.log10(math.e)-q
    return y

def derivative(x):
    y1 = ( (math.log( (1/math.e) * x) / math.log(10)) + (1/(2*math.log(10)*x)) + (1/math.log(10)))
    return y1

x0 = len(str(factorial_n))
arrX = [x0]
for i in range(0,3):
    tempx = arrX[-1]
    X_n = tempx - (function(tempx) / derivative(tempx))
    arrX.append(X_n)
print(math.floor(arrX[-1]))

