Calculate the inverse of factorial

Question

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.

Answer 1

JavaScript, 116 bytes

Black magic here! Gives a result in a few milliseconds.
Only elementary math functions are 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)

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 uses Euler Maclaurin formula. The digamma function is the derivative of gamma function divided by gamma function: \$ f'(y)=\Gamma(y) \cdot \psi(y) \$

Ungolfed :

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

//10000 iterations
for(i=0;i<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 => 3.390062988090518
120 => 4.99999939151027
720 => 6.00000187248195
40320 => 8.000003557030217
3628800 => 10.000003941731514

Answer 2

Mathematica - 74 54 49

The proper way will be

f[x_?NumberQ]:=NIntegrate[t^x E^-t,{t,0,∞}]
x/.FindRoot[f@x-Input[],{x,1}]

If we just drop the test ?NumberQ it would still work, but throw some nasty warnings, which would go away if we switch to symbolic integration Integrate, but this would be illegal (I suppose), because the function would automatically be converted to Gamma function. Also we can get rid of external function that way.

Anyway

x/.FindRoot[Integrate[t^x E^-t,{t,0,∞}]-Input[],{x,1}]

To heck with proper input, just function definition (can't let MatLab win)

x/.FindRoot[Integrate[t^x E^-t,{t,0,∞}]-#,{x,1}]&

If built-in factorial were allowed

N@InverseFunction[#!&]@Input[]

The above does not give an integer (which is the argument for a true factorial function). The following does:

Floor[InverseFunction[Gamma][n]-1]

Answer 3

ised: 72 46 characters

This is almost a perfect fit... there is a "language" out there that seems to be meant precisely for math golf: ised. Its obfuscated syntax makes for a very short code (no named variables, just integer memory slots and a lot of versatile single char operators). Defining the gamma function using an integral, I got it to 80 seemingly random characters

@4{:.1*@+{@3[.,.1,99]^x:*exp-$3}:}@6{:@{$4::@5avg${0,1}>$2}$5:}@0,0@1,99;$6:::.

Here, memory slot $4 is a factorial function, memory slot $6 bisection function and memory slot $2 is expected to be set to input (given before sourcing this code). Slots $0 and $1 are the bisection boundaries. Call example (assuming above code is in file inversefactorial.ised)

bash> ised '@2{556}' --f inversefactorial.ised
556
5.86118

Of course, you could use the builtin ! operator, in which case you get down to 45 characters

@6{:@{{@5avg${0,1}}!>$2}$5:}@0,0@1,99;$6:::.

Careful, operator precendence is weird sometimes.

Edit: remembered to inline the functions instead of saving them. Beat Mathematica with 72 characters!

@0,0@1,99;{:@{{:.1*@+{@3[.,.1,99]^x:*exp-$3}:}::@5avg${0,1}>$2}$5:}:::.

And using the ! builtin you get 41.

---

An year overdue update:

I just realized this was highly inefficient. Golfed-down to 60 characters:

@0#@1,99;{:@{.1*@3[.,.1,99]^@5avg${0,1}@:exp-$3>$2}$5:}:::.

If utf-8 is used (Mathematica does it too), we get to 57:

@0#@1,99;{:@{.1*@3[.,.1,99]^@5avg${0,1}·exp-$3>$2}$5:}∙.

A bit different rewrite can cut it down to 46 (or 27 if using builtin !):

{:x_S{.5@3[.,.1,99]^avgx·exp-$3*.1<$2}:}∙∓99_0

The last two characters can be removed if you are ok with having the answer printed twice.

Answer 4

MATLAB 54 47

If I pick the right challenges MATLAB is really nice for golfing :). In my code I find the solution to the equation (u-x!)=0 in which u is the user input, and x the variable to solve. This means that u=6 will lead to x=3, etc...

@(x)fsolve(@(y)u-quad(@(x)x.^y./exp(x),0,99),1)

The accuracy can be changed by altering the upper limit of the integral, which is set at 99. Lowering this will change the accuracy of the output as follows. For example for an input of 10:

upper limit = 99; answer = 3.390077650833145;
upper limit = 20; answer = 3.390082293675363;
upper limit = 10; answer = 3.402035336604546;
upper limit = 05; answer = 3.747303578099607;

etc.

Answer 5

Python - 199 chars

Ok, so you'll need a lot of stack space and a lot of time, but hey, it'll get there!

from random import *
from math import e
def f(x,n):
    q=randint(0,x)+random()
    z=0
    d=0.1**n
    y=d
    while y<100:
            z+=y**q*e**(-y)*d
            y+=d
    return q if round(z,n)==x else f(x,n)

Here's another approach with even more recursion.

from random import *
from math import e
def f(x,n):
    q=randint(0,x)+random()
    return q if round(h(q,0,0.1**n,0),n)==x else f(x,n)
def h(q,z,d,y):
    if y>100:return z
    else:return h(q,z+y**q*e**(-y)*d,d,y+d)

Both of these can be tested with >>>f(10,1) provided you set the recursion limit around 10000. More than one decimal place of accuracy will likely not complete with any realistic recursion limit.

Incorporating the comments and a few modifications, down to 199 chars.

from random import*
from math import*
def f(x,n):
    q=random()*x+random()
    z=y=0
    while y<100:
            z+=y**q*e**-y*0.1**n
            y+=0.1**n
    return q if round(z,n)==x else f(x,n)

Answer 6

Python 2.7 - 215 189 characters

f=lambda t:sum((x*.1)**t*2.71828**-(x*.1)*.1for x in range(999))
n=float(raw_input());x=1.;F=0;C=99
while 1:
 if abs(n-f(x))<1e-5:print x;break
 F,C,x=f(x)<n and(x,C,(x+C)/2)or(F,x,(x+F)/2)

Usage:

# echo 6 | python invfact_golf.py
2.99999904633
# echo 10 | python invfact_golf.py
3.39007514715
# echo 3628800 | python invfact_golf.py
9.99999685376

To change precision: change 1e-5 to a smaller number for greater precision, larger number for worse precision. For better precision you probably want to give a better value for e.

This just implements the factorial function as f, and then does a binary search to hone in on the most accurate value of the inverse of the input. Assumes the answer is less than or equal to 99 (it wouldn't work for an answer of 365 for sure, I get a math overflow error). Very reasonable space and time usage, always terminates.

Alternatively, replace if abs(n-f(x))<=10**-5: print x;break with print x to shave off 50 characters. It'll loop forever, giving you a more and more accurate estimate. Not sure if this would fit with the rules though.

Answer 7

dg - 131 133 bytes

o,d,n=0,0.1,float$input!
for w in(-2..9)=>while(sum$map(i->d*(i*d)**(o+ 10**(-w))/(2.718281**(i*d)))(0..999))<n=>o+=10**(-w)
print o

Since dg produces CPython bytecode this should count for Python as well, but oh... Some examples:

$ dg gam.dg 
10
3.3900766499999984
$ dg gam.dg 
24
3.9999989799999995
$ dg gam.dg 
100
4.892517629999997
$ dg gam.dg 
12637326743
13.27087070999999
$ dg gam.dg  # i'm not really sure about this one :P it's instantaneous though
28492739842739428347929842398472934929234239432948923
42.800660880000066
$ dg gam.dg  # a float example
284253.232359
8.891269689999989

EDIT: Added two bytes because I didn't remember that it should accept floats as well!

Answer 8

R, 92 bytes

A function, g, which takes input z and outputs the inverse factorial of that number

_{There is almost certainly more to be golfed out of this, so if you see something that I can improve, please let me know.}

library(pryr)
g=f(z,uniroot(f(a,integrate(f(x,x^a*exp(-x)),0,Inf)$v-z),c(0,z+1),tol=1e-9)$r)

Try it online!

Ungolfed and Commented

library(pryr)                     # Add pryr to workspace
inv.factorial = f(z,              # Declare function which  
  uniroot(                        # Finds the root of
    f(a, integrate(               # The integral of 
      f(x, x^a*exp(-x))           # The gamma function
        ,0 ,Inf                   # From 0 to Infinity
      )$value-z                   # Minus the input value, `z`
    ), c(0, z+1),                 # On the bound of 0 to z+1
    tol = 1e-323                  # With a specified tolerance
  )$root                          # And outputs the root
)                                 # End function

Try it online!

Answer 9

PARI/GP, 32 bytes

If we can use solve() and \$\Gamma\$ function gamma():

f(y)=solve(x=0,100,gamma(x+1)-y)

Try it online!

Answer 10

Javascript (without using loops!)

To do this, I used a well-known numerical approximation of the inverse of the Stirling Factorial approximation, (and also got inspiration by this ..cough.. cough.. code of someone else...)

function f(n){
    if(n==1) return 1;
    else if(n==2) return 2;
    else if(n==6) return 3;
    else if(n==24) return 4;
    else{
        return Math.round((((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))/Math.log((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))))/Math.log((((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))/Math.log((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))))))/Math.log((((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))/Math.log((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))))/Math.log((((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592))+(Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))/Math.log((Math.log(n)/Math.LN10 *  Math.log(10.) - .5 * Math.log(2.*3.141592)))))))))
    }
}