# Inverse Pi function

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

Examples:

120 -> 5.0000
10 -> 3.39008
3.14 -> 2.44815
2017 -> 6.53847
1.5 -> 1.66277

• 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. – orlp Jan 27 '17 at 19:12
• Wellp, that series expansion is enough to scare me... i.imgur.com/ttgzDSJ.gif – Magic Octopus Urn Jan 27 '17 at 19:22
• 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. – Greg Martin Jan 27 '17 at 19:28
• @GregMartin Added as well. – orlp Jan 27 '17 at 19:29
• 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. – Ruslan Jan 28 '17 at 17:23

# Mathematica, 1715 27 bytes

FindInstance[#==x!&&x>0,x]&


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

# Pyth, 4 bytes

.I.!


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

11e-5+tQYgG<


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!

### Explanation

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

NSolve[Gamma(x+1)=A1,x=1]


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: How it works

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

• Is it guaranteed to return a positive number, and does it work for 1.5, which has broken several answers? – Pavel Jan 28 '17 at 3:39
• @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. – TheBikingViking Jan 29 '17 at 17:46

# MATLAB, 59 bytes

@(x)fminsearch(@(t)(gamma(t+1)-x)^2,1,optimset('TolF',eps))


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 =
@(x)fminsearch(@(t)(gamma(t+1)-x)^2,1,optimset('TolF',eps))
>> f = ans; format long; f(120), f(10), f(3.14), f(2017)
ans =
5.000000000000008
ans =
3.390077650547032
ans =
2.448151165246967
ans =
6.538472664321318


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

# J, 86 33 bytes

((]-(-~^.@!)%[:^.@!D.1])^:_>:)@^.


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.

## Explanation

((]-(-~^.@!)%[:^.@!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[#-x!,{x,1}]&


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

FixedPoint[#-(#!-y)/Gamma'[#+1]&,Log[y=#]+1]&


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

Log[y=#]+1//.x_->x-(x!-y)/Gamma'[x+1]& # Jelly, 34 bytes

Ḋ!Æl_®
ȷİ‘×;µ!ÆlI÷I÷@Ç_@ḊḢ
Æl©‘ÇÐĿ


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

• I get an error with input 1.5 – Luis Mendo Jan 28 '17 at 0:35
• @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. – miles Jan 28 '17 at 0:53

## PARI/GP, 30 bytes

x->solve(t=1,x+1,gamma(t+1)-x)


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.

NArgMin[{(#-x!)^2,x>0},x]&


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;
for(;log(fabs(g-p))>-14;){
p=g;
g=(u+l)/2;
u=tgamma(g+1)>*z?g:(l=g,u);*z=g;
}
}