Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It only takes a minute to sign up.
Sign up to join this community
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
FindInstance[#==x!&&x>0,x]&
Output looks like {{x -> n}}, where n is the solution, which may not be allowed.
.I.!
A program that takes input of a number and prints the result.
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
1`1e-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!
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)
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.
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\$
Jan 29, 2017 at 17:46
@(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.
((]-(-~^.@!)%[:^.@!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.
((]-(-~^.@!)%[:^.@!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
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.
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.
Log[y=#]+1//.x_->x-(x!-y)/Gamma'[x+1]&
Ḋ!Æl_®
ȷİ‘×;µ!ÆlI÷I÷@Ç_@ḊḢ
Æl©‘ÇÐĿ
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)).
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.
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!)
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;
}
}
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\$