The Pi function is an extension of the factorial over the reals (or even complex numbers). For integers n, Π(n) = n!, but to get a definition over the reals we define it using an integral:

Pi(z) = integral t from 0 to infinity e^-t t^z dt

In this challenge we will invert the Π function.

Given a real number z ≥ 1, find positive x such that Π(x) = z. Your answer must be accurate for at least 5 decimal digits.


120 -> 5.0000
10 -> 3.39008
3.14 -> 2.44815
2017 -> 6.53847
1.5 -> 1.66277
  • 4
    \$\begingroup\$ Note that more often people use the Gamma (Γ) function. Π(x) = Γ(x+1). But IMO Γ is a shifted abomination, and Π is the true extension of the factorial. \$\endgroup\$
    – orlp
    Commented Jan 27, 2017 at 19:12
  • 1
    \$\begingroup\$ Wellp, that series expansion is enough to scare me... i.imgur.com/ttgzDSJ.gif \$\endgroup\$ Commented Jan 27, 2017 at 19:22
  • 1
    \$\begingroup\$ All of the examples you give have other solutions as well, for example 120 -> -0.991706. This is because Π(x) goes to infinity as x goes to -1 from the right. Perhaps you mean to insist that x>0 as well. \$\endgroup\$ Commented Jan 27, 2017 at 19:28
  • \$\begingroup\$ @GregMartin Added as well. \$\endgroup\$
    – orlp
    Commented Jan 27, 2017 at 19:29
  • 1
    \$\begingroup\$ There are some reasons to prefer the shifted version, despite it seems unnatural. See e.g. this answer on MathOverflow as well as others on that page. \$\endgroup\$
    – Ruslan
    Commented Jan 28, 2017 at 17:23

11 Answers 11


Mathematica, 17 15 27 bytes


Output looks like {{x -> n}}, where n is the solution, which may not be allowed.


Pyth, 4 bytes


A program that takes input of a number and prints the result.

Test suite

How it works

.I.!    Program. Input: Q
.I.!GQ  Implicit variable fill
.I      Find x such that:
  .!G    gamma(x+1)
     Q   == Q
        Implicitly print

MATL, 13 bytes


This uses linear seach in steps of 1e-5 starting at 1. So it is terribly slow, and times out in the online compiler.

To test it, the following link replaces the 1e-5 accuracy requirement by 1e-2. Try it online!


1        % Push 1 (initial value)
`        % Do...while
  1e-5   %   Push 1e-5
  +      %   Add
  t      %   Duplicate
  QYg    %   Pi function (increase by 1, apply gamma function)
  G<     %   Is it less than the input? If so: next iteration
         % End (implicit)
         % Display (implicit)

GeoGebra, 25 bytes


Entered in the CAS input, and expects input of a number in spreadsheet cell A1. Returns a one-element array of the form {x = <result>}.

Here is a gif of the execution:

Execution of progrma

How it works

NumericallySolve the following equation: Gamma(x+1)=A1, with starting value x=1.

  • \$\begingroup\$ Is it guaranteed to return a positive number, and does it work for 1.5, which has broken several answers? \$\endgroup\$
    – Pavel
    Commented Jan 28, 2017 at 3:39
  • \$\begingroup\$ @Pavel I can confirm that it works for 1.5. I haven't been able to find out which algorithm GeoGebra uses for numerical solving, but the initial value of x=1 has given purely positive answers for every value I have tried. \$\endgroup\$ Commented Jan 29, 2017 at 17:46

MATLAB, 59 bytes


This is an anonymous function that finds the minimizer of the squared difference betweeen the Pi function and its input, starting at 1, with very small tolerance (given by eps) to achieve the desired precision.

Test cases (run on Matlab R2015b):

>> @(x)fminsearch(@(t)(gamma(t+1)-x)^2,1,optimset('TolF',eps))
ans = 
>> f = ans; format long; f(120), f(10), f(3.14), f(2017)
ans =
ans =
ans =
ans =

You could try it online in Octave, but unfortunately some of the results lack the required precision.


J, 86 33 bytes


Uses Newton's method with log Pi to avoid overflow.

This is the previous version that computes log Gamma using Stirling's approximation. The step size (1e3) and number of terms in log Gamma (3) can be increased for possibly higher accuracy at the cost of performance.

3 :'(-(k-~g)%%&1e3(g=:((%~12 _360 1260 p.&:%*:)+-+^~-&^.%:@%&2p1)@>:)D:1])^:_>:k=:^.y'

Another version that computes the coefficient terms on the fly

3 :'(-((-^.y)+g)%%&1e3(g=:((%~(((%1-^@-)t:%]*<:)+:>:i.3)p.%@*:)+(*^.)-]+-:@^.@%&2p1)@>:)D:1])^:_>:^.y'

Try it online! and to see the terms converge.


((]-(-~^.@!)%[:^.@!D.1])^:_>:)@^.  Input: float y
(                            )@^.  Operate on log(y)
                           >:        Increment, the initial guess is log(y)+1
 (                     )^:_          Repeat until convergence starting with x = log(y)+1
                      ]                Get x
               ^.@!                    The log Pi verb
             [:    D.1                 Approximate its first derivative at x
       ^.@!                            Apply log Pi to x
     -~                                Subtract log(y) from it
            %                          Divide it by the derivative
  ]-                                   Subtract it from x and use as next value of x

Mathematica, 21 bytes


FindRoot applies Newton's method internally when there is an initial value.

The two methods below apply Newton's method directly.

Alternative using FixedPoint 45 bytes


A more precise implementation of Newton's method for solving this since Mathematica can compute the derivative directly instead of approximating it.

Using rules to replace repeatedly would be shorter, but there is a limit (65536) to how many iterations it can perform that might be hit whereas FixedPoint does not have a limit.

Alternative using rules, 38 bytes




Jelly, 34 bytes


Try it online! or View the intermediate values as they converge.

An implementation of J's combination of Newton's method and derivative approximation (secant method) to compute the inverse of Π(n).

It solves for the inverse of log(Π(n)) instead in order to avoid overflow.

It starts with an initial guess x0 = y+1 where y = log(Π(n)). Then it iterates until convergence using xn+1 = xn - (log(Π(xn)) - y) / (log((Π(1.001 * xn)) - log(Π(xn))) / (0.001 * xn)).

  • 3
    \$\begingroup\$ I get an error with input 1.5 \$\endgroup\$
    – Luis Mendo
    Commented Jan 28, 2017 at 0:35
  • \$\begingroup\$ @LuisMendo Wow that is a good catch! It occurs since one of the intermediate values is ~65807 which is a huge value after gamma is applied, and Python overflows. The same occurs in J since it relies on the same computation. \$\endgroup\$
    – miles
    Commented Jan 28, 2017 at 0:53

PARI/GP, 30 bytes


Finds solution between 1 andx+1. Unfortunately, x is not big enough as upper bound for input like 1.5.


Mathematica, 26 Bytes

Yet another Mathematica solution!

Equation solving can always be turned into a minimization problem.


Finds the argument that minimizes the difference between the left and right sides of the equation.

Using NArgMin rather than NMinimize forces the output to just be the desired result rather than the usual verbose rule-based output (and it saves a byte!)


C with libm, 111

Update - fixed for input 1.5.

f(double *z){double u=2**z,l=0,g=u,p=0;for(;log(fabs(g-p))>-14;)p=g,g=(u+l)/2,u=tgamma(g+1)>*z?g:(l=g,u);*z=g;}

gamma(x+1) is a monotonically increasing function over the range in question, shis is just a binary search until the difference between successive values is small. The starting lower bound is 0 and the starting upper bound is 2*x.

Input and output is via a pointer to a double passed to the function.

I'm pretty sure this can be golfed deeper - in particular I don't think I need 4 local doubles, but so far I'm not seeing an easy way to reduce this.

Try it online - Builds (linking with libm) and runs in a bash script.

Mildly ungolfed:

f(double *z){
    double u=2**z,l=0,g=u,p=0;

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.