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


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

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

10 Answers 10


JavaScript, 116 bytes

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


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

  //approximation for gamma

  //uncomment if more precision is needed

  //classic newton, gamma derivative is gamma*digamma


Test cases :

10 => 3.390062988090518
120 => 4.99999939151027
720 => 6.00000187248195
40320 => 8.000003557030217
3628800 => 10.000003941731514
  • \$\begingroup\$ Very nice answer altough it does not meet the required precision and it only works for numbers less than 706 \$\endgroup\$ Mar 5, 2014 at 22:06
  • \$\begingroup\$ @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 :) \$\endgroup\$
    – Michael M.
    Mar 5, 2014 at 22:20
  • \$\begingroup\$ Yes, now it's perfect, I voted it up :) \$\endgroup\$ Mar 5, 2014 at 22:22
  • 1
    \$\begingroup\$ You could save 1 char: n=prompt(M=Math) \$\endgroup\$
    – Florent
    Mar 6, 2014 at 13:33
  • \$\begingroup\$ Try running your code on a large number such as $10^{10^6}$, and ensuring that you get an integer result \$\endgroup\$ Apr 9, 2015 at 21:04

Mathematica - 74 54 49

The proper way will be

f[x_?NumberQ]:=NIntegrate[t^x E^-t,{t,0,∞}]

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.


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


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

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


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

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


Careful, operator precendence is weird sometimes.

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


And using the ! builtin you get 41.

An year overdue update:

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


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


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


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

  • \$\begingroup\$ I would be suprised if I would see someone beat this :o \$\endgroup\$ Mar 6, 2014 at 9:43
  • \$\begingroup\$ @JensRenders: I just did ;) \$\endgroup\$
    – mmumboss
    Mar 6, 2014 at 10:48
  • \$\begingroup\$ 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!). \$\endgroup\$
    – orion
    Mar 6, 2014 at 14:23
  • \$\begingroup\$ @mmumboss got you again :) \$\endgroup\$
    – orion
    Apr 13, 2015 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...


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;


  • \$\begingroup\$ 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" \$\endgroup\$ Mar 6, 2014 at 10:58
  • \$\begingroup\$ I don't see it in the ised and Mathematica solutions either? But I will look into it.. \$\endgroup\$
    – mmumboss
    Mar 6, 2014 at 11:04
  • 1
    \$\begingroup\$ I see the number 99 in the ised version, and the mathematica version is beaten anyways \$\endgroup\$ Mar 6, 2014 at 12:57
  • \$\begingroup\$ 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 ;) \$\endgroup\$
    – mmumboss
    Mar 6, 2014 at 13:03
  • \$\begingroup\$ But after changing inf to 99 does it meet the required precision? \$\endgroup\$
    – rubik
    Mar 6, 2014 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):
    while y<100:
    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):
    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):
    while y<100:
    return q if round(z,n)==x else f(x,n)
  • 2
    \$\begingroup\$ This is a code-golf question, so you need to provide the shortest answer, stating the length of your solution. \$\endgroup\$
    – VisioN
    Mar 5, 2014 at 18:02
  • \$\begingroup\$ 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. \$\endgroup\$ Mar 5, 2014 at 18:07
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – intx13
    Mar 5, 2014 at 18:13
  • \$\begingroup\$ 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 \$\endgroup\$ Mar 5, 2014 at 18:19
  • 2
    \$\begingroup\$ 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 ;) \$\endgroup\$
    – intx13
    Mar 5, 2014 at 20:15

Python 2.7 - 215 189 characters

f=lambda t:sum((x*.1)**t*2.71828**-(x*.1)*.1for x in range(999))
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)


# echo 6 | python invfact_golf.py
# echo 10 | python invfact_golf.py
# echo 3628800 | python invfact_golf.py

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.

  • \$\begingroup\$ I didn't know that site to count chars. I always use cat file | wc -c. \$\endgroup\$
    – rubik
    Mar 5, 2014 at 20:45
  • \$\begingroup\$ @rubik: oh nice, didn't think to use that. they both match =) \$\endgroup\$
    – Claudiu
    Mar 5, 2014 at 20:46

dg - 131 133 bytes

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 
$ dg gam.dg 
$ dg gam.dg 
$ dg gam.dg 
$ dg gam.dg  # i'm not really sure about this one :P it's instantaneous though
$ dg gam.dg  # a float example

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

  • \$\begingroup\$ Mine gives 42.8006566063, so they match within 5 digits of precision! \$\endgroup\$
    – Claudiu
    Mar 5, 2014 at 21:15
  • \$\begingroup\$ 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... \$\endgroup\$
    – rubik
    Mar 5, 2014 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.


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!


PARI/GP, 32 bytes

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


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;
        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)))))))))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.