Calculate the inverse of factorial Write the shortest code that will take any real number greater than 1 as input and will output its positive inverse factorial. In other words, it answers the question "what number factorial is equal to this number?". Use the Gamma function to extend the definition for factorial to any real number as described here.

For example:  

input=6 output=3 
input=10 output=3.390077654


because 3! = 6 and 3.390077654! = 10  

Rules  


It is forbidden to use built in factorial functions or gamma functions, or functions that rely on these functions.
The program should be able to calculate it to 5 decimal digits, with the theoretical ability to calculate it to any precision (It should contain a number that can be made arbitrary big or small to get arbitrary precision)
Any language is allowed, shortest code in characters wins.


I made a working example here. Have a look.
Javascript (116)

Black magics here ! Gives a result in few milliseconds.
Only elementary math functions used : ln, pow, exponential

x=9;n=prompt(M=Math);for(i=1e4;i--;)x+=(n/M.exp(-x)/M.pow(x,x-.5)/2.5066/(1+1/12/x+1/288/x/x)-1)/M.log(x);alert(x-1)


Too bad LaTeX is not supported on codegolf but basically, I coded a newton solver for f(y)=gamma(y)-n=0 and x=y-1 (because x! is gamma(x+1)) and approximations for gamma and digamma functions.

Gamma approximation is Stirling approximation
Digamma approximation use Euler Maclaurin formula
The digamma function is the derivative of gamma function divided by gamma function : f'(y)=gamma(y)*digamma(y)  

Ungolfed :

n = parseInt(prompt());
x = 9; //first guess, whatever but not too high (&lt;500 seems good)

//10000 iterations
for(i=0;i&lt;10000;i++) {

  //approximation for digamma
  d=Math.log(x);

  //approximation for gamma
  g=Math.exp(-x)*Math.pow(x,x-0.5)*Math.sqrt(Math.PI*2)*(1+1/12/x+1/288/x/x);

  //uncomment if more precision is needed
  //d=Math.log(x)-1/2/x-1/12/x/x+120/x/x/x/x;
  //g=Math.exp(-x)*Math.pow(x,x-0.5)*Math.sqrt(Math.PI*2)*(1+1/12/x+1/288/x/x-139/51840/x/x/x);

  //classic newton, gamma derivative is gamma*digamma
  x-=(g-n)/(g*d);
}

alert(x-1);


Test cases :

10 =&gt; 3.390062988090518
120 =&gt; 4.99999939151027
720 =&gt; 6.00000187248195
40320 =&gt; 8.000003557030217
3628800 =&gt; 10.000003941731514