Gamma function: Difference between revisions
Added Easylang
(Added Easylang) |
|||
(35 intermediate revisions by 23 users not shown) | |||
Line 21:
{{trans|Python}}
<
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
Line 40:
L(i) 1..10
print(‘#.14’.format(gamma(i / 3.0)))</
{{out}}
Line 58:
=={{header|360 Assembly}}==
For maximum compatibility, this program uses only the basic instruction set.
<
USING GAMMAT,R13
SAVEAR B STM-SAVEAR(R15)
Line 216:
EQUREGS
EQUREGS REGS=FPR
END GAMMAT</
{{out}}
<pre style="height:20ex">
Line 256:
The coefficients are taken from ''Mathematical functions and their approximations''
by [[wp:Yudell Luke|Yudell L. Luke]].
<
A : constant array (0..29) of Long_Float :=
( 1.00000_00000_00000_00000,
Line 296:
end loop;
return 1.0 / Sum;
end Gamma;</
Test program:
<
with Gamma;
Line 306:
Put_Line (Long_Float'Image (Gamma (Long_Float (I) / 3.0)));
end loop;
end Test_Gamma;</
{{Out}}
<pre>
Line 329:
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}}
<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release 1.8.8d.fc9.i386 - formatted transput is missing, gamma is an extension, and various LONG vs INT operators missing}} -->
<
[]LONG REAL p = ( LONG 0.99999 99999 99809 93,
LONG 676.52036 81218 851,
Line 459:
FOR x FROM 10 BY 10 TO 70 DO
printf((sample exp fmt, x, sample(x)))
OD</
{{out}}
<pre>
Line 493:
</pre>
=={{header|
<syntaxhighlight lang="rebol">A: @[
1.00000000000000000000 0.57721566490153286061 neg 0.65587807152025388108
neg 0.04200263503409523553 0.16653861138229148950 neg 0.04219773455554433675
neg 0.00962197152787697356 0.00721894324666309954 neg 0.00116516759185906511
neg 0.00021524167411495097 0.00012805028238811619 neg 0.00002013485478078824
neg 0.00000125049348214267 0.00000113302723198170 neg 0.00000020563384169776
0.00000000611609510448 0.00000000500200764447 neg 0.00000000118127457049
0.00000000010434267117 0.00000000000778226344 neg 0.00000000000369680562
0.00000000000051003703 neg 0.00000000000002058326 neg 0.00000000000000534812
0.00000000000000122678 neg 0.00000000000000011813 0.00000000000000000119
0.00000000000000000141 neg 0.00000000000000000023 0.00000000000000000002
]
ourGamma: function [x][
y: x - 1
result: (result*y) + get A n
result: 1 // result
return result
]
loop 1..10 'z [
v2: gamma z // 3
print [
pad (to :string z)++" =>" 10
pad (to :string v1)++" ~" 30
pad (to :string v2)++" :" 30
pad (to :string v1-v2) 30
]
]</syntaxhighlight>
{{out}}
<pre> 1 => 2.678938534707748 ~ 2.678938534707748 : 4.440892098500626e-16
2 => 1.3541179394264 ~ 1.3541179394264 : 0.0
3 => 1.0 ~ 1.0 : 0.0
4 => 0.8929795115692493 ~ 0.8929795115692493 : 0.0
5 => 0.9027452929509336 ~ 0.9027452929509336 : 0.0
6 => 1.0 ~ 1.0 : 0.0
7 => 1.190639348758999 ~ 1.190639348758999 : 0.0
8 => 1.504575488251335 ~ 1.504575488251556 : -2.204902926905561e-13
9 => 1.999999999908069 ~ 2.0 : -9.193090733106146e-11
10 => 2.778158462440097 ~ 2.778158480437665 : -1.799756743636749e-08</pre>
=={{header|AutoHotkey}}==
{{AutoHotkey case}}
Source: [http://www.autohotkey.com/forum/topic44657.html AutoHotkey forum] by Laszlo
<
Here is the upper incomplete Gamma function. Omitting or setting
the second parameter to 0 we get the (complete) Gamma function.
Line 607 ⟶ 625:
1.000000000e+000 1.190639348e+000 1.504575489e+000 2.000000000e+000 2.778158479e+000
1.386831185e+080 1.711224524e+098 8.946182131e+116 1.650795516e+136 9.332621544e+155
*/</
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# syntax: GAWK -f GAMMA_FUNCTION.AWK
BEGIN {
e = (1+1/100000)^100000
pi = atan2(0,-1)
print("X Stirling")
for (i=1; i<=20; i++) {
Line 629 ⟶ 647:
return exp(b*log(a))
}
</syntaxhighlight>
{{out}}
<pre>
Line 655 ⟶ 673:
</pre>
=={{header|
==={{header|ANSI BASIC}}===
{{trans|BBC BASIC}} - Lanczos method.
{{works with|Decimal BASIC}}
<syntaxhighlight lang="basic">100 DECLARE EXTERNAL FUNCTION FNlngamma
110
120 DEF FNgamma(z) = EXP(FNlngamma(z))
130
140 FOR x = 0.1 TO 2.05 STEP 0.1
150 PRINT USING$("#.#",x), USING$("##.############", FNgamma(x))
160 NEXT x
170 END
180
190 EXTERNAL FUNCTION FNlngamma(z)
200 DIM lz(0 TO 6)
210 RESTORE
220 MAT READ lz
230 DATA 1.000000000190015, 76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.0012086509738662, -0.000005395239385
240 IF z < 0.5 THEN
250 LET FNlngamma = LOG(PI / SIN(PI * z)) - FNlngamma(1.0 - z)
260 EXIT FUNCTION
270 END IF
280 LET z = z - 1.0
290 LET b = z + 5.5
300 LET a = lz(0)
310 FOR i = 1 TO 6
320 LET a = a + lz(i) / (z + i)
330 NEXT i
340 LET FNlngamma = (LOG(SQR(2*PI)) + LOG(a) - b) + LOG(b) * (z+0.5)
350 END FUNCTION</syntaxhighlight>
{{out}}
<pre>
.1 9.513507698670
.2 4.590843712000
.3 2.991568987689
.4 2.218159543760
.5 1.772453850902
.6 1.489192248811
.7 1.298055332647
.8 1.164229713725
.9 1.068628702119
1.0 1.000000000000
1.1 .951350769867
1.2 .918168742400
1.3 .897470696306
1.4 .887263817503
1.5 .886226925453
1.6 .893515349288
1.7 .908638732853
1.8 .931383770980
1.9 .961765831907
2.0 1.000000000000
</pre>
==={{header|BASIC256}}===
{{trans|FreeBASIC}} - Stirling method.
{{trans|Phix}} - Lanczos method.
<syntaxhighlight lang="freebasic">print " x Stirling Lanczos"
print
for i = 1 to 20
d = i / 10.0
print d;
print chr(9); ljust(string(gammaStirling(d)), 13, "0");
print chr(9); ljust(string(gammaLanczos(d)), 13, "0")
next i
end
function gammaStirling (x)
e = exp(1) # e is not predefined in BASIC256
return sqr(2.0 * pi / x) * ((x / e) ^ x)
end function
function gammaLanczos (x)
dim p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7}
g = 7
if x < 0.5 then return pi / (sin(pi * x) * gammaLanczos(1-x))
x -= 1
a = p[0]
t = x + g + 0.5
for i = 1 to 8
a += p[i] / (x + i)
next i
return sqr(2.0 * pi) * (t ^ (x + 0.5)) * exp(-t) * a
end function</syntaxhighlight>
==={{header|BBC BASIC}}===
{{works with|BBC BASIC for Windows}}
Uses the Lanczos approximation.
<
INSTALL @lib$+"FNUSING"
Line 680 ⟶ 786:
a += lz(i%) / (z + i%)
NEXT
= (LNSQR(2*PI) + LN(a) - b) + LN(b) * (z+0.5)</
'''Output:'''
<pre>
Line 710 ⟶ 816:
This implements [[wp:Stirling's approximation|Stirling's approximation]] and [[wp:Spouge's approximation|Spouge's approximation]].
<
#include <stdlib.h>
#include <math.h>
Line 761 ⟶ 867:
}
return 0;
}</
{{out}}
<pre> Stirling Spouge GSL libm
Line 778 ⟶ 884:
=={{header|C sharp}}==
This is just rewritten from the Wikipedia Lanczos article. Works with complex numbers as well as reals.
<
using System.Numerics;
Line 805 ⟶ 911:
}
}
</syntaxhighlight>
=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <math.h>
#include <numbers>
#include <stdio.h>
#include <vector>
// Calculate the coefficients used by Spouge's approximation (based on the C
// implemetation)
std::vector<double> CalculateCoefficients(int numCoeff)
{
std::vector<double> c(numCoeff);
double k1_factrl = 1.0;
c[0] = sqrt(2.0 * std::numbers::pi);
for(size_t k=1; k < numCoeff; k++)
{
c[k] = exp(numCoeff-k) * pow(numCoeff-k, k-0.5) / k1_factrl;
k1_factrl *= -(double)k;
}
return c;
}
// The Spouge approximation
double Gamma(const std::vector<double>& coeffs, double x)
{
const size_t numCoeff = coeffs.size();
double accm = coeffs[0];
for(size_t k=1; k < numCoeff; k++)
{
accm += coeffs[k] / ( x + k );
}
accm *= exp(-(x+numCoeff)) * pow(x+numCoeff, x+0.5);
return accm/x;
}
int main()
{
// estimate the gamma function with 1, 4, and 10 coefficients
const auto coeff1 = CalculateCoefficients(1);
const auto coeff4 = CalculateCoefficients(4);
const auto coeff10 = CalculateCoefficients(10);
const auto inputs = std::vector<double>{
0.001, 0.01, 0.1, 0.5, 1.0,
1.461632145, // minimum of the gamma function
2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100,
150 // causes overflow for this implemetation
};
printf("%16s%16s%16s%16s%16s\n", "gamma( x ) =", "Spouge 1", "Spouge 4", "Spouge 10", "built-in");
for(auto x : inputs)
{
printf("gamma(%7.3f) = %16.10g %16.10g %16.10g %16.10g\n",
x,
Gamma(coeff1, x),
Gamma(coeff4, x),
Gamma(coeff10, x),
std::tgamma(x)); // built-in gamma function
}
}
</syntaxhighlight>
{{out}}
<pre>
gamma( x ) = Spouge 1 Spouge 4 Spouge 10 built-in
gamma( 0.001) = 921.6767466 999.4237321 999.4237725 999.4237725
gamma( 0.010) = 91.76063453 99.43258106 99.43258512 99.43258512
gamma( 0.100) = 8.834899532 9.513507269 9.513507699 9.513507699
gamma( 0.500) = 1.677913105 1.772453737 1.772453851 1.772453851
gamma( 1.000) = 0.9595021757 0.9999999124 1 1
gamma( 1.462) = 0.8562774501 0.8856030992 0.8856031944 0.8856031944
gamma( 2.000) = 0.9727015986 0.9999998717 1 1
gamma( 2.500) = 1.298145499 1.329340195 1.329340388 1.329340388
gamma( 3.000) = 1.958847928 1.999999681 2 2
gamma( 4.000) = 5.900958398 5.999998905 6 6
gamma( 5.000) = 23.66927282 23.99999523 24 24
gamma( 6.000) = 118.5808531 119.9999749 120 120
gamma( 7.000) = 712.5427759 719.9998442 720 720
gamma( 8.000) = 4993.567678 5039.99889 5040 5040
gamma( 9.000) = 39985.50687 40319.99104 40320 40320
gamma( 10.000) = 360142.0459 362879.9193 362880 362880
gamma( 50.000) = 6.072887637e+62 6.082817933e+62 6.08281864e+62 6.08281864e+62
gamma(100.000) = 9.324924563e+155 9.332620912e+155 9.332621544e+155 9.332621544e+155
gamma(150.000) = inf inf inf 3.808922638e+260
</pre>
=={{header|Clojure}}==
<
"Returns Gamma(z + 1 = number) using Lanczos approximation."
[number]
Line 822 ⟶ 1,012:
(Math/exp (- (+ n 7 0.5)))
(+ (first c)
(apply + (map-indexed #(/ %2 (+ n %1 1)) (next c))))))))</
{{out}}
<
<pre>
0.1 9.5135
Line 859 ⟶ 1,049:
=={{header|Common Lisp}}==
<
(defconstant tcoeff
'( 1.00000000000000000000 0.57721566490153286061 -0.65587807152025388108
Line 884 ⟶ 1,074:
(loop for i from 1 to 10
do (
format t "~12,10f~%" (gamma (/ i 3.0))))</
{{out|Produces}}
<pre>
Line 902 ⟶ 1,092:
====Taylor Series | Lanczos Method | Builtin Function====
{{trans|Taylor Series from Ruby; Lanczos Method from C#}}
<
def a
[ 1.00000_00000_00000_00000, 0.57721_56649_01532_86061, -0.65587_80715_20253_88108,
Line 944 ⟶ 1,134:
puts " Taylor Series Lanczos Method Builtin Function"
(1..27).each { |i| n = i/3.0; puts "gamma(%.2f) = %.14e %.14e %.14e" % [n, taylor_gamma(n), lanczos_gamma(n), Math.gamma(n)] }
</syntaxhighlight>
{{out}}
<pre>
Line 978 ⟶ 1,168:
=={{header|D}}==
<
real taylorGamma(in real x) pure nothrow @safe @nogc {
Line 1,039 ⟶ 1,229:
x.taylorGamma, x.lanczosGamma, x.gamma);
}
}</
{{out}}
<pre>0.333333: 2.6789385347077476335e+00 2.6789385347077470551e+00 2.6789385347077476339e+00
Line 1,051 ⟶ 1,241:
3.000000: 1.9999999999992207405e+00 2.0000000000000015575e+00 2.0000000000000000000e+00
3.333333: 2.7781584802531739378e+00 2.7781584804376666336e+00 2.7781584804376642124e+00</pre>
=={{header|Delphi}}==
{{libheader| System.SysUtils}}
{{libheader| System.Math}}
{{Trans|Go}}
<syntaxhighlight lang="delphi">
program Gamma_function;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
System.Math;
function lanczos7(z: double): Double;
begin
var t := z + 6.5;
var x := 0.99999999999980993 + 676.5203681218851 / z - 1259.1392167224028 / (z
+ 1) + 771.32342877765313 / (z + 2) - 176.61502916214059 / (z + 3) +
12.507343278686905 / (z + 4) - 0.13857109526572012 / (z + 5) +
9.9843695780195716e-6 / (z + 6) + 1.5056327351493116e-7 / (z + 7);
Result := Sqrt(2) * Sqrt(pi) * Power(t, z - 0.5) * exp(-t) * x;
end;
begin
var xs: TArray<double> := [-0.5, 0.1, 0.5, 1, 1.5, 2, 3, 10, 140, 170];
writeln(' x Lanczos7');
for var x in xs do
writeln(format('%5.1f %24.16g', [x, lanczos7(x)]));
readln;
end.</syntaxhighlight>
{{out}}
<pre> x Lanczos7
-0,5 -3,544907701811089
0,1 9,513507698668747
0,5 1,772453850905517
1,0 1
1,5 0,8862269254527583
2,0 1
3,0 2,000000000000002
10,0 362880,0000000007
140,0 9,615723196940235E238
170,0 4,269068009004271E304</pre>
=={{header|EasyLang}}==
{{trans|AWK}}
<syntaxhighlight>
e = 2.718281828459
func stirling x .
return sqrt (2 * pi / x) * pow (x / e) x
.
print " X Stirling"
for i to 20
d = i / 10
numfmt 2 4
write d & " "
numfmt 3 4
print stirling d
.
</syntaxhighlight>
=={{header|Elixir}}==
{{trans|Ruby}}
<
@a [ 1.00000_00000_00000_00000, 0.57721_56649_01532_86061, -0.65587_80715_20253_88108,
-0.04200_26350_34095_23553, 0.16653_86113_82291_48950, -0.04219_77345_55544_33675,
Line 1,074 ⟶ 1,325:
Enum.each(Enum.map(1..10, &(&1/3)), fn x ->
:io.format "~f ~18.15f~n", [x, Gamma.taylor(x)]
end)</
{{out}}
Line 1,094 ⟶ 1,345:
Solved using the Lanczos Coefficients described in Numerical Recipes (Press et al)
<syntaxhighlight lang="f sharp">
open System
Line 1,111 ⟶ 1,362:
seq { for i in 1 .. 10 do yield ((double)i*10.0) } |> Seq.iter ( fun v -> System.Console.WriteLine("{0} : {1}", v, gamma v ) )
</syntaxhighlight>
{{out}}
Line 1,145 ⟶ 1,396:
90 : 1.65079551625067E+136
100 : 9.33262154536104E+155
</pre>
{{Trans|C#}}
The C# version can be translated to F# to support complex numbers:
<syntaxhighlight lang="f sharp">
open System.Numerics
open System
let rec gamma (z: Complex) =
let mutable z = z
let lanczosCoefficients = [| 676.520368121885; -1259.1392167224; 771.323428777653; -176.615029162141; 12.5073432786869; -0.13857109526572; 9.98436957801957E-06; 1.50563273514931E-07 |]
if z.Real < 0.5 then
Math.PI / (sin (Math.PI * z) * gamma (1.0 - z))
else
let mutable x = Complex.One
z <- z - 1.0
for i = 0 to lanczosCoefficients.Length - 1 do
x <- x + lanczosCoefficients.[i] / (z + Complex(i, 0) + 1.0)
let t = z + float lanczosCoefficients.Length - 0.5
sqrt (2.0 * Math.PI) * (t ** (z + 0.5)) * exp (-t) * x
Seq.iter (fun i -> printfn "Gamma(%f) = %A" i (gamma (Complex(i, 0)))) [ 0 .. 100 ]
Seq.iter2 (fun i j -> printfn "Gamma(%f + i%f) = %A" i j (gamma (Complex(i, j)))) [ 0 .. 100 ] [ 0 .. 100 ]
</syntaxhighlight>
{{out}}
<pre>
Gamma(0.000000) = (NaN, NaN)
Gamma(1.000000) = (1,0000000000000049, 0)
Gamma(2.000000) = (1,0000000000000115, 0)
Gamma(3.000000) = (2,0000000000000386, 0)
Gamma(4.000000) = (6,000000000000169, 0)
Gamma(5.000000) = (24,00000000000084, 0)
Gamma(6.000000) = (120,00000000000514, 0)
Gamma(7.000000) = (720,0000000000364, 0)
Gamma(8.000000) = (5040,000000000289, 0)
Gamma(9.000000) = (40320,000000002554, 0)
Gamma(10.000000) = (362880,00000002526, 0)
Gamma(11.000000) = (3628800,0000002664, 0)
Gamma(12.000000) = (39916800,000003114, 0)
Gamma(13.000000) = (479001600,00004023, 0)
Gamma(14.000000) = (6227020800,000546, 0)
Gamma(15.000000) = (87178291200,00801, 0)
Gamma(16.000000) = (1307674368000,1267, 0)
Gamma(17.000000) = (20922789888002,08, 0)
Gamma(18.000000) = (355687428096034,7, 0)
Gamma(19.000000) = (6402373705728658, 0)
Gamma(20.000000) = (1,2164510040884514E+17, 0)
Gamma(21.000000) = (2,4329020081769083E+18, 0)
Gamma(22.000000) = (5,109094217171516E+19, 0)
Gamma(23.000000) = (1,1240007277777331E+21, 0)
Gamma(24.000000) = (2,5852016738887927E+22, 0)
Gamma(25.000000) = (6,20448401733312E+23, 0)
Gamma(26.000000) = (1,5511210043332811E+25, 0)
Gamma(27.000000) = (4,032914611266542E+26, 0)
Gamma(28.000000) = (1,0888869450419632E+28, 0)
Gamma(29.000000) = (3,048883446117507E+29, 0)
Gamma(30.000000) = (8,841761993740793E+30, 0)
Gamma(31.000000) = (2,652528598122235E+32, 0)
Gamma(32.000000) = (8,222838654178949E+33, 0)
Gamma(33.000000) = (2,631308369337265E+35, 0)
Gamma(34.000000) = (8,683317618812963E+36, 0)
Gamma(35.000000) = (2,9523279903964145E+38, 0)
Gamma(36.000000) = (1,0333147966387422E+40, 0)
Gamma(37.000000) = (3,719933267899472E+41, 0)
Gamma(38.000000) = (1,3763753091228064E+43, 0)
Gamma(39.000000) = (5,230226174666675E+44, 0)
Gamma(40.000000) = (2,0397882081200028E+46, 0)
Gamma(41.000000) = (8,159152832480012E+47, 0)
Gamma(42.000000) = (3,3452526613168034E+49, 0)
Gamma(43.000000) = (1,4050061177530564E+51, 0)
Gamma(44.000000) = (6,041526306338149E+52, 0)
Gamma(45.000000) = (2,658271574788788E+54, 0)
Gamma(46.000000) = (1,1962222086549537E+56, 0)
Gamma(47.000000) = (5,502622159812779E+57, 0)
Gamma(48.000000) = (2,586232415112008E+59, 0)
Gamma(49.000000) = (1,2413915592537631E+61, 0)
Gamma(50.000000) = (6,082818640343433E+62, 0)
Gamma(51.000000) = (3,04140932017172E+64, 0)
Gamma(52.000000) = (1,551118753287575E+66, 0)
Gamma(53.000000) = (8,065817517095389E+67, 0)
Gamma(54.000000) = (4,27488328406056E+69, 0)
Gamma(55.000000) = (2,308436973392699E+71, 0)
Gamma(56.000000) = (1,2696403353659833E+73, 0)
Gamma(57.000000) = (7,109985878049497E+74, 0)
Gamma(58.000000) = (4,0526919504882125E+76, 0)
Gamma(59.000000) = (2,350561331283167E+78, 0)
Gamma(60.000000) = (1,386831185457067E+80, 0)
Gamma(61.000000) = (8,3209871127424E+81, 0)
Gamma(62.000000) = (5,075802138772854E+83, 0)
Gamma(63.000000) = (3,1469973260391715E+85, 0)
Gamma(64.000000) = (1,9826083154046777E+87, 0)
Gamma(65.000000) = (1,2688693218589942E+89, 0)
Gamma(66.000000) = (8,247650592083449E+90, 0)
Gamma(67.000000) = (5,443449390775078E+92, 0)
Gamma(68.000000) = (3,647111091819299E+94, 0)
Gamma(69.000000) = (2,48003554243712E+96, 0)
Gamma(70.000000) = (1,711224524281613E+98, 0)
Gamma(71.000000) = (1,1978571669971308E+100, 0)
Gamma(72.000000) = (8,504785885679606E+101, 0)
Gamma(73.000000) = (6,123445837689312E+103, 0)
Gamma(74.000000) = (4,470115461513189E+105, 0)
Gamma(75.000000) = (3,307885441519755E+107, 0)
Gamma(76.000000) = (2,4809140811398187E+109, 0)
Gamma(77.000000) = (1,8854947016662648E+111, 0)
Gamma(78.000000) = (1,451830920283022E+113, 0)
Gamma(79.000000) = (1,1324281178207572E+115, 0)
Gamma(80.000000) = (8,946182130783977E+116, 0)
Gamma(81.000000) = (7,1569457046271725E+118, 0)
Gamma(82.000000) = (5,797126020748008E+120, 0)
Gamma(83.000000) = (4,753643337013366E+122, 0)
Gamma(84.000000) = (3,945523969721089E+124, 0)
Gamma(85.000000) = (3,31424013456571E+126, 0)
Gamma(86.000000) = (2,8171041143808564E+128, 0)
Gamma(87.000000) = (2,4227095383675335E+130, 0)
Gamma(88.000000) = (2,1077572983797526E+132, 0)
Gamma(89.000000) = (1,8548264225741817E+134, 0)
Gamma(90.000000) = (1,650795516091023E+136, 0)
Gamma(91.000000) = (1,4857159644819212E+138, 0)
Gamma(92.000000) = (1,3520015276785438E+140, 0)
Gamma(93.000000) = (1,2438414054642616E+142, 0)
Gamma(94.000000) = (1,1567725070817618E+144, 0)
Gamma(95.000000) = (1,0873661566568553E+146, 0)
Gamma(96.000000) = (1,032997848824012E+148, 0)
Gamma(97.000000) = (9,916779348710516E+149, 0)
Gamma(98.000000) = (9,619275968249195E+151, 0)
Gamma(99.000000) = (9,42689044888421E+153, 0)
Gamma(100.000000) = (9,33262154439535E+155, 0)
Gamma(0.000000 + i0.000000) = (NaN, NaN)
Gamma(1.000000 + i1.000000) = (0,49801566811835923, -0,15494982830180806)
Gamma(2.000000 + i2.000000) = (0,11229424234632254, 0,3236128855019324)
Gamma(3.000000 + i3.000000) = (-0,4401134076370088, -0,0636372431263299)
Gamma(4.000000 + i4.000000) = (0,7058649325913451, -0,49673908399741584)
Gamma(5.000000 + i5.000000) = (-0,9743952418053669, 2,0066898827226805)
Gamma(6.000000 + i6.000000) = (1,0560845455210948, -7,123931816061554)
Gamma(7.000000 + i7.000000) = (-0,26095668519941206, 27,88827411508434)
Gamma(8.000000 + i8.000000) = (1,8442848156317595, -125,96060801752867)
Gamma(9.000000 + i9.000000) = (-94,00399991734474, 643,3621714431141)
Gamma(10.000000 + i10.000000) = (1423,851941789479, -3496,081973308168)
Gamma(11.000000 + i11.000000) = (-16211,00700465313, 18168,810510285286)
Gamma(12.000000 + i12.000000) = (158471,8890918886, -68793,30331463458)
Gamma(13.000000 + i13.000000) = (-1329505,1052081874, -142199,12520863872)
Gamma(14.000000 + i14.000000) = (8576976,67312178, 7218722,503716219)
Gamma(15.000000 + i15.000000) = (-20768001,573587183, -99064686,32101583)
Gamma(16.000000 + i16.000000) = (-490395650,85195, 847486174,9207268)
Gamma(17.000000 + i17.000000) = (9782747798,66319, -2523864726,357996)
Gamma(18.000000 + i18.000000) = (-91408144092,80728, -62548333665,79536)
Gamma(19.000000 + i19.000000) = (80368797570,63837, 1283152922013,1064)
Gamma(20.000000 + i20.000000) = (12322153606702,379, -9813622771583,531)
Gamma(21.000000 + i21.000000) = (-191651224429571,5, -67416801166719,305)
Gamma(22.000000 + i22.000000) = (476610838765573,75, 2709614130691551,5)
Gamma(23.000000 + i23.000000) = (31282423285710508, -23340622982977492)
Gamma(24.000000 + i24.000000) = (-5,062346571412891E+17, -2,8075410996386413E+17)
Gamma(25.000000 + i25.000000) = (-1,1135374386470528E+18, 8,889271476011264E+18)
Gamma(26.000000 + i26.000000) = (1,4103207063357242E+20, -3,1111347801608966E+19)
Gamma(27.000000 + i27.000000) = (-1,20394153792971E+21, -2,100813378567903E+21)
Gamma(28.000000 + i28.000000) = (-2,9772583255514483E+22, 2,9845666640271078E+22)
Gamma(29.000000 + i29.000000) = (6,474322395366405E+23, 4,002110222682179E+23)
Gamma(30.000000 + i30.000000) = (4,982468347052982E+24, -1,3332730971666784E+25)
Gamma(31.000000 + i31.000000) = (-2,701233953257487E+26, -5,3335652308619844E+25)
Gamma(32.000000 + i32.000000) = (-3,5481927258667466E+26, 5,492417944693437E+27)
Gamma(33.000000 + i33.000000) = (1,13499763334942E+29, -3,936060260026878E+27)
Gamma(34.000000 + i34.000000) = (-2,497617380231244E+29, -2,4036706229026682E+30)
Gamma(35.000000 + i35.000000) = (-5,244047476964302E+31, 7,550247174479453E+30)
Gamma(36.000000 + i36.000000) = (1,8326953464053433E+32, 1,1815797667494789E+33)
Gamma(37.000000 + i37.000000) = (2,7496330950708185E+34, -3,7903032739968576E+33)
Gamma(38.000000 + i38.000000) = (-6,153702881069039E+34, -6,593476804299637E+35)
Gamma(39.000000 + i39.000000) = (-1,622188754583588E+37, 3,702408035066972E+35)
Gamma(40.000000 + i40.000000) = (-2,9787072201603815E+37, 4,069596487794187E+38)
Gamma(41.000000 + i41.000000) = (1,0327223226588694E+40, 2,028448840265614E+39)
Gamma(42.000000 + i42.000000) = (9,258591867044077E+40, -2,623834405446389E+41)
Gamma(43.000000 + i43.000000) = (-6,5825696143104E+42, -3,6674407207212435E+42)
Gamma(44.000000 + i44.000000) = (-1,3470021394042014E+44, 1,5970917674147318E+44)
Gamma(45.000000 + i45.000000) = (3,611191294816286E+45, 4,700641025433376E+45)
Gamma(46.000000 + i46.000000) = (1,572050028597055E+47, -6,978410499639601E+46)
Gamma(47.000000 + i47.000000) = (-8,064316457005957E+47, -5,037500759309252E+48)
Gamma(48.000000 + i48.000000) = (-1,5346852327090097E+50, -1,8750204416673013E+49)
Gamma(49.000000 + i49.000000) = (-1,963123049571723E+51, 4,3640542056764855E+51)
Gamma(50.000000 + i50.000000) = (1,112141672863102E+53, 1,0242389193853564E+53)
Gamma(51.000000 + i51.000000) = (4,3124175139534486E+54, -2,2723736598857867E+54)
Gamma(52.000000 + i52.000000) = (-1,9794169465380243E+55, -1,5907071359978385E+56)
Gamma(53.000000 + i53.000000) = (-5,198770735397539E+57, -1,364053254830618E+57)
Gamma(54.000000 + i54.000000) = (-1,120153853109903E+59, 1,4557034579325553E+59)
Gamma(55.000000 + i55.000000) = (3,0502206279504673E+60, 5,6213787579421754E+60)
Gamma(56.000000 + i56.000000) = (2,263901598439912E+62, -1,3869357395249425E+61)
Gamma(57.000000 + i57.000000) = (3,1337184675942048E+63, -7,566798379913671E+63)
Gamma(58.000000 + i58.000000) = (-1,9547779277051327E+65, -2,289059342042218E+65)
Gamma(59.000000 + i59.000000) = (-1,0990589957539216E+67, 2,43674092519427E+66)
Gamma(60.000000 + i60.000000) = (-1,2138821648921022E+68, 4,1070828497360523E+68)
Gamma(61.000000 + i61.000000) = (1,1429060333301634E+70, 1,199617402383565E+70)
Gamma(62.000000 + i62.000000) = (6,356301051435098E+71, -1,4385777048397157E+71)
Gamma(63.000000 + i63.000000) = (8,496122222382356E+72, -2,4629448359526944E+73)
Gamma(64.000000 + i64.000000) = (-6,555666234065589E+74, -8,308825832655404E+74)
Gamma(65.000000 + i65.000000) = (-4,35298630131942E+76, 3,566503620394996E+75)
Gamma(66.000000 + i66.000000) = (-9,093720755985274E+77, 1,5886817012950856E+78)
Gamma(67.000000 + i67.000000) = (3,276704637172729E+79, 7,067540192000104E+79)
Gamma(68.000000 + i68.000000) = (3,3000031716552424E+81, 6,606403155114497E+80)
Gamma(69.000000 + i69.000000) = (1,1027558271573558E+83, -9,8053446601259E+82)
Gamma(70.000000 + i70.000000) = (-4,322638823338879E+83, -6,551055369074648E+84)
Gamma(71.000000 + i71.000000) = (-2,4300476213174805E+86, -1,695901195606269E+86)
Gamma(72.000000 + i72.000000) = (-1,3066738432408389E+88, 3,647419995967721E+87)
Gamma(73.000000 + i73.000000) = (-2,631035116210543E+89, 5,722329713478048E+89)
Gamma(74.000000 + i74.000000) = (1,2168690154152985E+91, 2,7033309001464033E+91)
Gamma(75.000000 + i75.000000) = (1,3482397421660745E+93, 4,28037329878204E+92)
Gamma(76.000000 + i76.000000) = (5,937687180258459E+94, -3,3970662984923895E+94)
Gamma(77.000000 + i77.000000) = (7,877518886977203E+95, -3,258429077036584E+96)
Gamma(78.000000 + i78.000000) = (-8,893580405683355E+97, -1,4068641589865653E+98)
Gamma(79.000000 + i79.000000) = (-8,171120040126758E+99, -1,8182361788970335E+99)
Gamma(80.000000 + i80.000000) = (-3,620672062418723E+101, 2,252383956928435E+101)
Gamma(81.000000 + i81.000000) = (-5,470323981863164E+102, 2,130475015794634E+103)
Gamma(82.000000 + i82.000000) = (5,5003853678128614E+104, 1,0085770500872239E+105)
Gamma(83.000000 + i83.000000) = (5,740274728712029E+106, 1,986132310865518E+106)
Gamma(84.000000 + i84.000000) = (3,002634263016605E+108, -1,245707037705857E+108)
Gamma(85.000000 + i85.000000) = (7,865153067057636E+109, -1,575289481058232E+110)
Gamma(86.000000 + i86.000000) = (-2,2903303039804873E+111, -9,374401780068583E+111)
Gamma(87.000000 + i87.000000) = (-4,291880277570836E+113, -3,196185984412844E+113)
Gamma(88.000000 + i88.000000) = (-2,9995707183640408E+115, 1,1845722184191957E+114)
Gamma(89.000000 + i89.000000) = (-1,2943245318203042E+117, 1,1073023439932073E+117)
Gamma(90.000000 + i90.000000) = (-2,0007723042777158E+118, 9,568042990491017E+118)
Gamma(91.000000 + i91.000000) = (2,4147687416625443E+120, 5,132960890087842E+120)
Gamma(92.000000 + i92.000000) = (2,9270673985442265E+122, 1,5846423875074053E+122)
Gamma(93.000000 + i93.000000) = (1,9583971561168918E+124, -2,5166936071920716E+123)
Gamma(94.000000 + i94.000000) = (8,735240498823332E+125, -7,993043748530299E+125)
Gamma(95.000000 + i95.000000) = (1,6159775690644545E+127, -6,9922194609237035E+127)
Gamma(96.000000 + i96.000000) = (-1,569012926564151E+129, -4,10649302643008E+129)
Gamma(97.000000 + i97.000000) = (-2,2080981691082815E+131, -1,5902953016632018E+131)
Gamma(98.000000 + i98.000000) = (-1,6995881449793942E+133, -8,982249136857376E+131)
Gamma(99.000000 + i99.000000) = (-9,394719018633069E+134, 5,234766160326565E+134)
Gamma(100.000000 + i100.000000) = (-3,3597454530316526E+136, 5,986962556433683E+136)
</pre>
=={{header|Factor}}==
<
USING: picomath prettyprint ;
0.1 gamma . ! 9.513507698668723
2.0 gamma . ! 1.0
10. gamma . ! 362880.0</
=={{header|Forth}}==
Cristinel Mortici describes this method in Applied Mathematics Letters. "A substantial improvement of the Stirling formula". This algorithm is said to give about 7 good digits, but becomes more inaccurate close to zero. Therefore, a "shift" is performed to move the value returned into the more accurate domain.
<
: (mortici) ( f1 -- f2)
Line 1,170 ⟶ 1,654:
fdup (gamma-shift) s>f f+ (mortici) fswap
1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
;</
<pre>
0.1e gamma f. 9.51348888533932 ok
Line 1,178 ⟶ 1,662:
</pre>
This is a word, based on a formula of Ramanujan's famous "lost notebook", which was rediscovered in 1976. His formula contained a constant, which had a value between 1/100 and 1/30. In 2008, E.A. Karatsuba described the function, which determines the value of this constant. Since it contains the gamma function itself, it can't be used in a word calculating the gamma function, so here it is emulated by two symmetrical sigmoidals.
<
\ an approximation of the d(x) function
: ~d(x) ( f1 -- f2)
Line 1,202 ⟶ 1,686:
fdup (gamma-shift) 1- s>f f+ (ramanujan) fswap
1 s>f (gamma-shift) 0 do fover i s>f f+ f* loop fswap fdrop f/
;</
<pre>
0.1e gamma f. 9.51351721918848 ok
Line 1,214 ⟶ 1,698:
{{works with|Fortran|2008}}
{{works with|Fortran|95 with extensions}}
<
implicit none
Line 1,311 ⟶ 1,795:
end function lacz_gamma
end program ComputeGammaInt</
{{out}}
<pre> Simpson Lanczos Builtin
Line 1,327 ⟶ 1,811:
=={{header|FreeBASIC}}==
{{trans|Java}}
<
Const pi = 3.1415926535897932
Line 1,373 ⟶ 1,857:
Print
Print "Press any key to quit"
Sleep</
{{out}}
Line 1,402 ⟶ 1,886:
=={{header|Go}}==
<
import (
Line 1,428 ⟶ 1,912:
1.5056327351493116e-7/(z+7)
return math.Sqrt2 * math.SqrtPi * math.Pow(t, z-.5) * math.Exp(-t) * x
}</
{{out}}
<pre>
Line 1,446 ⟶ 1,930:
=={{header|Groovy}}==
{{trans|Ada}}
<
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
Line 1,460 ⟶ 1,944:
(1..10).each{ printf("% 1.9e\n", gamma(it / 3.0)) }
</syntaxhighlight>
{{out}}
<pre> 2.678938535e+00
Line 1,476 ⟶ 1,960:
Based on [http://www.haskell.org/haskellwiki/?title=Gamma_and_Beta_function&oldid=25546 HaskellWiki] ([http://www.haskell.org/haskellwiki/HaskellWiki:Copyrights compatible license]):
:The Gamma and Beta function as described in 'Numerical Recipes in C++', the approximation is taken from [Lanczos, C. 1964 SIAM Journal on Numerical Analysis, ser. B, vol. 1, pp. 86-96]
<
cof =
[ 76.18009172947146
Line 1,496 ⟶ 1,980:
main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]</
Or equivalently, as a point-free applicative expression:
<
cof :: [Double]
Line 1,520 ⟶ 2,004:
main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]</
{{Out}}
<pre>2.252712651734255
Line 1,537 ⟶ 2,021:
This works in Unicon. Changing the <tt>!10</tt> into <tt>(1 to 10)</tt> would enable it
to work in Icon.
<
every write(left(i := !10/10.0,5),gamma(.i))
end
Line 1,543 ⟶ 2,027:
procedure gamma(z) # Stirling's approximation
return (2*&pi/z)^0.5 * (z/&e)^z
end</
{{Out}}
Line 1,563 ⟶ 2,047:
=={{header|J}}==
This code shows the built-in method, which works for any value (positive, negative and complex numbers -- but note that negative integer arguments give infinite results).
<
The following direct coding of the task comes from the [[J:Essays/Stirling's%20Approximation|Stirling's approximation essay]] on the J wiki:
<
scorr =: 1 1r12 1r288 _139r51840 _571r2488320&p.@%
stirlg=: sbase * scorr</
Checking against <code>!@<:</code> we can see that this approximation loses accuracy for small arguments
<
10 362880 362880
3.14159 2.28803 2.28804
2.71828 1.56746 1.56747
1.5 0.886155 0.886227
1 0.999499 1</
(Column 1 is the argument, column 2 is the stirling approximation and column 3 uses the builtin support for gamma.)
=={{header|Java}}==
Implementation of Stirling's approximation and Lanczos approximation.
<
public double st_gamma(double x){
Line 1,611 ⟶ 2,095:
}
}
}</
{{out}}
<pre>
Line 1,639 ⟶ 2,123:
=={{header|JavaScript}}==
Implementation of Lanczos approximation.
<
var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
Line 1,658 ⟶ 2,142:
return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}</
=={{header|jq}}==
Line 1,664 ⟶ 2,148:
====Taylor Series====
{{trans|Ada}}
<
[
1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108, -0.04200263503409523553,
Line 1,679 ⟶ 2,163:
| reduce range(2; 1+$n) as $an
($a[$n-1]; (. * $y) + $a[$n - $an])
| 1 / . ;</
====Lanczos Approximation====
<
def gamma_by_lanczos:
def pow(x): if x == 0 then 1 elif x == 1 then . else x * log | exp end;
Line 1,696 ⟶ 2,180:
($p[0]; . + ($p[$i]/($x + $i) ))
* ((2*$pi) | sqrt) * ($t | pow($x+0.5)) * ((-$t)|exp)
end;</
====Stirling's Approximation====
<
def pow(x): if x == 0 then 1 elif x == 1 then . else x * log | exp end;
((1|atan) * 8) as $twopi
| . as $x
| (($twopi/$x) | sqrt) * ( ($x / (1|exp)) | pow($x));</
====Examples====
Stirling's method produces poor results, so to save space, the examples
contrast the Taylor series and Lanczos methods with built-in tgamma:
<
" i: gamma lanczos tgamma",
(range(1;11)
| . / 3.0
| "\(pad(18)): \(gamma|pad(18)) : \(gamma_by_lanczos|pad(18)) : \(tgamma)")</
{{Out}}
<
i: gamma lanczos tgamma
0.3333333333333333: 2.6789385347077483 : 2.6789385347077483 : 2.678938534707748
Line 1,724 ⟶ 2,208:
2.6666666666666665: 1.5045754882515399 : 1.5045754882515576 : 1.5045754882515558
3 : 1.9999999999939684 : 2.0000000000000013 : 2
3.3333333333333335: 2.778158479338573 : 2.778158480437665 : 2.7781584804376647</
=={{header|Jsish}}==
{{trans|Javascript}}
<
/* Gamma function, in Jsish, using the Lanczos approximation */
function gamma(x) {
Line 1,783 ⟶ 2,267:
5.5 +5.234278e+01
=!EXPECTEND!=
*/</
{{out}}
Line 1,818 ⟶ 2,302:
'''Built-in function''':
<syntaxhighlight lang
'''By adaptive Gauss-Kronrod integration''':
<
gammaquad(t::Float64) = first(quadgk(x -> x ^ (t - 1) * exp(-x), zero(t), Inf, reltol = 100eps(t)))
@show gammaquad(1.0)</
{{out}}
Line 1,831 ⟶ 2,315:
{{works with|Julia|1.0}}
'''Library function''':
<
gamma(1/2) - sqrt(pi)</
{{out}}
<pre>2.220446049250313e-16</pre>
=={{header|Koka}}==
Based on OCaml implementation
<syntaxhighlight lang="koka">
import std/num/float64
fun gamma-lanczos(x)
val g = 7.0
// Coefficients used by the GNU Scientific Library
val c = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
fun ag(z: float64, d: int)
if d == 0 then c[0].unjust + ag(z, 1)
elif d < 8 then c[d].unjust / (z + d.float64) + ag(z, d.inc)
else c[d].unjust / (z + d.float64)
fun gamma(z)
val z' = z - 1.0
val p = z' + g + 0.5
sqrt(2.0 * pi) * pow(p, (z' + 0.5)) * exp(0.0 - p) * ag(z', 0)
gamma(x)
val e = exp(1.0)
fun gamma-stirling(x)
sqrt(2.0 * pi / x) * pow(x / e, x)
fun gamma-stirling2(x')
// Extended Stirling method seen in Abramowitz and Stegun
val d = [1.0/12.0, 1.0/288.0, -139.0/51840.0, -571.0/2488320.0]
fun corr(z, x, n)
if n < d.length - 1 then d[n].unjust / x + corr(z, x*z, n.inc)
else d[n].unjust / x
fun gamma(z)
gamma-stirling(z)*(1.0 + corr(z, z, 0))
gamma(x')
fun mirror(gma, z)
if z > 0.5 then gma(z) else pi / sin(pi * z) / gma(1.0 - z)
fun main()
println("z\tLanczos\t\t\tStirling\t\tStirling2")
for(1, 20) fn(i)
val z = i.float64 / 10.0
println(z.show(1) ++ "\t" ++ mirror(gamma-lanczos, z).show ++ "\t" ++
mirror(gamma-stirling, z).show ++ "\t" ++ mirror(gamma-stirling2, z).show)
for(1, 7) fn(i)
val z = 10.0 * i.float64
println(z.show ++ "\t" ++ gamma-lanczos(z).show ++ "\t" ++
gamma-stirling(z).show ++ "\t" ++ gamma-stirling2(z).show)
</syntaxhighlight>
{{out}}
<pre>
z Lanczos Stirling Stirling2
0.1 9.5135076986687359 10.405084329555955 9.5210418318004439
0.2 4.5908437119988017 5.0739927535371763 4.596862295030256
0.3 2.9915689876875904 3.3503395433773222 2.9984402802949961
0.4 2.218159543757686 2.5270578096699556 2.2277588907113128
0.5 1.7724538509055157 2.0663656770612464 1.7883901437260497
0.6 1.4891922488128184 1.3071588574483559 1.4827753636029286
0.7 1.2980553326475577 1.1590532921139201 1.2950806801024195
0.8 1.1642297137253037 1.0533709684256085 1.1627054102439229
0.9 1.068628702119319 0.97706150787769541 1.0677830813298756
1.0 1.0000000000000002 0.92213700889578909 0.99949946853364036
1.1 0.95135076986687361 0.88348995316870382 0.95103799705518899
1.2 0.91816874239976076 0.85775533539659088 0.91796405783487933
1.3 0.89747069630627774 0.84267825944839203 0.8973312868034562
1.4 0.88726381750307537 0.8367445486370817 0.88716548484542823
1.5 0.88622692545275827 0.83895655252649626 0.88615538430170204
1.6 0.89351534928769061 0.8486932421525738 0.89346184003019224
1.7 0.90863873285329122 0.86562147179388405 0.90859770150945562
1.8 0.93138377098024272 0.8896396352879945 0.93135158986107858
1.9 0.96176583190738729 0.92084272189422933 0.96174006762796482
2.0 1.0000000000000002 0.95950217574449159 0.99997898067003532
10 362880.00000000105 359869.56187410367 362879.99717458693
20 1.2164510040883245e+17 1.2113934233805675e+17 1.2164510037907557e+17
30 8.841761993739658e+30 8.8172365307655063e+30 8.8417619934546387e+30
40 2.0397882081197221e+46 2.0355431612365591e+46 2.0397882081041343e+46
50 6.0828186403425409e+62 6.0726891878763362e+62 6.0828186403274418e+62
60 1.3868311854568534e+80 1.3849063858294502e+80 1.3868311854555093e+80
70 1.7112245242813438e+98 1.7091885781910795e+98 1.711224524280615e+98
</pre>
=={{header|Kotlin}}==
<
fun gammaStirling(x: Double): Double = Math.sqrt(2.0 * Math.PI / x) * Math.pow(x / Math.E, x)
Line 1,872 ⟶ 2,438:
println("%17.15f".format(gammaLanczos(d)))
}
}</
{{out}}
Line 1,899 ⟶ 2,465:
2.00 0.959502175744492 1.000000000000000
</pre>
=={{header|Lambdatalk}}==
Following Javascript, with Lanczos approximation.
<syntaxhighlight lang="Scheme">
{def gamma.p
{A.new 0.99999999999980993
676.5203681218851
-1259.1392167224028
771.32342877765313
-176.61502916214059
12.507343278686905
-0.13857109526572012
9.9843695780195716e-6
1.5056327351493116e-7
}}
-> gamma.p
{def gamma.rec
{lambda {:x :a :i}
{if {< :i {A.length {gamma.p}}}
then {gamma.rec :x
{+ :a {/ {A.get :i {gamma.p}} {+ :x :i}} }
{+ :i 1}}
else :a
}}}
-> gamma.rec
{def gamma
{lambda {:x}
{if {< :x 0.5}
then {/ {PI}
{* {sin {* {PI} :x}}
{gamma {- 1 :x}}}}
else {let { {:x {- :x 1}}
{:t {+ {- :x 1} 7 0.5}}
} {* {sqrt {* 2 {PI}}}
{pow :t {+ :x 0.5}}
{exp -:t}
{gamma.rec :x {A.first {gamma.p}} 1}}
}}}}
-> gamma
{S.map {lambda {:i}
{div} Γ(:i) = {gamma :i}}
{S.serie -5.5 5.5 0.5}}
Γ(-5.5) = 0.010912654781909836
Γ(-5) = -42755084646679.17
Γ(-4.5) = -0.06001960130050417
Γ(-4) = 267219279041745.34
Γ(-3.5) = 0.27008820585226917
Γ(-3) = -1425169488222640
Γ(-2.5) = -0.9453087204829418
Γ(-2) = 6413262697001885
Γ(-1.5) = 2.363271801207352
Γ(-1) = -25653050788007544
Γ(-0.5) = -3.5449077018110295
Γ(0) = Infinity
Γ(0.5) = 1.7724538509055159
Γ(1) = 0.9999999999999998
Γ(1.5) = 0.8862269254527586
Γ(2) = 1.0000000000000002
Γ(2.5) = 1.3293403881791384
Γ(3) = 2.000000000000001
Γ(3.5) = 3.3233509704478426
Γ(4) = 6.000000000000007
Γ(4.5) = 11.631728396567446
Γ(5) = 23.999999999999996
Γ(5.5) = 52.34277778455358
</syntaxhighlight>
=={{header|Limbo}}==
Line 1,905 ⟶ 2,543:
A fairly straightforward port of the Go code. (It could almost have been done with sed). A few small differences are in the use of a tuple as a return value for the builtin gamma function, and we import a few functions from the math library so that we don't have to qualify them.
<
include "sys.m"; sys: Sys;
Line 1,949 ⟶ 2,587:
return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;
}
</syntaxhighlight>
{{output}}
Line 1,967 ⟶ 2,605:
=={{header|Lua}}==
Uses the [[wp:Reciprocal gamma function]] to calculate small values.
<
function recigamma(z)
return z + gamma * z^2 + coeff * z^3 + quad * z^4 + qui * z^5 + set * z^6
Line 1,977 ⟶ 2,615:
else return (z - 1) * gammafunc(z-1)
end
end</
=={{header|M2000 Interpreter}}==
<syntaxhighlight lang="m2000 interpreter">
Module PrepareLambdaFunctions {
Const e = 2.7182818284590452@
Line 2,020 ⟶ 2,658:
Print $("")
clipboard doc$
</syntaxhighlight>
χ Stirling Lanczos
0.10 5.697187148977170 9.5135076986687024462927178610
Line 2,045 ⟶ 2,683:
=={{header|Maple}}==
Built-in method that accepts any value.
<
GAMMA(7*I);
M := Matrix(2, 3, 'fill' = -3.6);
MTM:-gamma(M);</
{{Out|Output}}
<pre>2027025*sqrt(Pi)*(1/256)
Line 2,054 ⟶ 2,692:
Matrix(2, 3, [[.2468571430, .2468571430, .2468571430], [.2468571430, .2468571430, .2468571430]])</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
This code shows the built-in method, which works for any value (positive, negative and complex numbers).
<syntaxhighlight lang
Output integers and half-integers (a space is multiplication in Mathematica):
<pre>
Line 2,103 ⟶ 2,741:
=={{header|Maxima}}==
<
gamma_coeff(n) := block([a: makelist(1, n)],
Line 2,123 ⟶ 2,761:
gamma_approx(12.3b0) - gamma(12.3b0);
/* -9.25224705314470500985141176997b-15 */</
=={{header|МК-61/52}}==
Line 2,136 ⟶ 2,774:
=={{header|Modula-3}}==
{{trans|Ada}}
<
FROM IO IMPORT Put;
Line 2,172 ⟶ 2,810:
Put(Extended(Taylor(FLOAT(i, EXTENDED) / 3.0X0), style := Style.Sci) & "\n");
END;
END Gamma.</
{{out}}
<pre>
Line 2,190 ⟶ 2,828:
{{trans|Ada}}
The algorithm is a translation of that from the Ada solution. We have added a comparison with the gamma function provided by the “math” module from Nim standard library (which is, in fact, the C “tgamma” function).
<syntaxhighlight lang
const A = [
1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108,
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
Line 2,203 ⟶ 2,843:
proc gamma(x: float): float =
let y = x
result =
for n in countdown(A.high
result = result * y +
result = 1
echo "Our gamma function Nim gamma function Difference"
echo "------------------ ------------------ ----------"
for i in 1..10:
let val2 = math.gamma(i.toFloat / 3)
echo &"{val1:18.16f} {val2:18.16f} {val1 - val2:11.4e}"</syntaxhighlight>
{{Out}}
<pre>Our gamma function Nim gamma function Difference
------------------ ------------------ ----------
2.6789385347077483 2.6789385347077479 4.4409e-16
1.3541179394264005 1.3541179394264005 0.0000e+00
1.0000000000000000 1.0000000000000000 0.0000e+00
0.8929795115692493 0.8929795115692493 0.0000e+00
0.9027452929509336 0.9027452929509336 0.0000e+00
1.0000000000000000 1.0000000000000000 0.0000e+00
1.1906393487589990 1.1906393487589990 0.0000e+00
1.5045754882515399 1.5045754882515558 -1.5987e-14
1.9999999999939684 2.0000000000000000 -6.0316e-12
2.7781584793385732 2.7781584804376647 -1.0991e-09</pre>
=={{header|OCaml}}==
<
let pi = 4. *. atan 1.
let sqrttwopi = sqrt (2. *. pi)
Line 2,286 ⟶ 2,933:
(Stirling.gamma z)
(Stirling2.gamma z)
done</
{{out}}
<pre>
Line 2,320 ⟶ 2,967:
=={{header|Octave}}==
<
p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \
771.32342877765313, -176.61502916214059, 12.507343278686905, \
Line 2,341 ⟶ 2,988:
for i = 1:10
printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));
endfor</
{{out}}
<pre>2.678939 2.678939
Line 2,357 ⟶ 3,004:
=={{header|Oforth}}==
<
[
Line 2,373 ⟶ 3,020:
2 Pi * sqrt *
t z 0.5 + powf *
t neg exp * ;</
{{out}}
Line 2,400 ⟶ 3,047:
=={{header|PARI/GP}}==
===Built-in===
<syntaxhighlight lang
===Double-exponential integration===
<code>[[+oo],k]</code> means that the function approaches <math>+\infty</math> as <math>\exp(-kx).</math>
<
===Romberg integration===
<
===Stirling approximation===
<
=={{header|Pascal}}==
A console application in Free Pascal, created with the Lazarus IDE.
Based on the algorithm for ln(Gamma(x)) (x > 0) in Press et al., Numerical Recipes, 3rd edition, pp. 256-7. For x >= 1/2, we simply take the exponential of their value; for x < 1/2 we calculate Gamma(1 - x) and use the reflection formula
Gamma(x)*Gamma(1 - x) = pi/sin(pi*x).
<syntaxhighlight lang="pascal">
program GammaTest;
{$mode objfpc}{$H+}
uses SysUtils;
function Gamma( x : extended) : extended;
const COF : array [0..14] of extended =
( 0.999999999999997092, // may as well include this in the array
57.1562356658629235,
-59.5979603554754912,
14.1360979747417471,
-0.491913816097620199,
0.339946499848118887e-4,
0.465236289270485756e-4,
-0.983744753048795646e-4,
0.158088703224912494e-3,
-0.210264441724104883e-3,
0.217439618115212643e-3,
-0.164318106536763890e-3,
0.844182239838527433e-4,
-0.261908384015814087e-4,
0.368991826595316234e-5);
const
K = 2.5066282746310005;
PI_OVER_K = PI / K;
var
j : integer;
tmp, w, ser : extended;
reflect : boolean;
begin
reflect := (x < 0.5);
if reflect then w := 1.0 - x else w := x;
tmp := w + 5.2421875;
tmp := (w + 0.5)*Ln(tmp) - tmp;
ser := COF[0];
for j := 1 to 14 do ser := ser + COF[j]/(w + j);
try
if reflect then
result := PI_OVER_K * w * Exp(-tmp) / (Sin(PI*x) * ser)
else
result := K * Exp(tmp) * ser / w;
except
raise SysUtils.Exception.CreateFmt(
'Gamma(%g) is undefined or out of floating-point range', [x]);
end;
end;
// Main routine for testing the Gamma function
var
x, k : extended;
begin
WriteLn( 'Is it seamless at x = 1/2 ?');
x := 0.49999999999999;
WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
x := 0.50000000000001;
WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
WriteLn( 'Test a few values:');
WriteLn( SysUtils.Format( 'Gamma(1) = %g', [Gamma(1)]));
WriteLn( SysUtils.Format( 'Gamma(2) = %g', [Gamma(2)]));
WriteLn( SysUtils.Format( 'Gamma(3) = %g', [Gamma(3)]));
WriteLn( SysUtils.Format( 'Gamma(10) = %g', [Gamma(10)]));
WriteLn( SysUtils.Format( 'Gamma(101) = %g', [Gamma(101)]));
WriteLn( ' 100! = 9.33262154439442E157');
WriteLn( SysUtils.Format( 'Gamma(1/2) = %g', [Gamma(0.5)]));
WriteLn( SysUtils.Format( 'Sqrt(pi) = %g', [Sqrt(PI)]));
WriteLn( SysUtils.Format( 'Gamma(-7/2) = %g', [Gamma(-3.5)]));
(*
Note here a bug or misfeature in Lazarus (doesn't occur in Delphi):
Putting (16.0/105.0)*Sqrt(PI) does not give the required precision.
We have to explicitly define the integers as extended floating-point.
*)
k := extended(16.0)/extended(105.0);
WriteLn( SysUtils.Format( ' (16/105)Sqrt(pi) = %g', [k*Sqrt(PI)]));
WriteLn( SysUtils.Format( 'Gamma(-9/2) = %g', [Gamma(-4.5)]));
k := extended(32.0)/extended(945.0);
WriteLn( SysUtils.Format( '-(32/945)Sqrt(pi) = %g', [-k*Sqrt(PI)]));
end.
</syntaxhighlight>
{{out}}
<pre>
Is it seamless at x = 1/2 ?
Gamma(0.49999999999999) = 1.77245385090555
Gamma(0.50000000000001) = 1.77245385090548
Test a few values:
Gamma(1) = 1
Gamma(2) = 1
Gamma(3) = 2
Gamma(10) = 362880
Gamma(101) = 9.33262154439452E157
100! = 9.33262154439442E157
Gamma(1/2) = 1.77245385090552
Sqrt(pi) = 1.77245385090552
Gamma(-7/2) = 0.270088205852269
(16/105)Sqrt(pi) = 0.270088205852269
Gamma(-9/2) = -0.0600196013005043
-(32/945)Sqrt(pi) = -0.0600196013005042
</pre>
=={{header|Perl}}==
<
use warnings;
use constant pi => 4*atan2(1, 1);
Line 2,482 ⟶ 3,234:
print join(' ', map { sprintf "%.12f", Gamma($_/3, $method) } 1 .. 10);
print "\n";
}</
{{out}}
<pre> MPFR: 2.678938534708 1.354117939426 1.000000000000 0.892979511569 0.902745292951 1.000000000000 1.190639348759 1.504575488252 2.000000000000 2.778158480438
Line 2,490 ⟶ 3,242:
=={{header|Phix}}==
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">12</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">spouge_gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">accm</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">accm</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">k1_factrl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #000080;font-style:italic;">-- (k - 1)!*(-1)^k with 0!==1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">12</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1.5</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">k1_factrl</span>
<span style="color: #000000;">k1_factrl</span> <span style="color: #0000FF;">*=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">12</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">c</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]/(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">*=</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">12</span><span style="color: #0000FF;">))*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">12</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- Gamma(z+1)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">accm</span><span style="color: #0000FF;">/</span><span style="color: #000000;">z</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">taylor_gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- (good for values between 0 and 1, apparently)</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">1.00000_00000_00000_00000</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.57721_56649_01532_86061</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.65587_80715_20253_88108</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.04200_26350_34095_23553</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.16653_86113_82291_48950</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.04219_77345_55544_33675</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00962_19715_27876_97356</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00721_89432_46663_09954</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00116_51675_91859_06511</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00021_52416_74114_95097</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00012_80502_82388_11619</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00002_01348_54780_78824</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_12504_93482_14267</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_11330_27231_98170</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_02056_33841_69776</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00061_16095_10448</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00050_02007_64447</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00011_81274_57049</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00001_04342_67117</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_07782_26344</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00000_03696_80562</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_00510_03703</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00000_00020_58326</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00000_00005_34812</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_00001_22678</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00000_00000_11813</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_00000_00119</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_00000_00141</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.00000_00000_00000_00023</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">0.00000_00000_00000_00002</span> <span style="color: #0000FF;">}</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">[$]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">*</span><span style="color: #000000;">y</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">n</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">s</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">lanczos_gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">z</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0.5</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">return</span> <span style="color: #004600;">PI</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">*</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">lanczos_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">z</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000080;font-style:italic;">-- use a lanczos approximation:</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.99999999999980993</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">z</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">6.5</span><span style="color: #0000FF;">;</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span> <span style="color: #000000;">676.5203681218851</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">1259.1392167224028</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">771.32342877765313</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">176.61502916214059</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">12.507343278686905</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">-</span><span style="color: #000000;">0.13857109526572012</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">9.9843695780195716e-6</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">1.5056327351493116e-7</span> <span style="color: #0000FF;">}</span>
<span style="color: #000000;">z</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">z</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">+</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">exp</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">x</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">sqPI</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">sq</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">zm</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span><span style="color: #0000FF;">=</span><span style="color: #008000;">"%19.16f"</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mul</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">zm</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">mul</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">spouge_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">taylor_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">lanczos_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)},</span>
<span style="color: #000000;">error</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_abs</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e</span><span style="color: #0000FF;">))</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">join</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #008000;">", "</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">:=</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">)&</span><span style="color: #008000;">", "</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">bdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">smallest</span><span style="color: #0000FF;">(</span><span style="color: #000000;">error</span><span style="color: #0000FF;">,</span><span style="color: #000000;">return_index</span><span style="color: #0000FF;">:=</span><span style="color: #004600;">true</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">best</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">bdx</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">bdx</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">mul</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> <span style="color: #000080;font-style:italic;">-- (potentially mark >1)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">best</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">22</span><span style="color: #0000FF;">-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..</span><span style="color: #000000;">i</span><span style="color: #0000FF;">*</span><span style="color: #000000;">22</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"*,"</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">es</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">e</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%5g: %s %s, %19.16f\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">es</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">*</span><span style="color: #000000;">p</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" z ------ spouge ----- ----- taylor ------ ----- lanczos ----- ---- expected ----- %19.16f\n"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">({{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">/</span><span style="color: #000000;">4</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">/</span><span style="color: #000000;">3</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">7</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">8</span><span style="color: #0000FF;">/</span><span style="color: #000000;">15</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">6</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">sq</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">0.001</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">999.4237725</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.15f"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">0.01</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">99.43258512</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.16f"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">0.1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">9.513507699</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.16f"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">362880</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.12f"</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">100</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">9.332621544e155</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.13g"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">machine_bits</span><span style="color: #0000FF;">()=</span><span style="color: #000000;">64</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">sq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">150</span><span style="color: #0000FF;">,</span><span style="color: #000000;">sqPI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">3.808922638e260</span><span style="color: #0000FF;">},</span><span style="color: #008000;">"%19.13g"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (fatal power overflow error on 32 bits)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<!--</syntaxhighlight>-->
{{out}}
The closest to the expected result for each z (row) is marked with a trailing asterisk.<br>
The final column is the value of PI (to 16dp) we would get from that best/starred result.
<pre>
z ------ spouge ----- ----- taylor ------ ----- lanczos ----- ---- expected ----- 3.1415926535897932
-1.5: 2.3632718012073547*, 2.3632718095606211, 2.3632718012073532, 2.3632718012073547, 3.1415926535897932
-0.5: -3.5449077018110320*, -3.5449077018110306, -3.5449077018110308, -3.5449077018110321, 3.1415926535897932
0.5: 1.7724538509055158, 1.7724538509055160*, 1.7724538509055166, 1.7724538509055160, 3.1415926535897932
1: 0.9999999999999998, 1.0000000000000000*, 1.0000000000000002, 1.0000000000000000, 3.1415926535897932
1.5: 0.8862269254527577, 0.8862269254527580*, 0.8862269254527583, 0.8862269254527580, 3.1415926535897932
2: 0.9999999999999994, 1.0000000000000000*, 1.0000000000000005, 1.0000000000000000, 3.1415926535897932
2.5: 1.3293403881791359, 1.3293403881791365*, 1.3293403881791379, 1.3293403881791370, 3.1415926535897906
3: 1.9999999999999978, 1.9999999999939679, 2.0000000000000016*, 2.0000000000000000, 3.1415926535897981
3.5: 3.3233509704478376, 3.3233509583896768, 3.3233509704478456*, 3.3233509704478426, 3.1415926535897990
4: 5.9999999999999884, 5.9999914100724727, 6.0000000000000063*, 6.0000000000000000, 3.1415926535897998
0.001: 999.423772484595421*, 999.423772484595404, 999.423772484595254, 999.423772500000000, 3.1415926534929476
0.01: 99.4325851191505990, 99.4325851191506035*, 99.4325851191505828, 99.4325851200000000, 3.1415926535361195
0.1: 9.5135076986687313, 9.5135076986687318*, 9.5135076986687300, 9.5135076990000000, 3.1415926533710076
10: 362879.999999996094, 0.000000029163, 362880.000000000725*, 362880.000000000000, 3.1415926535898058
100: 9.332621544394e+155, 7.510232292979e-39, 9.332621544394e+155*, 9.332621544e+155, 3.1415926538549222
150: 3.80892263763e+260*, 5.128102530869e-44, 3.80892263763e+260, 3.808922638e+260, 3.1415926529801383
</pre>
=== mpfr version ===
Above translated to mpfr,
{{libheader|Phix/mpfr}}
<!--<syntaxhighlight lang="phix">-->
<span style="color: #008080;">without</span> <span style="color: #008080;">javascript_semantics</span> <span style="color: #000080;font-style:italic;">-- (no mpfr_exp(), mpfr_gamma() in pwa/p2js)</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">dp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">30</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"%5s: %33s, %33s, %32s\n"</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.2"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (mpfr_get_fixed(maxlen), mpfr_gamma)</span>
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #7060A8;">mpfr_set_default_precision</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">87</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- 87 decimal places. </span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">mc</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_inits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">40</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">mpfr_spouge_gamma</span><span style="color: #0000FF;">(</span><span style="color: #004080;">mpfr</span> <span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">mpfr_cmp_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #000080;font-style:italic;">-- mc[1] := sqrt(2*PI)</span>
<span style="color: #7060A8;">mpfr_const_pi</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- k1_factrl = (k - 1)!*(-1)^k with 0!==1</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">k1_factrl</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">tmk</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000080;font-style:italic;">-- mc[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl</span>
<span style="color: #7060A8;">mpfr_set_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">tmk</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">k</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mpfr_exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">tmk</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_set_d</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1.5</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mpfr_pow</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">tmk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k1_factrl</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- k1_factrl *= -(k-1)</span>
<span style="color: #7060A8;">mpfr_mul_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k1_factrl</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k1_factrl</span><span style="color: #0000FF;">,-(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">accm</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000080;font-style:italic;">-- accm += mc[k]/(z+k-1)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">ck</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]),</span>
<span style="color: #000000;">zk</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">zk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">zk</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ck</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ck</span><span style="color: #0000FF;">,</span><span style="color: #000000;">zk</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ck</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- atom zc = z+length(mc)
-- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">ez</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">zh</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mc</span><span style="color: #0000FF;">))</span>
<span style="color: #7060A8;">mpfr_neg</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ez</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mpfr_exp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">ez</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ez</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">zh</span><span style="color: #0000FF;">,</span><span style="color: #000000;">zh</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mpfr_pow</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">zh</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ez</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- return accm/z</span>
<span style="color: #7060A8;">mpfr_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">accm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">accm</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">mPI</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">mqPI</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpfr_const_pi</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mPI</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">pistr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mPI</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">makempfr</span><span style="color: #0000FF;">(</span><span style="color: #004080;">object</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">if</span> <span style="color: #004080;">string</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #008000;">'/'</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #7060A8;">mpfr_div_si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">to_integer</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]))</span>
<span style="color: #008080;">elsif</span> <span style="color: #004080;">sequence</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #0000FF;">{</span><span style="color: #004080;">mpfr</span> <span style="color: #000000;">x1</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">object</span> <span style="color: #000000;">x2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #004080;">string</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">d</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">else</span>
<span style="color: #7060A8;">mpfr_div_d</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #000000;">x2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">else</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">mq</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">zm</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">=</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">mul</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">apply</span><span style="color: #0000FF;">(</span><span style="color: #000000;">zm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">makempfr</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mpfr_spouge_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">maxlen</span><span style="color: #0000FF;">:=</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">mpfr_gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">e</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">mpfr</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_init_set</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">mul</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">mpfr_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">zs</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">es</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">e</span><span style="color: #0000FF;">,</span><span style="color: #000000;">d</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">maxlen</span><span style="color: #0000FF;">:=</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">ps</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">zs</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">es</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ps</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" z %s %s %s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #7060A8;">pad</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" spouge "</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"BOTH"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">'-'</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">pad</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" expected "</span><span style="color: #0000FF;">,</span><span style="color: #000000;">dp</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"BOTH"</span><span style="color: #0000FF;">,</span><span style="color: #008000;">'-'</span><span style="color: #0000FF;">),</span><span style="color: #000000;">pistr</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">({{-</span><span style="color: #000000;">1.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">0.75</span><span style="color: #0000FF;">},{-</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">0.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}},{</span><span style="color: #000000;">1.5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}},{</span><span style="color: #000000;">2.5</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"4/3"</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">}},{</span><span style="color: #000000;">3.5</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"8/15"</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">}}},</span><span style="color: #000000;">mq</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #008000;">"1/1000"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"999.4237724845954661149822012996"</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">28</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #008000;">"1/100"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"99.43258511915060371353298887051"</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">29</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #008000;">"1/10"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"9.513507698668731836292487177265"</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">30</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">362880</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">25</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">100</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"9.332621544394415268169923885627e155"</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">mq</span><span style="color: #0000FF;">({</span><span style="color: #000000;">150</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">mqPI</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"3.80892263763056972698595524350735e260"</span><span style="color: #0000FF;">}},</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
z ------------ spouge ------------ ----------- expected ----------- 3.141592653589793238462643383279
-1.5: 2.363271801207354703064223311121, 2.363271801207354703064223311121, 3.141592653589793238462643383279
-0.5: -3.544907701811032054596334966e0, -3.544907701811032054596334966e0, 3.141592653589793238462643383279
0.5: 1.772453850905516027298167483341, 1.772453850905516027298167483341, 3.141592653589793238462643383279
1: 1.000000000000000000000000000000, 1, 3.141592653589793238462643383279
1.5: 0.886226925452758013649083741670, 0.886226925452758013649083741670, 3.141592653589793238462643383279
2: 1.000000000000000000000000000000, 1, 3.141592653589793238462643383279
2.5: 1.329340388179137020473625612505, 1.329340388179137020473625612505, 3.141592653589793238462643383279
3: 2.000000000000000000000000000000, 2, 3.141592653589793238462643383279
3.5: 3.323350970447842551184064031264, 3.323350970447842551184064031264, 3.141592653589793238462643383279
4: 6.000000000000000000000000000000, 6, 3.141592653589793238462643383279
0.001: 999.4237724845954661149822012996, 999.4237724845954661149822012996, 3.141592653589793238462643383279
0.009: 99.43258511915060371353298887051, 99.43258511915060371353298887051, 3.141592653589793238462643383279
0.1: 9.513507698668731836292487177265, 9.513507698668731836292487177265, 3.141592653589793238462643383279
10: 362880.0000000000000000000000000, 362880, 3.141592653589793238462643383279
100: 9.33262154439441526816992388e155, 9.33262154439441526816992388e155, 3.141592653589793238462643383279
150: 3.80892263763056972698595524e260, 3.80892263763056972698595524e260, 3.141592653589793238462643383279
</pre>
=={{header|Phixmonti}}==
<
-0.65587807152056 var coeff
-0.042002635033944 var quad
Line 2,707 ⟶ 3,554:
for
dup print " = " print gammafunc print nl
endfor</
=={{header|PicoLisp}}==
{{trans|Ada}}
<
(de *A
Line 2,730 ⟶ 3,577:
(for A (cdr *A)
(setq Sum (+ A (*/ Sum Y 1.0))) )
(*/ 1.0 1.0 Sum) ) )</
{{out}}
<pre>: (for I (range 1 10)
Line 2,746 ⟶ 3,593:
=={{header|PL/I}}==
<
test: procedure options (main);
Line 2,790 ⟶ 3,637:
end lanczos_gamma;
end test;</
{{out}}
<pre>
Line 2,808 ⟶ 3,655:
=={{header|PowerShell}}==
I would download the Math.NET Numerics dll(s). Documentation and download at: http://cyber-defense.sans.org/blog/2015/06/27/powershell-for-math-net-numerics/comment-page-1/
<syntaxhighlight lang="powershell">
Add-Type -Path "C:\Program Files (x86)\Math\MathNet.Numerics.3.12.0\lib\net40\MathNet.Numerics.dll"
1..20 | ForEach-Object {[MathNet.Numerics.SpecialFunctions]::Gamma($_ / 10)}
</syntaxhighlight>
{{Out}}
Line 2,836 ⟶ 3,683:
0.961765831907388
1
</pre>
=={{header|Prolog}}==
This version matches Wolfram Alpha to within a few digits at the end, so the last few digits are a bit off. There's an early check to stop evaluating coefficients once the desired accuracy is reached. Removing this check does not improve accuracy vs. Wolfram Alpha.
<syntaxhighlight lang="prolog">
gamma_coefficients(
[ 1.00000000000000000000000, 0.57721566490153286060651, -0.65587807152025388107701,
-0.04200263503409523552900, 0.16653861138229148950170, -0.04219773455554433674820,
-0.00962197152787697356211, 0.00721894324666309954239, -0.00116516759185906511211,
-0.00021524167411495097281, 0.00012805028238811618615, -0.00002013485478078823865,
-0.00000125049348214267065, 0.00000113302723198169588, -0.00000020563384169776071,
0.00000000611609510448141, 0.00000000500200764446922, -0.00000000118127457048702,
0.00000000010434267116911, 0.00000000000778226343990, -0.00000000000369680561864,
0.00000000000051003702874, -0.00000000000002058326053, -0.00000000000000534812253,
0.00000000000000122677862, -0.00000000000000011812593, 0.00000000000000000118669,
0.00000000000000000141238, -0.00000000000000000022987, 0.00000000000000000001714
]).
tolerance(1e-17).
gamma(X, _) :- X =< 0.0, !, fail.
gamma(X, Y) :-
X < 1.0, small_gamma(X, Y), !.
gamma(1, 1) :- !.
gamma(1.0, 1) :- !.
gamma(X, Y) :-
X1 is X - 1,
gamma(X1, Y1),
Y is X1 * Y1.
small_gamma(X, Y) :-
gamma_coefficients(Cs),
recip_gamma(X, 1.0, Cs, 1.0, 0.0, Y0),
Y is 1 / Y0.
recip_gamma(_, _, [], _, Y, Y) :- !.
recip_gamma(_, _, [], X0, X1, Y) :- tolerance(Tol), abs(X1 - X0) < Tol, Y = X1, !. % early exit
recip_gamma(X, PrevPow, [C|Cs], _, X1, Y) :-
Power is PrevPow * X,
X2 is X1 + C*Power,
recip_gamma(X, Power, Cs, X1, X2, Y).
</syntaxhighlight>
{{Out}}
<pre>
% see how close gamma(0.5) is to the square root of pi.
?- gamma(0.5,X), Y is sqrt(pi), Err is abs(X - Y).
X = 1.772453850905516,
Y = 1.7724538509055159,
Err = 2.220446049250313e-16.
?- gamma(1.5,X).
X = 0.886226925452758.
?- gamma(4.9,X).
X = 20.667385961857857.
?- gamma(5,X).
X = 24.
?- gamma(5.01,X).
X = 24.364473447872836.
?- gamma(6.9,X).
X = 597.4941281573107.
?- gamma(6.95,X).
X = 655.7662628554252.
?- gamma(7,X).
X = 720.
% 100!
?- gamma(101.0,X).
X = 9.33262154439441e+157.
% Note when passed integer, gamma(101) returns full big int precision
?- gamma(101,X).
X = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000.
?- gamma(100.98,X).
X = 8.510619261391532e+157.
</pre>
Line 2,843 ⟶ 3,771:
* Natural Logarithm of the Complete Gamma function
* Factorial function
<
; Complete Gamma function for x>0 and x<2 (approx)
; Extended outside this range via: Gamma(x+1) = x*Gamma(x)
Line 2,924 ⟶ 3,852:
Procedure Factorial(x) ; AKJ 01-May-10
ProcedureReturn Gamma(x+1)
EndProcedure</
;Examples
<
For i = 12 To 15
Debug StrD(i/3.0, 3)+" "+StrD(Gamma(i/3.0))
Line 2,933 ⟶ 3,861:
Debug "Ln(Gamma(5.0)) = "+StrD(GammLn(5.0), 16) ; Ln(24)
Debug ""
Debug "Factorial 6 = "+StrD(Factorial(6), 0) ; 72</
{{out}}
<pre>[Debug] Gamma():
Line 2,948 ⟶ 3,876:
===Procedural===
{{trans|Ada}}
<
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
Line 2,970 ⟶ 3,898:
for i in range(1,11):
print " %20.14e" % gamma(i/3.0)
</syntaxhighlight>
{{out}}
<pre> 2.67893853470775e+00
Line 2,986 ⟶ 3,914:
In terms of fold/reduce:
{{Works with|Python|3.7}}
<
from functools import reduce
Line 3,067 ⟶ 3,995:
# MAIN ---
if __name__ == '__main__':
main()</
{{Out}}
<pre> i -> gamma(i/3):
Line 3,085 ⟶ 4,013:
Lanczos' approximation is loosely converted from the Octave code.
{{trans|Octave}}
<
nemes <- function(z) sqrt(2*pi/z) * (exp(-1)*(z + (12*z - (10*z)^-1)^-1))^z
Line 3,107 ⟶ 4,035:
{
z <- z - 1
x <- p[1]
for (i in 1:8) {
x <- x+p[i+1]/(z+i)
}
tt <- z + g + 0.5
sqrt(2*pi) * tt^(z+0.5) * exp(-tt) * x
Line 3,134 ⟶ 4,065:
all.equal(as.complex(gamma(z)), lanczos(z)) # TRUE
all.equal(gamma(z), spouge(z)) # TRUE
data.frame(z=z, stirling=stirling(z), nemes=nemes(z), lanczos=lanczos(z), spouge=spouge(z), builtin=gamma(z))</
{{out}}
z stirling nemes lanczos spouge builtin
Line 3,149 ⟶ 4,080:
=={{header|Racket}}==
<
(define (gamma number)
(if (> 1/2 number)
Line 3,177 ⟶ 4,108:
; 1.1642297137253037
; 1.068628702119319
; 1.0)</
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku"
constant g = 9;
z < .5 ?? pi/ sin(pi * z) / Γ(1 - z) !!
Line 3,202 ⟶ 4,133:
}
say Γ($_) for 1/3, 2/3 ... 10/3;</
{{out}}
<pre>2.67893853470775
Line 3,222 ⟶ 4,153:
Note: The Taylor series isn't much good above values of <big>'''6½'''</big>.
<
/*The GAMMA function symbol is the Greek capital letter: Γ */
numeric digits 90 /*be able to handle extended precision.*/
Line 3,289 ⟶ 4,220:
sum=sum*xm + #.k
end /*k*/
return 1/sum</
{{out|output|text= when using the input of: <tt> 0.5 </tt>}}
<pre>
Line 3,313 ⟶ 4,244:
Many of the "normal" high-level mathematical functions aren't available in REXX, so some of them (RYO) are included here.
<
e=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138
numeric digits length(e) - length(.) /*use the number of decimal digits in E*/
Line 3,381 ⟶ 4,312:
do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/
numeric digits d; return g/1</
{{out|output|text= when using the default input:}}
<pre>
Line 3,397 ⟶ 4,328:
=={{header|Ring}}==
<
decimals(3)
gamma = 0.577
Line 3,416 ⟶ 4,347:
but fabs(z) <= 0.5 return 1 / recigamma(z)
else return (z - 1) * gammafunc(z-1) ok
</syntaxhighlight>
=={{header|RLaB}}==
Line 3,427 ⟶ 4,358:
:<math>RecGamma(a) = \frac{1}{\Gamma(a)}</math>;
:<math>Pochhammer(a,x) = \frac{\Gamma(a+x)}{\Gamma(x)}</math>.
=={{header|RPL}}==
{{trans|Ada}}
≪ 1 -
{ 1.00000000000000000000 0.57721566490153286061 -0.65587807152025388108
-0.04200263503409523553 0.16653861138229148950 -0.04219773455554433675
-0.00962197152787697356 0.00721894324666309954 -0.00116516759185906511
-0.00021524167411495097 0.00012805028238811619 -0.00002013485478078824
-0.00000125049348214267 0.00000113302723198170 -0.00000020563384169776
0.00000000611609510448 0.00000000500200764447 -0.00000000118127457049
0.00000000010434267117 0.00000000000778226344 -0.00000000000369680562
0.00000000000051003703 -0.00000000000002058326 -0.00000000000000534812
0.00000000000000122678 -0.00000000000000011813 0.00000000000000000119
0.00000000000000000141 -0.00000000000000000023 0.00000000000000000002 }
→ y a
≪ a DUP SIZE GET
a SIZE 1 - 1 '''FOR''' n
y * a n GET +
-1 '''STEP'''
INV
≫ ≫ '<span style="color:blue>GAMMA</span>' STO
.3 <span style="color:blue>GAMMA</span>
The built-in FACT instruction is obviously based on a similar Taylor formula, since it returns same results:
.3 1 - FACT
{{out}}
<pre>
2: 2.99156898769
1: 2.99156898769
</pre>
=={{header|Ruby}}==
====Taylor series====
{{trans|Ada}}
<
-0.04200_26350_34095_23553, 0.16653_86113_82291_48950, -0.04219_77345_55544_33675,
-0.00962_19715_27876_97356, 0.00721_89432_46663_09954, -0.00116_51675_91859_06511,
Line 3,447 ⟶ 4,408:
end
(1..10).each {|i| puts format("%.14e", gamma(i/3.0))}</
{{out}}
<pre>2.67893853470775e+00
Line 3,460 ⟶ 4,421:
2.77815847933857e+00</pre>
====Built in====
<
{{out}}
<pre>2.678938534707748
Line 3,472 ⟶ 4,433:
2.0
2.7781584804376647
</pre>
=={{header|Rust}}==
====Stirling====
{{trans|AWK}}
<syntaxhighlight lang="rust">
use std::f64::consts;
fn main() {
let gamma = |x: f64| { assert_ne!(x, 0.0); (2.0*consts::PI/x).sqrt() * (x * (x/consts::E).ln()).exp()};
(1..=20).for_each(|x| {
let x = f64::from(x) / 10.0;
println!("{:.02} => {:.10}", x, gamma(x));
});
}
</syntaxhighlight>
{{out}}
<pre>
0.10 => 5.6971871490
0.20 => 3.3259984240
0.30 => 2.3625300363
0.40 => 1.8414763359
0.50 => 1.5203469011
0.60 => 1.3071588574
0.70 => 1.1590532921
0.80 => 1.0533709684
0.90 => 0.9770615079
1.00 => 0.9221370089
1.10 => 0.8834899532
1.20 => 0.8577553354
1.30 => 0.8426782594
1.40 => 0.8367445486
1.50 => 0.8389565525
1.60 => 0.8486932422
1.70 => 0.8656214718
1.80 => 0.8896396353
1.90 => 0.9208427219
2.00 => 0.9595021757
</pre>
=={{header|Scala}}==
<
object Gamma {
Line 3,500 ⟶ 4,501:
println("%.1f -> %.16f %.16f".formatLocal(ENGLISH, x, stGamma(x), laGamma(x)))
}
}</
{{out}}
<pre>Gamma Stirling Lanczos
Line 3,529 ⟶ 4,530:
{{trans|Ruby}} for Taylor
<
(import (scheme base)
(scheme inexact)
Line 3,592 ⟶ 4,593:
(number->string (gamma-taylor i))
"\n")))
</syntaxhighlight>
{{out}}
Line 3,679 ⟶ 4,680:
=={{header|Scilab}}==
<syntaxhighlight lang="text">function x=gammal(z) // Lanczos'
lz=[ 1.000000000190015 ..
76.18009172947146 -86.50532032941677 24.01409824083091 ..
Line 3,700 ⟶ 4,701:
x=i/10
printf("%4.1f %9f %9f\n",x,gamma(x),gammal(x))
end</
{{out}}
<pre style="height:20ex">
Line 3,738 ⟶ 4,739:
=={{header|Seed7}}==
{{trans|Ada}}
<
include "float.s7i";
Line 3,779 ⟶ 4,780:
writeln((gamma(flt(I) / 3.0)) digits 15);
end for;
end func;</
{{out}}
<pre>2.678937911987305
Line 3,794 ⟶ 4,795:
=={{header|Sidef}}==
{{trans|Ruby}}
<
-0.04200_26350_34095_23553, 0.16653_86113_82291_48950, -0.04219_77345_55544_33675,
-0.00962_19715_27876_97356, 0.00721_89432_46663_09954, -0.00116_51675_91859_06511,
Line 3,812 ⟶ 4,813:
for i in 1..10 {
say ("%.14e" % gamma(i/3))
}</
{{out}}
<pre>
Line 3,828 ⟶ 4,829:
Lanczos approximation:
<
var epsilon = 0.0000001
func withinepsilon(x) {
Line 3,864 ⟶ 4,865:
for i in 1..10 {
say ("%.14e" % gamma(i/3))
}</
{{out}}
Line 3,881 ⟶ 4,882:
A simpler implementation:
<
define τ = Num.tau
Line 3,891 ⟶ 4,892:
for i in (1..10) {
say ("%.14e" % Γ(i/3))
}</
{{out}}
Line 3,911 ⟶ 4,912:
The results are compared to Mata's '''[https://www.stata.com/help.cgi?mf_gamma gamma]''' function for each real between 1/100 and 100, by steps of 1/100.
<
_gamma_coef = 1.0,
5.772156649015328606065121e-1,
Line 3,971 ⟶ 4,972:
v=map(&gamma_(),x)
max(abs((v-u):/u))
end</
'''Output'''
Line 3,979 ⟶ 4,980:
Here follows the Python program to compute coefficients.
<
mp.dps = 50
Line 4,009 ⟶ 5,010:
for x in gc:
print(mp.nstr(x, 25))</
=={{header|Tcl}}==
Line 4,015 ⟶ 5,016:
{{tcllib|math}}
{{tcllib|math::calculus}}
<
package require math::calculus
Line 4,049 ⟶ 5,050:
set f2 [expr $f2 * $x]
}
}</
{{out}}
<pre> 1.0 1.000000 1.000000 1.000000
Line 4,076 ⟶ 5,077:
As far as Gamma(n)=(n-1)! , we have everything needed.
===Stirling's approximation===
<
I/2→X
X^(X-1/2)e^(-X)√(2π)→Y
Disp X,(X-1)!,Y
Pause
End</
{{out}}
The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) .
Line 4,119 ⟶ 5,120:
===Lanczos' approximation===
<
I/2→X
e^(ln((1.0
Line 4,132 ⟶ 5,133:
Disp X,(X-1)!,Y
Pause
End</
{{out}}
The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) .
Line 4,168 ⟶ 5,169:
24
</pre>
=={{header|TXR}}==
===Taylor Series===
{{trans|Ada}}
Separator commas in numeric tokens are supported only as of TXR 283.
<syntaxhighlight lang="txrlisp">(defun gamma (x)
(/ (rpoly (- x 1.0)
#( 1.00000,00000,00000,00000 0.57721,56649,01532,86061
-0.65587,80715,20253,88108 -0.04200,26350,34095,23553
0.16653,86113,82291,48950 -0.04219,77345,55544,33675
-0.00962,19715,27876,97356 0.00721,89432,46663,09954
-0.00116,51675,91859,06511 -0.00021,52416,74114,95097
0.00012,80502,82388,11619 -0.00002,01348,54780,78824
-0.00000,12504,93482,14267 0.00000,11330,27231,98170
-0.00000,02056,33841,69776 0.00000,00061,16095,10448
0.00000,00050,02007,64447 -0.00000,00011,81274,57049
0.00000,00001,04342,67117 0.00000,00000,07782,26344
-0.00000,00000,03696,80562 0.00000,00000,00510,03703
-0.00000,00000,00020,58326 -0.00000,00000,00005,34812
0.00000,00000,00001,22678 -0.00000,00000,00000,11813
0.00000,00000,00000,00119 0.00000,00000,00000,00141
-0.00000,00000,00000,00023 0.00000,00000,00000,00002))))
(each ((i 1..11))
(put-line (pic "##.######" (gamma (/ i 3.0)))))</syntaxhighlight>
{{out}}
<pre> 2.678939
1.354118
1.000000
0.892980
0.902745
1.000000
1.190639
1.504575
2.000000
2.778158</pre>
===Stirling===
<syntaxhighlight lang="txrlisp">(defun gamma (x)
(* (sqrt (/ (* 2 %pi%)
x))
(expt (/ x %e%) x)))
(each ((i 1..11))
(put-line (pic "##.######" (gamma (/ i 3.0)))))</syntaxhighlight>
{{out}}
<pre> 2.156976
1.202851
0.922137
0.839743
0.859190
0.959502
1.149106
1.458490
1.945403
2.709764</pre>
===Lanczos===
{{trans|Haskell}}
The Haskell version calculates the natural log of the gamma function, which is why the function is called <code>gammaln</code>; we correct that here by calling <code>exp</code>:
<syntaxhighlight lang="txrlisp">(defun gamma (x)
(let* ((cof #(76.18009172947146 -86.50532032941677
24.01409824083091 -1.231739572450155
0.001208650973866179 -0.000005395239384953))
(ser0 1.000000000190015)
(x55 (+ x 5.5))
(tmp (- x55 (* (+ x 0.5) (log x55))))
(ser (+ ser0 (sum [mapcar / cof (succ x)]))))
(exp (- (log (/ (* 2.5066282746310005 ser) x)) tmp))))
(each ((i (rlist 0.1..1.0..0.1 2..10)))
(put-line (pic "##.# ######.######" i (gamma i))))</syntaxhighlight>
{{out}}
<pre> 0.1 9.513508
0.2 4.590844
0.3 2.991569
0.4 2.218160
0.5 1.772454
0.6 1.489192
0.7 1.298055
0.8 1.164230
0.9 1.068629
1.0 1.000000
2.0 1.000000
3.0 2.000000
4.0 6.000000
5.0 24.000000
6.0 120.000000
7.0 720.000000
8.0 5040.000000
9.0 40320.000000
10.0 362880.000000
</pre>
From Wikipedia Python code. Output is identical to above.
<syntaxhighlight lang="txrlisp">(defun gamma (x)
(if (< x 0.5)
(/ %pi%
(* (sin (* %pi% x))
(gamma (- 1 x))))
(let* ((cof #(676.5203681218851 -1259.1392167224028
771.32342877765313 -176.61502916214059
12.507343278686905 -0.13857109526572012
9.9843695780195716e-6 1.5056327351493116e-7))
(ser0 0.99999999999980993)
(z (pred x))
(tmp (+ z (len cof) -0.5))
(ser (+ ser0 (sum [mapcar / cof (succ z)]))))
(* (sqrt (* 2 %pi%))
(expt tmp (+ z 0.5))
(exp (- tmp))
ser))))
(each ((i (rlist 0.1..1.0..0.1 2..10)))
(put-line (pic "##.# ######.######" i (gamma i))))</syntaxhighlight>
=={{header|Visual FoxPro}}==
Translation of BBC Basic but with OOP extensions. Also some ideas from Numerical Methods (Press ''et al'').
<
LOCAL i As Integer, x As Double, o As lanczos
CLOSE DATABASES ALL
Line 4,232 ⟶ 5,362:
ENDDEFINE
</syntaxhighlight>
{{out}}
<pre>
Line 4,257 ⟶ 5,387:
2.0 1.000000000000000
</pre>
=={{header|V (Vlang)}}==
{{trans|go}}
<syntaxhighlight lang="v (vlang)">import math
fn main() {
println(" x math.Gamma Lanczos7")
for x in [-.5, .1, .5, 1, 1.5, 2, 3, 10, 140, 170] {
println("${x:5.1f} ${math.gamma(x):24.16} ${lanczos7(x):24.16}")
}
}
fn lanczos7(z f64) f64 {
t := z + 6.5
x := .99999999999980993 +
676.5203681218851/z -
1259.1392167224028/(z+1) +
771.32342877765313/(z+2) -
176.61502916214059/(z+3) +
12.507343278686905/(z+4) -
.13857109526572012/(z+5) +
9.9843695780195716e-6/(z+6) +
1.5056327351493116e-7/(z+7)
return math.sqrt2 * math.sqrt_pi * math.pow(t, z-.5) * math.exp(-t) * x
}</syntaxhighlight>
{{out}}
<pre> x math.Gamma Lanczos7
-0.5 -3.544907701811032 -3.544907701811087
0.1 9.513507698668732 9.513507698668752
0.5 1.772453850905516 1.772453850905517
1.0 1 1
1.5 0.8862269254527579 0.8862269254527587
2.0 1 1
3.0 2 2
10.0 362880 362880.0000000015
140.0 9.61572319694107e+238 9.615723196940201e+238
170.0 4.269068009004746e+304 +Inf</pre>
=={{header|Wren}}==
Line 4,263 ⟶ 5,429:
{{libheader|Wren-math}}
The ''gamma'' method in the Math class is based on the Lanczos approximation.
<
import "./math" for Math
var stirling = Fn.new { |x| (2 * Num.pi / x).sqrt * (x / Math.e).pow(x) }
Line 4,271 ⟶ 5,437:
for (i in 1..20) {
var d = i / 10
Fmt.print("$4.2f\t$16.14f\t$16.14f", d, stirling.call(d), Math.gamma(d))
}</syntaxhighlight>
{{out}}
Line 4,282 ⟶ 5,446:
0.10 5.69718714897717 9.51350769866875
0.20 3.32599842402239 4.59084371199881
0.30 2.36253003626962 2.
0.40 1.84147633593624 2.21815954375769
0.50 1.52034690106628 1.77245385090552
Line 4,300 ⟶ 5,464:
1.90 0.92084272189423 0.96176583190739
2.00 0.95950217574449 1.00000000000000
</pre>
=={{header|XPL0}}==
{{trans|Ada}}
<syntaxhighlight lang "XPL0">function real Gamma (X);
real X, A, Y, Sum;
integer N;
begin
A \constant array (0..29) of Long_Float\ :=
[ 1.00000_00000_00000_00000,
0.57721_56649_01532_86061,
-0.65587_80715_20253_88108,
-0.04200_26350_34095_23553,
0.16653_86113_82291_48950,
-0.04219_77345_55544_33675,
-0.00962_19715_27876_97356,
0.00721_89432_46663_09954,
-0.00116_51675_91859_06511,
-0.00021_52416_74114_95097,
0.00012_80502_82388_11619,
-0.00002_01348_54780_78824,
-0.00000_12504_93482_14267,
0.00000_11330_27231_98170,
-0.00000_02056_33841_69776,
0.00000_00061_16095_10448,
0.00000_00050_02007_64447,
-0.00000_00011_81274_57049,
0.00000_00001_04342_67117,
0.00000_00000_07782_26344,
-0.00000_00000_03696_80562,
0.00000_00000_00510_03703,
-0.00000_00000_00020_58326,
-0.00000_00000_00005_34812,
0.00000_00000_00001_22678,
-0.00000_00000_00000_11813,
0.00000_00000_00000_00119,
0.00000_00000_00000_00141,
-0.00000_00000_00000_00023,
0.00000_00000_00000_00002
];
Y := X - 1.0;
Sum := A (29);
for N:= 29-1 downto 0 do
Sum := Sum * Y + A (N);
return 1.0 / Sum;
end \Gamma\;
\Test program:
integer I;
begin
Format(0, 14);
for I:= 1 to 10 do
[RlOut(0, Gamma (Float (I) / 3.0)); CrLf(0)];
end</syntaxhighlight>
{{out}}
<pre>
2.67893853470775E+000
1.35411793942640E+000
1.00000000000000E+000
8.92979511569249E-001
9.02745292950934E-001
1.00000000000000E+000
1.19063934875900E+000
1.50457548825154E+000
1.99999999999397E+000
2.77815847933857E+000
</pre>
=={{header|Yabasic}}==
{{trans|Phix}}
<
sub gamma(z)
Line 4,333 ⟶ 5,563:
for i = 0.1 to 2.1 step .1
print i, " = "; : si(gamma(i))
next</
=={{header|zkl}}==
{{trans|D}} but without a built in gamma function.
<
var table = T(
0x1p+0, 0x1.2788cfc6fb618f4cp-1,
Line 4,359 ⟶ 5,589:
return(1.0 / sm);
}</
<
// Coefficients used by the GNU Scientific Library.
// http://en.wikipedia.org/wiki/Lanczos_approximation
Line 4,387 ⟶ 5,617:
return((2.0 * PI).sqrt() * t.pow(z + 0.5) * exp(-t) * x);
}
}</
{{out}}
<pre>
|