从a mathforum post开始：
设k表示二γ函数的正零点，约为1.461632...
c =平方（2pi）/e -伽玛（k）...
设L（x）= ln（（x+c）/Sqrt（2pi）），则伽马近似的逆为
近似值为反伽玛或AIG（x）= L（x）/ W（L（x）/ e）+ 1/2。

k <- 1.461632
cc <- sqrt(2*pi)/exp(1)-gamma(k)
L <- function(x) {
    log((x+cc)/sqrt(2*pi))
}
AIG <- function(x) {
    Lx <- L(x)
    Lx/(emdbook::lambertW(Lx*exp(-1))) + 1/2
}

par(las=1,bty="l")
curve(1-AIG(gamma(x))/x,from=2,to=20,
      ylab="relative error of approximation")

或者，您可以使用uniroot()：

AIG(5)
ufun <- function(x=5) {
   uniroot(function(z) gamma(z)-x,c(1.00001,10))$root
}
ufun(5) ## 3.852341
AIG(5)  ## 3.848149