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

• This could use some more test cases, in particular to cover zero and negative inputs. – Peter Taylor Mar 5 '14 at 16:43
• I edited that the input should be greater than 1 because otherwise there could be muliple answers. – Jens Renders Mar 5 '14 at 16:45
• There can be multiple answers anyway unless you also add a requirement that the output must be greater than 1. – Peter Taylor Mar 5 '14 at 17:11
• Your working example gives 3.99999 when inputted 24. So is such a solution acceptable? – rubik Mar 5 '14 at 20:07
• yes because this can be seen as 4, to 5 decimal places correct – Jens Renders Mar 5 '14 at 20:11

## 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 (<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

• Very nice answer altough it does not meet the required precision and it only works for numbers less than 706 – Jens Renders Mar 5 '14 at 22:06
• @JensRenders, well, I've added some iterations of the newton solver, changed the initial guess and a better approximation for gamma function. That should fits the rules now. Let me now if it's ok :) – Michael M. Mar 5 '14 at 22:20
• Yes, now it's perfect, I voted it up :) – Jens Renders Mar 5 '14 at 22:22
• You could save 1 char: n=prompt(M=Math) – Florent Mar 6 '14 at 13:33
• Try running your code on a large number such as $10^{10^6}$, and ensuring that you get an integer result – David G. Stork Apr 9 '15 at 21:04

## Mathematica - 7454 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]

• Ahh all those built-in functions! I don't think this is beatable except in a similar way. – rubik Mar 5 '14 at 21:15
• Mathematica is so unfair for math stuff ! :D – Michael M. Mar 5 '14 at 21:24
• from the name itself MATHematica – Dadan Mar 6 '14 at 10:38
• Is NumberQ pattern test required? Or parens in E^(-t)? Is it cheating to turn NIntegrate to Integrate? Probably... :) – orion Mar 7 '14 at 9:31
• It's turning into a real challenge ;) – mmumboss Mar 7 '14 at 13:03

# 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. • I would be suprised if I would see someone beat this :o – Jens Renders Mar 6 '14 at 9:43 • @JensRenders: I just did ;) – mmumboss Mar 6 '14 at 10:48 • To clarify the discussion about accuracy: it's set by .1 (integration step) and 99 (integration limit). Bisection goes to machine precision. The bisection limit @1,99 can be kept at 99 unless you want to input numbers above (99!). – orion Mar 6 '14 at 14:23 • @mmumboss got you again :) – orion Apr 13 '15 at 11:43 ## 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. • you should specify the option for accuracy as this is required in the rules! "It should contain a number that can be made arbitrary big or small to get arbitrary precision" – Jens Renders Mar 6 '14 at 10:58 • I don't see it in the ised and Mathematica solutions either? But I will look into it.. – mmumboss Mar 6 '14 at 11:04 • I see the number 99 in the ised version, and the mathematica version is beaten anyways – Jens Renders Mar 6 '14 at 12:57 • Given the position in the code, this is probably the upper limit for the integral. In my code this is inf. So yeah, if I change this inf into 99, my answer becomes less accurate, which means that this number influences the precision, and therefore I meet the rules. If I change it to 99 I even save a char ;) – mmumboss Mar 6 '14 at 13:03 • But after changing inf to 99 does it meet the required precision? – rubik Mar 6 '14 at 13:11 # 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)  • This is a code-golf question, so you need to provide the shortest answer, stating the length of your solution. – VisioN Mar 5 '14 at 18:02 • A nice method but the problem is that you can not guarrantee that this will ever find the answer... Also, this is codegolf zo you might try to minimalise character use. – Jens Renders Mar 5 '14 at 18:07 • Python's random() uses a Mersenne Twister which I believe walks the space of Python's floats, so it should always terminate if there is an answer within the tolerance. – intx13 Mar 5 '14 at 18:13 • Do you mean that it returns every float value before repeating one? if that's the case than this code would be valid if you where able to overcome the stack overflow – Jens Renders Mar 5 '14 at 18:19 • The code is capable, it's just that you and I may not have the time nor the computer resources to execute it to completion ;) – intx13 Mar 5 '14 at 20:15 # Python 2.7 - 215189 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. • I didn't know that site to count chars. I always use cat file | wc -c. – rubik Mar 5 '14 at 20:45 • @rubik: oh nice, didn't think to use that. they both match =) – Claudiu Mar 5 '14 at 20:46 # 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! • Mine gives 42.8006566063, so they match within 5 digits of precision! – Claudiu Mar 5 '14 at 21:15 • That's great! I don't know where the upper bound is, but it should break somewhere. For 1e100 it gives: 69.95780520000001, for 1e150 it outputs 96.10586423000002, whereas for 1e200 it blows up. But really I don't know whether those results are reliable... – rubik Mar 5 '14 at 21:18 # 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!

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