Gamma function: Difference between revisions
Line 1,345: | Line 1,345: | ||
:<math>Gamma(a,x) = \frac{1}{\Gamma(a)} \int_x^{\infty} dt \, t^{a-1} \, exp(-t) </math>, the regularized Gamma function which is also known as the normalized incomplete Gamma function,; |
:<math>Gamma(a,x) = \frac{1}{\Gamma(a)} \int_x^{\infty} dt \, t^{a-1} \, exp(-t) </math>, the regularized Gamma function which is also known as the normalized incomplete Gamma function,; |
||
:<math>GammaRegularizedC(a,x) = \frac{1}{\Gamma(a)} \int_0^x dt \, t^{a-1} \, exp(-t)</math>, which the GSL calls the complementary normalized Gamma function; |
:<math>GammaRegularizedC(a,x) = \frac{1}{\Gamma(a)} \int_0^x dt \, t^{a-1} \, exp(-t)</math>, which the GSL calls the complementary normalized Gamma function; |
||
:<math>LogGamma(a) = \ln \Gamma(a)</math>; |
|||
:<math>RecGamma(a) = \frac{1}{\Gamma(a)}; |
|||
:<math>Pochhammer(a,x) = \frac{\Gamma(a+x)}{\Gamma(x)}. |
|||
=={{header|Ruby}}== |
|||
{{trans|Ada}} |
|||
<lang 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, |
|||
-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 ] |
|||
def gamma(x) |
|||
y = Float(x) - 1 |
|||
1.0 / $a.reverse.inject {|sum, an| sum = sum * y + an} |
|||
end |
|||
(1..10).each {|i| puts format("%.14e", gamma(i/3.0))}</lang> |
|||
Output: |
|||
<pre>2.67893853470775e+00 |
|||
1.35411793942640e+00 |
|||
1.00000000000000e+00 |
|||
8.92979511569249e-01 |
|||
9.02745292950934e-01 |
|||
1.00000000000000e+00 |
|||
1.19063934875900e+00 |
|||
1.50457548825154e+00 |
|||
1.99999999999397e+00 |
|||
2.77815847933857e+00</pre> |
|||
=={{header|Seed7}}== |
|||
{{trans|Ada}} |
|||
<lang seed7>$ include "seed7_05.s7i"; |
|||
include "float.s7i"; |
|||
const func float: gamma (in float: X) is func |
|||
result |
|||
var float: result is 0.0; |
|||
local |
|||
const array float: A is [] ( |
|||
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); |
|||
var float: Y is 0.0; |
|||
var float: Sum is A[maxIdx(A)]; |
|||
var integer: N is 0; |
|||
begin |
|||
Y := X - 1.0; |
|||
for N range pred(maxIdx(A)) downto minIdx(A) do |
|||
Sum := Sum * Y + A[N]; |
|||
end for; |
|||
result := 1.0 / Sum; |
|||
end func; |
|||
const proc: main is func |
|||
local |
|||
var integer: I is 0; |
|||
begin |
|||
for I range 1 to 10 do |
|||
writeln((gamma(flt(I) / 3.0)) digits 15); |
|||
end for; |
|||
end func;</lang> |
|||
Sample output: |
|||
<pre>2.678937911987305 |
|||
1.354117870330811 |
|||
1.000000000000000 |
|||
0.892979443073273 |
|||
0.902745306491852 |
|||
1.000000000000000 |
|||
1.190639257431030 |
|||
1.504575252532959 |
|||
1.999999523162842 |
|||
2.778157949447632</pre> |
|||
=={{header|Tcl}}== |
|||
{{works with|Tcl|8.5}}<br>{{libheader|tcllib}} |
|||
<lang tcl>package require math |
|||
package require math::calculus |
|||
# gamma(1) and gamma(1.5) |
|||
set f 1.0 |
|||
set f2 [expr {sqrt(acos(-1.))/2.}] |
|||
for {set x 1.0} {$x <= 10.0} {set x [expr {$x + 0.5}]} { |
|||
# method 1 - numerical integration, Romberg's method, special |
|||
# case for an improper integral |
|||
set g1 [math::calculus::romberg \ |
|||
[list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \ |
|||
0 1 -relerror 1e-8] |
|||
set g2 [math::calculus::romberg_infinity \ |
|||
[list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \ |
|||
1 Inf -relerror 1e-8] |
|||
set gamma [expr {[lindex $g1 0] + [lindex $g2 0]}] |
|||
# method 2 - library function |
|||
set libgamma [expr {exp([math::ln_Gamma $x])}] |
|||
# method 3 - special forms for integer and half-integer arguments |
|||
if {$x == entier($x)} { |
|||
puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f] |
|||
set f [expr $f * $x] |
|||
} else { |
|||
puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f2] |
|||
set f2 [expr $f2 * $x] |
|||
} |
|||
}</lang> |
|||
Output: |
|||
<pre> 1.0 1.000000 1.000000 1.000000 |
|||
1.5 0.886228 0.886227 0.886227 |
|||
2.0 1.000000 1.000000 1.000000 |
|||
2.5 1.329340 1.329340 1.329340 |
|||
3.0 2.000000 2.000000 2.000000 |
|||
3.5 3.323351 3.323351 3.323351 |
|||
4.0 6.000000 6.000000 6.000000 |
|||
4.5 11.631731 11.631728 11.631728 |
|||
5.0 24.000009 24.000000 24.000000 |
|||
5.5 52.342778 52.342778 52.342778 |
|||
6.0 120.000000 120.000000 120.000000 |
|||
6.5 287.885278 287.885278 287.885278 |
|||
7.0 720.000001 720.000000 720.000000 |
|||
7.5 1871.254311 1871.254305 1871.254306 |
|||
8.0 5040.000032 5039.999999 5040.000000 |
|||
8.5 14034.298267 14034.407291 14034.407293 |
|||
9.0 40320.000705 40319.999992 40320.000000 |
|||
9.5 119292.464880 119292.461971 119292.461995 |
|||
10.0 362880.010950 362879.999927 362880.000000</pre> |
Revision as of 18:37, 20 September 2010
You are encouraged to solve this task according to the task description, using any language you may know.
Implement one algorithm (or more) to compute the Gamma (Γ) function (in the real field only). If your language has the function as builtin or you know a library which has it, compare your implementation's results with the results of the builtin/library function. The Gamma function can be defined as
This suggests a straightforward (but inefficient) way of computing the Γ through numerical integration. Better suggested methods:
Ada
The implementation uses Taylor series coefficients of . The coefficients are taken from Mathematical functions and their approximations by Yudell L. Luke. <lang ada>function Gamma (X : Long_Float) return Long_Float is
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 : constant Long_Float := X - 1.0; Sum : Long_Float := A (A'Last);
begin
for N in reverse A'First..A'Last - 1 loop Sum := Sum * Y + A (N); end loop; return 1.0 / Sum;
end Gamma;</lang> Test program: <lang ada>with Ada.Text_IO; use Ada.Text_IO; with Gamma;
procedure Test_Gamma is begin
for I in 1..10 loop Put_Line (Long_Float'Image (Gamma (Long_Float (I) / 3.0))); end loop;
end Test_Gamma;</lang> Sample output:
2.67893853470775E+00 1.35411793942640E+00 1.00000000000000E+00 8.92979511569249E-01 9.02745292950934E-01 1.00000000000000E+00 1.19063934875900E+00 1.50457548825154E+00 1.99999999999397E+00 2.77815847933858E+00
ALGOL 68
- Stirling & Spouge methods.
- Lanczos method.
<lang algol68># Coefficients used by the GNU Scientific Library # []LONG REAL p = ( LONG 0.99999 99999 99809 93,
LONG 676.52036 81218 851, -LONG 1259.13921 67224 028, LONG 771.32342 87776 5313, -LONG 176.61502 91621 4059, LONG 12.50734 32786 86905, -LONG 0.13857 10952 65720 12, LONG 9.98436 95780 19571 6e-6, LONG 1.50563 27351 49311 6e-7);
PROC lanczos gamma = (LONG REAL in z)LONG REAL: (
# Reflection formula # LONG REAL z := in z; IF z < LONG 0.5 THEN long pi / (long sin(long pi*z)*lanczos gamma(1-z)) ELSE z -:= 1; LONG REAL x := p[1]; FOR i TO UPB p - 1 DO x +:= p[i+1]/(z+i) OD; LONG REAL t = z + UPB p - LONG 1.5; long sqrt(2*long pi) * t**(z+LONG 0.5) * long exp(-t) * x FI
);
PROC taylor gamma = (LONG REAL x)LONG REAL: BEGIN # good for values between 0 and 1 #
[]LONG REAL a = ( LONG 1.00000 00000 00000 00000, LONG 0.57721 56649 01532 86061, -LONG 0.65587 80715 20253 88108, -LONG 0.04200 26350 34095 23553, LONG 0.16653 86113 82291 48950, -LONG 0.04219 77345 55544 33675, -LONG 0.00962 19715 27876 97356, LONG 0.00721 89432 46663 09954, -LONG 0.00116 51675 91859 06511, -LONG 0.00021 52416 74114 95097, LONG 0.00012 80502 82388 11619, -LONG 0.00002 01348 54780 78824, -LONG 0.00000 12504 93482 14267, LONG 0.00000 11330 27231 98170, -LONG 0.00000 02056 33841 69776, LONG 0.00000 00061 16095 10448, LONG 0.00000 00050 02007 64447, -LONG 0.00000 00011 81274 57049, LONG 0.00000 00001 04342 67117, LONG 0.00000 00000 07782 26344, -LONG 0.00000 00000 03696 80562, LONG 0.00000 00000 00510 03703, -LONG 0.00000 00000 00020 58326, -LONG 0.00000 00000 00005 34812, LONG 0.00000 00000 00001 22678, -LONG 0.00000 00000 00000 11813, LONG 0.00000 00000 00000 00119, LONG 0.00000 00000 00000 00141, -LONG 0.00000 00000 00000 00023, LONG 0.00000 00000 00000 00002 ); LONG REAL y = x - 1; LONG REAL sum := a [UPB a]; FOR n FROM UPB a - 1 DOWNTO LWB a DO sum := sum * y + a [n] OD; 1/sum
END # taylor gamma #;
LONG REAL long e = long exp(1);
PROC sterling gamma = (LONG REAL n)LONG REAL: ( # improves for values much greater then 1 #
long sqrt(2*long pi/n)*(n/long e)**n
);
PROC factorial = (LONG INT n)LONG REAL: (
IF n=0 OR n=1 THEN 1 ELIF n=2 THEN 2 ELSE n*factorial(n-1) FI
);
REF[]LONG REAL fm := NIL;
PROC sponge gamma = (LONG REAL x)LONG REAL: (
INT a = 12; # alter to get required precision # REF []LONG REAL fm := NIL; LONG REAL res; IF fm :=: REF[]LONG REAL(NIL) THEN fm := HEAP[0:a-1]LONG REAL; fm[0] := long sqrt(LONG 2*long pi); FOR k TO a-1 DO fm[k] := (((k-1) MOD 2=0) | 1 | -1) * long exp(a-k) *
(a-k) **(k-LONG 0.5) / factorial(k-1)
OD FI; res := fm[0]; FOR k TO a-1 DO res +:= fm[k] / ( x + k ) OD; res *:= long exp(-(x+a)) * (x+a)**(x + LONG 0.5); res/x
);
FORMAT real fmt = $g(-real width, real width - 2)$; FORMAT long real fmt16 = $g(-17, 17 - 2)$; # accurate to about 16 decimal places #
[]STRING methods = ("Genie", "Lanczos", "Sponge", "Taylor","Stirling");
printf(($11xg12xg12xg13xg13xgl$,methods));
FORMAT sample fmt = $"gamma("g(-3,1)")="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$; FORMAT sqr sample fmt = $"gamma("g(-3,1)")**2="f(real fmt)n(UPB methods-1)(", "f(long real fmt16))l$; FORMAT sample exp fmt = $"gamma("g(-3)")="g(-15,11,0)n(UPB methods-1)(","g(-18,14,0))l$;
PROC sample = (LONG REAL x)[]LONG REAL:
(gamma(SHORTEN x), lanczos gamma(x), sponge gamma(x), taylor gamma(x), sterling gamma(x));
FOR i FROM 1 TO 20 DO
LONG REAL x = i / LONG 10; printf((sample fmt, x, sample(x))); IF i = 5 THEN # insert special case of a half # printf((sqr sample fmt, x, gamma(SHORTEN x)**2, lanczos gamma(x)**2, sponge gamma(x)**2, taylor gamma(x)**2, sterling gamma(x)**2)) FI
OD; FOR x FROM 10 BY 10 TO 70 DO
printf((sample exp fmt, x, sample(x)))
OD</lang> Output:
Genie Lanczos Sponge Taylor Stirling gamma(0.1)=9.5135076986687, 9.513507698668730, 9.513507698668731, 9.513509522249043, 5.697187148977169 gamma(0.2)=4.5908437119988, 4.590843711998802, 4.590843711998803, 4.590843743037192, 3.325998424022393 gamma(0.3)=2.9915689876876, 2.991568987687590, 2.991568987687590, 2.991568988322729, 2.362530036269620 gamma(0.4)=2.2181595437577, 2.218159543757688, 2.218159543757688, 2.218159543764845, 1.841476335936235 gamma(0.5)=1.7724538509055, 1.772453850905517, 1.772453850905516, 1.772453850905353, 1.520346901066281 gamma(0.5)**2=3.1415926535898, 3.141592653589795, 3.141592653589793, 3.141592653589216, 2.311454699581843 gamma(0.6)=1.4891922488128, 1.489192248812817, 1.489192248812817, 1.489192248812758, 1.307158857448356 gamma(0.7)=1.2980553326476, 1.298055332647558, 1.298055332647558, 1.298055332647558, 1.159053292113920 gamma(0.8)=1.1642297137253, 1.164229713725304, 1.164229713725303, 1.164229713725303, 1.053370968425609 gamma(0.9)=1.0686287021193, 1.068628702119320, 1.068628702119319, 1.068628702119319, 0.977061507877695 gamma(1.0)=1.0000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 0.922137008895789 gamma(1.1)=0.9513507698669, 0.951350769866873, 0.951350769866873, 0.951350769866873, 0.883489953168704 gamma(1.2)=0.9181687423998, 0.918168742399761, 0.918168742399760, 0.918168742399761, 0.857755335396591 gamma(1.3)=0.8974706963063, 0.897470696306277, 0.897470696306277, 0.897470696306277, 0.842678259448392 gamma(1.4)=0.8872638175031, 0.887263817503076, 0.887263817503075, 0.887263817503064, 0.836744548637082 gamma(1.5)=0.8862269254528, 0.886226925452758, 0.886226925452758, 0.886226925452919, 0.838956552526496 gamma(1.6)=0.8935153492877, 0.893515349287691, 0.893515349287690, 0.893515349288799, 0.848693242152574 gamma(1.7)=0.9086387328533, 0.908638732853291, 0.908638732853290, 0.908638732822421, 0.865621471793884 gamma(1.8)=0.9313837709802, 0.931383770980243, 0.931383770980242, 0.931383769950169, 0.889639635287994 gamma(1.9)=0.9617658319074, 0.961765831907388, 0.961765831907387, 0.961765815012982, 0.920842721894229 gamma(2.0)=1.0000000000000, 1.000000000000000, 0.999999999999999, 1.000000010045742, 0.959502175744492 gamma( 10)= 3.6288000000e5, 3.6288000000000e5, 3.6288000000000e5, 4.051218760300e-7, 3.5986956187410e5 gamma( 20)= 1.216451004e17, 1.216451004088e17, 1.216451004088e17, 1.07701514977e-18, 1.211393423381e17 gamma( 30)= 8.841761994e30, 8.841761993740e30, 8.841761993739e30, 7.98891286318e-23, 8.817236530765e30 gamma( 40)= 2.039788208e46, 2.039788208120e46, 2.039788208120e46, 6.97946184592e-25, 2.035543161237e46 gamma( 50)= 6.082818640e62, 6.082818640343e62, 6.082818640342e62, 1.81016585713e-26, 6.072689187876e62 gamma( 60)= 1.386831185e80, 1.386831185457e80, 1.386831185457e80, 9.27306839649e-28, 1.384906385829e80 gamma( 70)= 1.711224524e98, 1.711224524281e98, 1.711224524281e98, 7.57303907062e-29, 1.709188578191e98
AutoHotkey
Search autohotkey.com: function
Source: AutoHotkey forum by Laszlo
<lang autohotkey>
/*
Here is the upper incomplete Gamma function. Omitting or setting
the second parameter to 0 we get the (complete) Gamma function.
The code is based on: "Computation of Special Functions" Zhang and Jin,
John Wiley and Sons, 1996
- /
SetFormat FloatFast, 0.9e
Loop 10
MsgBox % GAMMA(A_Index/3) "`n" GAMMA(A_Index*10)
GAMMA(a,x=0) { ; upper incomplete gamma: Integral(t**(a-1)*e**-t, t = x..inf)
If (a > 171 || x < 0) Return 2.e308 ; overflow
xam := x > 0 ? -x+a*ln(x) : 0 If (xam > 700) Return 2.e308 ; overflow
If (x > 1+a) { ; no need for gamma(a) t0 := 0, k := 60 Loop 60 t0 := (k-a)/(1+k/(x+t0)), --k Return exp(xam) / (x+t0) }
r := 1, ga := 1.0 ; compute ga = gamma(a) ... If (a = round(a)) ; if integer: factorial If (a > 0) Loop % a-1 ga *= A_Index Else ; negative integer ga := 1.7976931348623157e+308 ; Dmax Else { ; not integer If (abs(a) > 1) { z := abs(a) m := floor(z) Loop %m% r *= (z-A_Index) z -= m } Else z := a
gr := ((((((((((((((((((((((( 0.14e-14 *z - 0.54e-14) *z - 0.206e-13) *z + 0.51e-12) *z - 0.36968e-11) *z + 0.77823e-11) *z + 0.1043427e-9) *z - 0.11812746e-8) *z + 0.50020075e-8) *z + 0.6116095e-8) *z - 0.2056338417e-6) *z + 0.1133027232e-5) *z - 0.12504934821e-5) *z - 0.201348547807e-4) *z + 0.1280502823882e-3) *z - 0.2152416741149e-3) *z - 0.11651675918591e-2) *z + 0.7218943246663e-2) *z - 0.9621971527877e-2) *z - 0.421977345555443e-1) *z + 0.1665386113822915) *z - 0.420026350340952e-1) *z - 0.6558780715202538) *z + 0.5772156649015329) *z + 1
ga := 1.0/(gr*z) * r If (a < -1) ga := -3.1415926535897931/(a*ga*sin(3.1415926535897931*a)) }
If (x = 0) ; complete gamma requested Return ga
s := 1/a ; here x <= 1+a r := s Loop 60 { r *= x/(a+A_Index) s += r If (abs(r/s) < 1.e-15) break } Return ga - exp(xam)*s
}
/* The 10 results shown: 2.678938535e+000 1.354117939e+000 1.0 8.929795115e-001 9.027452930e-001 3.628800000e+005 1.216451004e+017 8.841761994e+030 2.039788208e+046 6.082818640e+062
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
- /
</lang>
C
This implements Stirling's approximation and Spouge's approximation.
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- include <gsl/gsl_sf_gamma.h>
- ifndef M_PI
- define M_PI 3.14159265358979323846
- endif
/* very simple approximation */ double st_gamma(double x) {
return sqrt(2.0*M_PI/x)*pow(x/M_E, x);
}
- define A 12
double sp_gamma(double z) {
const int a = A; static double c_space[A]; static double *c = NULL; int k; double accm; if ( c == NULL ) { double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/ c = c_space; c[0] = sqrt(2.0*M_PI); for(k=1; k < a; k++) { c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
k1_factrl *= -k;
} } accm = c[0]; for(k=1; k < a; k++) { accm += c[k] / ( z + k ); } accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */ return accm/z;
}
int main() {
double x;
printf("%15s%15s%15s\n", "Stirling", "Spouge", "GSL"); for(x=1.0; x <= 10.0; x+=1.0) { printf("%15.8lf%15.8lf%15.8lf\n", st_gamma(x/3.0), sp_gamma(x/3.0),
gsl_sf_gamma(x/3.0));
} return 0;
}</lang>
Output:
Stirling Spouge GSL 2.15697602 2.67893853 2.67893853 1.20285073 1.35411794 1.35411794 0.92213701 1.00000000 1.00000000 0.83974270 0.89297951 0.89297951 0.85919025 0.90274529 0.90274529 0.95950218 1.00000000 1.00000000 1.14910642 1.19063935 1.19063935 1.45849038 1.50457549 1.50457549 1.94540320 2.00000000 2.00000000 2.70976382 2.77815848 2.77815848
Common Lisp
<lang lisp>
- Taylor series coefficients
(defconstant tcoeff
'( 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))
- number of coefficients
(defconstant numcoeff (length tcoeff))
(defun gamma (x)
(let ((y (- x 1.0)) (sum (nth (- numcoeff 1) tcoeff))) (loop for i from (- numcoeff 2) downto 0 do (setf sum (+ (* sum y) (nth i tcoeff)))) (/ 1.0 sum)))
(loop for i from 1 to 10
do ( format t "~12,10f~%" (gamma (/ i 3.0))))
</lang>
Produces:
2.6789380000 1.3541179000 1.0000000000 0.8929794500 0.9027453000 1.0000000000 1.1906393000 1.5045753000 1.9999995000 2.7781580000
D
<lang d>import std.stdio: writefln;
real gamma(const real x) {
static enum real[30] a = [ 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];
immutable real y = x - 1.0; auto sm = a[$ - 1]; foreach_reverse (an; a[0 .. $-1]) sm = sm * y + an; return 1.0 / sm;
}
void main() {
foreach (i; 1 .. 11) writefln("%f: %20.20e", i / 3.0L, gamma(i / 3.0L));
}</lang> Output:
0.333333: 2.67893853470774763300e+00 0.666667: 1.35411793942640041690e+00 1.000000: 1.00000000000000000000e+00 1.333333: 8.92979511569249211240e-01 1.666667: 9.02745292950933611320e-01 2.000000: 1.00000000000000000010e+00 2.333333: 1.19063934875899893320e+00 2.666667: 1.50457548825154009640e+00 3.000000: 1.99999999999396790280e+00 3.333333: 2.77815847933857519470e+00
Fortran
This code shows two methods: Numerical Integration through Simpson formula, and Lanczos approximation. The results of testing are printed altogether with the values given by the function gamma; this function is defined in the Fortran 2008 standard, and supported by GNU Fortran (and other vendors) as extension; if not present in your compiler, you can remove the last part of the print in order to get it compiled with any Fortran 95 compliant compiler.
<lang fortran>program ComputeGammaInt
implicit none
integer :: i
write(*, "(3A15)") "Simpson", "Lanczos", "Builtin" do i=1, 10 write(*, "(3F15.8)") my_gamma(i/3.0), lacz_gamma(i/3.0), gamma(i/3.0) end do
contains
pure function intfuncgamma(x, y) result(z) real :: z real, intent(in) :: x, y z = x**(y-1.0) * exp(-x) end function intfuncgamma
function my_gamma(a) result(g) real :: g real, intent(in) :: a
real, parameter :: small = 1.0e-4 integer, parameter :: points = 100000
real :: infty, dx, p, sp(2, points), x integer :: i logical :: correction
x = a
correction = .false. ! value with x<1 gives \infty, so we use ! \Gamma(x+1) = x\Gamma(x) ! to avoid the problem if ( x < 1.0 ) then correction = .true. x = x + 1 end if
! find a "reasonable" infinity... ! we compute this integral indeed ! \int_0^M dt t^{x-1} e^{-t} ! where M is such that M^{x-1} e^{-M} ≤ \epsilon infty = 1.0e4 do while ( intfuncgamma(infty, x) > small ) infty = infty * 10.0 end do
! using simpson dx = infty/real(points) sp = 0.0 forall(i=1:points/2-1) sp(1, 2*i) = intfuncgamma(2.0*(i)*dx, x) forall(i=1:points/2) sp(2, 2*i - 1) = intfuncgamma((2.0*(i)-1.0)*dx, x) g = (intfuncgamma(0.0, x) + 2.0*sum(sp(1,:)) + 4.0*sum(sp(2,:)) + & intfuncgamma(infty, x))*dx/3.0
if ( correction ) g = g/a
end function my_gamma
recursive function lacz_gamma(a) result(g) real, intent(in) :: a real :: g
real, parameter :: pi = 3.14159265358979324 integer, parameter :: cg = 7
! these precomputed values are taken by the sample code in Wikipedia, ! and the sample itself takes them from the GNU Scientific Library real, dimension(0:8), parameter :: p = & (/ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, & 771.32342877765313, -176.61502916214059, 12.507343278686905, & -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 /)
real :: t, w, x integer :: i
x = a
if ( x < 0.5 ) then g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) ) else x = x - 1.0 t = p(0) do i=1, cg+2 t = t + p(i)/(x+real(i)) end do w = x + real(cg) + 0.5 g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t end if end function lacz_gamma
end program ComputeGammaInt</lang>
Output:
Simpson Lanczos Builtin 2.65968132 2.67893744 2.67893839 1.35269761 1.35411859 1.35411787 1.00000060 1.00000024 1.00000000 0.88656044 0.89297968 0.89297950 0.90179849 0.90274525 0.90274531 0.99999803 1.00000036 1.00000000 1.19070935 1.19063985 1.19063926 1.50460517 1.50457609 1.50457561 2.00000286 2.00000072 2.00000000 2.77815390 2.77816010 2.77815843
Groovy
<lang groovy>a = [ 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].reverse()
def gamma = { 1.0 / a.inject(0) { sm, a_i -> sm * (it - 1) + a_i } }
(1..10).each{ printf("% 1.9e\n", gamma(it / 3.0)) } </lang> Sample output:
2.678938535e+00 1.354117939e+00 1.000000000e+00 8.929795116e-01 9.027452930e-01 1.000000000e+00 1.190639349e+00 1.504575488e+00 2.000000000e+00 2.778158479e+00
Haskell
From HaskellWiki (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]
<lang haskell>cof :: [Double] cof = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953]
ser :: Double ser = 1.000000000190015
gammaln :: Double -> Double gammaln xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5)
ser' = ser + sum (zipWith (/) cof [xx+1..]) in -tmp' + log(2.5066282746310005 * ser' / xx)</lang>
J
This code shows the built-in method, which works for any value (positive, negative and complex numbers).
<lang j>gamma=: !@<:</lang>
The following direct coding of the task comes from the Stirling's approximation essay on the J wiki:
<lang j>sbase =: %:@(2p1&%) * %&1x1 ^ ] scorr =: 1 1r12 1r288 _139r51840 _571r2488320&p.@% stirlg=: sbase * scorr</lang>
Checking against !@<:
we can see that this approximation loses accuracy for small arguments
<lang j> (,. stirlg ,. gamma) 10 1p1 1x1 1.5 1
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</lang>
Java
Implementation of Stirling's approximation and Lanczos approximation. <lang java> import java.lang.Math;
public class GammaFunction {
public double st_gamma(double x){ return Math.sqrt(2*Math.PI/x)*Math.pow((x/Math.E), x); }
public double la_gamma(double x){ double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7}; int g = 7; if(x < 0.5) return Math.PI / (Math.sin(Math.PI * x)*la_gamma(1-x)); else{ x -= 1; double a = p[0]; double t = x+g+0.5; for(int i = 1; i < p.length; i++){ a += p[i]/(x+i); }
return Math.sqrt(2*Math.PI)*Math.pow(t, x+0.5)*Math.exp(-t)*a; } }
public static void main(String[] args) { GammaFunction test = new GammaFunction(); System.out.println("Gamma \t\tStirling \t\tLanczos"); for(double i = 1.0; i <= 20.0 ; i += 1.0){ System.out.println("" + i/10.0 + "\t\t" + test.st_gamma(i/10.0) + "\t" + test.la_gamma(i/10.0)); } } } </lang>
Here are some sample Outputs:
Gamma Stirling Lanczos 0.1 5.697187148977169 9.513507698668734 0.2 3.3259984240223925 4.590843711998803 0.3 2.3625300362696198 2.9915689876875904 0.4 1.8414763359362354 2.218159543757687 0.5 1.5203469010662807 1.7724538509055159 0.6 1.307158857448356 1.489192248812818 0.7 1.15905329211392 1.2980553326475577 0.8 1.0533709684256085 1.1642297137253035 0.9 0.9770615078776954 1.0686287021193193 1.0 0.9221370088957891 0.9999999999999998 1.1 0.8834899531687038 0.9513507698668735 1.2 0.8577553353965909 0.9181687423997607 1.3 0.8426782594483921 0.8974706963062777 1.4 0.8367445486370817 0.8872638175030757 1.5 0.8389565525264963 0.8862269254527586 1.6 0.8486932421525738 0.8935153492876909 1.7 0.865621471793884 0.9086387328532916 1.8 0.8896396352879945 0.9313837709802425 1.9 0.9208427218942293 0.9617658319073877 2.0 0.9595021757444916 1.0000000000000002
Lua
Uses the wp:Reciprocal gamma function to calculate small values. <lang lua> gamma, coeff, quad, qui, set = 0.577215664901, -0.65587807152056, -0.042002635033944, 0.16653861138228, -0.042197734555571 function recigamma(z)
return z + gamma * z^2 + coeff * z^3 + quad * z^4 + qui * z^5 + set * z^6
end
function gammafunc(z)
if z == 1 then return 1 elseif math.abs(z) <= 0.5 then return 1 / recigamma(z) else return (z - 1) * gammafunc(z-1) end
end </lang>
Mathematica
This code shows the built-in method, which works for any value (positive, negative and complex numbers). <lang mathematica>
Gamma[x]
</lang>
Output integers and half-integers (a space is multiplication in Mathematica):
1/2 Sqrt[pi] 1 1 3/2 Sqrt[pi]/2 2 1 5/2 (3 Sqrt[pi])/4 3 2 7/2 (15 Sqrt[pi])/8 4 6 9/2 (105 Sqrt[pi])/16 5 24 11/2 (945 Sqrt[pi])/32 6 120 13/2 (10395 Sqrt[pi])/64 7 720 15/2 (135135 Sqrt[pi])/128 8 5040 17/2 (2027025 Sqrt[pi])/256 9 40320 19/2 (34459425 Sqrt[pi])/512 10 362880
Output approximate numbers:
0.1 9.51351 0.2 4.59084 0.3 2.99157 0.4 2.21816 0.5 1.77245 0.6 1.48919 0.7 1.29806 0.8 1.16423 0.9 1.06863 1. 1.
Output complex numbers:
I -0.15495-0.498016 I 2 I 0.00990244-0.075952 I 3 I 0.0112987-0.00643092 I 4 I 0.00173011+0.00157627 I 5 I -0.000271704+0.000339933 I
Modula-3
<lang modula3>MODULE Gamma EXPORTS Main;
FROM IO IMPORT Put; FROM Fmt IMPORT Extended, Style;
PROCEDURE Taylor(x: EXTENDED): EXTENDED =
CONST a = ARRAY [0..29] OF EXTENDED { 1.00000000000000000000X0, 0.57721566490153286061X0, -0.65587807152025388108X0, -0.04200263503409523553X0, 0.16653861138229148950X0, -0.04219773455554433675X0, -0.00962197152787697356X0, 0.00721894324666309954X0, -0.00116516759185906511X0, -0.00021524167411495097X0, 0.00012805028238811619X0, -0.00002013485478078824X0, -0.00000125049348214267X0, 0.00000113302723198170X0, -0.00000020563384169776X0, 0.00000000611609510448X0, 0.00000000500200764447X0, -0.00000000118127457049X0, 0.00000000010434267117X0, 0.00000000000778226344X0, -0.00000000000369680562X0, 0.00000000000051003703X0, -0.00000000000002058326X0, -0.00000000000000534812X0, 0.00000000000000122678X0, -0.00000000000000011813X0, 0.00000000000000000119X0, 0.00000000000000000141X0, -0.00000000000000000023X0, 0.00000000000000000002X0 }; VAR y := x - 1.0X0; sum := a[LAST(a)];
BEGIN FOR i := LAST(a) - 1 TO FIRST(a) BY -1 DO sum := sum * y + a[i]; END; RETURN 1.0X0 / sum; END Taylor;
BEGIN
FOR i := 1 TO 10 DO Put(Extended(Taylor(FLOAT(i, EXTENDED) / 3.0X0), style := Style.Sci) & "\n"); END;
END Gamma.</lang> Output:
2.6789385347077490e+000 1.3541179394264005e+000 1.0000000000000000e+000 8.9297951156924930e-001 9.0274529295093360e-001 1.0000000000000000e+000 1.1906393487589992e+000 1.5045754882515399e+000 1.9999999999939684e+000 2.7781584793385790e+000
OCaml
<lang ocaml>
let e = exp 1. let pi = 4. *. atan 1. let sqrttwopi = sqrt (2. *. pi)
module Lanczos = struct
(* Lanczos method *) (* Coefficients used by the GNU Scientific Library *) let g = 7. let c = [|0.99999999999980993; 676.5203681218851; -1259.1392167224028;
771.32342877765313; -176.61502916214059; 12.507343278686905; -0.13857109526572012; 9.9843695780195716e-6; 1.5056327351493116e-7|]
let rec ag z d = if d = 0 then c.(0) +. ag z 1 else if d < 8 then c.(d) /. (z +. float d) +. ag z (succ d) else c.(d) /. (z +. float d)
let gamma z = let z = z -. 1. in let p = z +. g +. 0.5 in sqrttwopi *. p ** (z +. 0.5) *. exp (-. p) *. ag z 0
end
module Stirling = struct
(* Stirling method *) let gamma z = sqrttwopi /. sqrt z *. (z /. e) ** z
end
module Stirling2 = struct
(* Extended Stirling method seen in Abramowitz and Stegun *) let d = [|1./.12.; 1./.288.; -139./.51840.; -571./.2488320.|]
let rec corr z x n = if n < Array.length d - 1 then d.(n) /. x +. corr z (x *. z) (succ n) else d.(n) /. x
let gamma z = Stirling.gamma z *. (1. +. corr z z 0)
end
let mirror gma z =
if z > 0.5 then gma z else pi /. sin (pi *. z) /. gma (1. -. z)
let _ =
Printf.printf "z\t\tLanczos\t\tStirling\tStirling2\n"; for i = 1 to 20 do let z = float i /. 10. in Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n" z
(mirror Lanczos.gamma z) (mirror Stirling.gamma z) (mirror Stirling2.gamma z)
done; for i = 1 to 7 do let z = 10. *. float i in Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n" z
(Lanczos.gamma z) (Stirling.gamma z) (Stirling2.gamma z)
done
</lang>
Output
z Lanczos Stirling Stirling2 0.1 9.51350770e+00 1.04050843e+01 9.52104183e+00 0.2 4.59084371e+00 5.07399275e+00 4.59686230e+00 0.3 2.99156899e+00 3.35033954e+00 2.99844028e+00 0.4 2.21815954e+00 2.52705781e+00 2.22775889e+00 0.5 1.77245385e+00 2.06636568e+00 1.78839014e+00 0.6 1.48919225e+00 1.30715886e+00 1.48277536e+00 0.7 1.29805533e+00 1.15905329e+00 1.29508068e+00 0.8 1.16422971e+00 1.05337097e+00 1.16270541e+00 0.9 1.06862870e+00 9.77061508e-01 1.06778308e+00 1 1.00000000e+00 9.22137009e-01 9.99499469e-01 1.1 9.51350770e-01 8.83489953e-01 9.51037997e-01 1.2 9.18168742e-01 8.57755335e-01 9.17964058e-01 1.3 8.97470696e-01 8.42678259e-01 8.97331287e-01 1.4 8.87263818e-01 8.36744549e-01 8.87165485e-01 1.5 8.86226925e-01 8.38956553e-01 8.86155384e-01 1.6 8.93515349e-01 8.48693242e-01 8.93461840e-01 1.7 9.08638733e-01 8.65621472e-01 9.08597702e-01 1.8 9.31383771e-01 8.89639635e-01 9.31351590e-01 1.9 9.61765832e-01 9.20842722e-01 9.61740068e-01 2 1.00000000e+00 9.59502176e-01 9.99978981e-01 10 3.62880000e+05 3.59869562e+05 3.62879997e+05 20 1.21645100e+17 1.21139342e+17 1.21645100e+17 30 8.84176199e+30 8.81723653e+30 8.84176199e+30 40 2.03978821e+46 2.03554316e+46 2.03978821e+46 50 6.08281864e+62 6.07268919e+62 6.08281864e+62 60 1.38683119e+80 1.38490639e+80 1.38683119e+80 70 1.71122452e+98 1.70918858e+98 1.71122452e+98
Octave
<lang octave>function g = lacz_gamma(a, cg=7)
p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \ 771.32342877765313, -176.61502916214059, 12.507343278686905, \ -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ]; x=a; if ( x < 0.5 ) g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) ); else x = x - 1.0; t = p(1); for i=1:(cg+1) t = t + p(i+1)/(x+double(i)); endfor w = x + double(cg) + 0.5; g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t; endif
endfunction
for i = 1:10
printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));
endfor</lang>
Output
2.678939 2.678939 1.354118 1.354118 1.000000 1.000000 0.892980 0.892980 0.902745 0.902745 1.000000 1.000000 1.190639 1.190639 1.504575 1.504575 2.000000 2.000000 2.778158 2.778158
Which suggests that the built-in gamma uses the same approximation.
Perl 6
<lang perl6>my @coeff = reverse <
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
>;
sub gamma($x) {
my $y = $x - 1; 1 / @coeff.reduce: { $^a * $y + $^b };
}
for 1/3, 2/3 ... 10/3 -> $n {
say $n.fmt("%.5f"), " ", gamma($n).fmt("%.14e");
}</lang>
PL/I
<lang PL/I> /* From Rosetta Fortran */ test: procedure options (main);
declare i fixed binary;
on underflow ;
put skip list ('Lanczos', 'Builtin' ); do i = 1 to 10; put skip list (lanczos_gamma(i/3.0q0), gamma(i/3.0q0) ); end;
lanczos_gamma: procedure (a) returns (float (18)) recursive;
declare a float (18); declare pi float (18) value (3.14159265358979324E0); declare cg fixed binary initial ( 7 );
/* these precomputed values are taken by the sample code in Wikipedia, */ /* and the sample itself takes them from the GNU Scientific Library */ declare p(0:8) float (18) static initial ( 0.99999999999980993e0, 676.5203681218851e0, -1259.1392167224028e0, 771.32342877765313e0, -176.61502916214059e0, 12.507343278686905e0, -0.13857109526572012e0, 9.9843695780195716e-6, 1.5056327351493116e-7 );
declare ( t, w, x ) float (18); declare i fixed binary;
x = a;
if x < 0.5 then return ( pi / ( sin(pi*x) * lanczos_gamma(1.0-x) ) ); else do; x = x - 1.0; t = p(0); do i = 1 to cg+2; t = t + p(i)/(x+i); end; w = x + float(cg) + 0.5; return ( sqrt(2*pi) * w**(x+0.5) * exp(-w) * t ); end; end lanczos_gamma;
end test; </lang> Output: <lang> Lanczos Builtin
2.67893853470774706E+0000 2.678938534707747630E+0000 1.35411793942640071E+0000 1.354117939426400420E+0000 1.00000000000000021E+0000 1.000000000000000000E+0000 8.92979511569249470E-0001 8.929795115692492110E-0001 9.02745292950933961E-0001 9.027452929509336110E-0001 1.00000000000000048E+0000 1.000000000000000000E+0000 1.19063934875899964E+0000 1.190639348758998950E+0000 1.50457548825155704E+0000 1.504575488251556020E+0000 2.00000000000000154E+0000 2.000000000000000000E+0000 2.77815848043766660E+0000 2.778158480437664210E+0000
</lang>
PicoLisp
<lang PicoLisp>(scl 28)
(de *A
~(flip (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 ) ) )
(de gamma (X)
(let (Y (- X 1.0) Sum (car *A)) (for A (cdr *A) (setq Sum (+ A (*/ Sum Y 1.0))) ) (*/ 1.0 1.0 Sum) ) )</lang>
Output:
: (for I (range 1 10) (prinl (round (gamma (*/ I 1.0 3)) 14)) ) 2.67893853470775 1.35411793942640 1.00000000000000 0.89297951156925 0.90274529295093 1.00000000000000 1.19063934875900 1.50457548825154 1.99999999999397 2.77815847933858
PureBasic
Below is PureBasic code for: Complete Gamma function Natural Logarithm of the Complete Gamma function Factorial function
<lang PureBasic>Procedure.d Gamma(x.d) ; AKJ 01-May-10
- Complete Gamma function for x>0 and x<2 (approx)
- Extended outside this range via
- Gamma(x+1) = x*Gamma(x)
- Based on http://rosettacode.org/wiki/Gamma_function [Ada]
Protected Dim A.d(28) A(0) = 1.0 A(1) = 0.5772156649015328606 A(2) =-0.6558780715202538811 A(3) =-0.0420026350340952355 A(4) = 0.1665386113822914895 A(5) =-0.0421977345555443368 ; was ...33675 A(6) =-0.0096219715278769736 A(7) = 0.0072189432466630995 A(8) =-0.0011651675918590651 A(9) =-0.0002152416741149510 A(10) = 0.0001280502823881162 A(11) =-0.0000201348547807882 A(12) =-0.0000012504934821427 A(13) = 0.0000011330272319817 A(14) =-0.0000002056338416978 A(15) = 0.0000000061160951045 A(16) = 0.0000000050020076445 A(17) =-0.0000000011812745705 A(18) = 0.0000000001043426712 A(19) = 0.0000000000077822634 A(20) =-0.0000000000036968056 A(21) = 0.0000000000005100370 A(22) =-0.0000000000000205833 A(23) =-0.0000000000000053481 A(24) = 0.0000000000000012268 A(25) =-0.0000000000000001181 A(26) = 0.0000000000000000012 A(27) = 0.0000000000000000014 A(28) =-0.0000000000000000002
- A(29) = 0.00000000000000000002
Protected y.d, Prod.d, Sum.d, N If x<=0.0: ProcedureReturn 0.0: EndIf ; Error y = x-1.0: Prod = 1.0 While y>=1.0
Prod*y: y-1.0 ; Recurse using Gamma(x+1) = x*Gamma(x)
Wend Sum= A(28) For N = 27 To 0 Step -1
Sum*y+A(N)
Next N If Sum=0.0: ProcedureReturn Infinity(): EndIf ProcedureReturn Prod / Sum EndProcedure
Procedure.d GammLn(x.d) ; AKJ 01-May-10
- Returns Ln(Gamma()) or 0 on error
- Based on Numerical Recipes gamma.h
Protected j, tmp.d, y.d, ser.d Protected Dim cof.d(13) cof(0) = 57.1562356658629235 cof(1) = -59.5979603554754912 cof(2) = 14.1360979747417471 cof(3) = -0.491913816097620199 cof(4) = 0.339946499848118887e-4 cof(5) = 0.465236289270485756e-4 cof(6) = -0.983744753048795646e-4 cof(7) = 0.158088703224912494e-3 cof(8) = -0.210264441724104883e-3 cof(9) = 0.217439618115212643e-3 cof(10) = -0.164318106536763890e-3 cof(11) = 0.844182239838527433e-4 cof(12) = -0.261908384015814087e-4 cof(13) = 0.368991826595316234e-5 If x<=0: ProcedureReturn 0: EndIf ; Bad argument y = x tmp = x+5.2421875 tmp = (x+0.5)*Log(tmp)-tmp ser = 0.999999999999997092 For j=0 To 13
y + 1: ser + cof(j)/y
Next j ProcedureReturn tmp+Log(2.5066282746310005*ser/x) EndProcedure
Procedure Factorial(x) ; AKJ 01-May-10
ProcedureReturn Gamma(x+1)
EndProcedure</lang> Examples: <lang PureBasic>Debug "Gamma()" For i = 12 To 15
Debug StrD(i/3.0, 3)+" "+StrD(Gamma(i/3.0))
Next i Debug "" Debug "Ln(Gamma(5.0)) = "+StrD(GammLn(5.0), 16) ; Ln(24) Debug "" Debug "Factorial 6 = "+StrD(Factorial(6), 0) ; 72</lang> Output:
[Debug] Gamma(): [Debug] 4.000 6.0000000000 [Debug] 4.333 9.2605282681 [Debug] 4.667 14.7114047740 [Debug] 5.000 24.0000000000 [Debug] [Debug] Ln(Gamma(5.0)) = 3.1780538303479458 [Debug] [Debug] Factorial 6 = 720
Python
<lang python>_a = ( 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 )
def gamma (x):
y = float(x) - 1.0; sm = _a[-1]; for an in _a[-2::-1]: sm = sm * y + an; return 1.0 / sm;
if __name__ == '__main__':
for i in range(1,11): print " %20.14e" % gamma(i/3.0)
</lang> Sample output:
2.67893853470775e+00 1.35411793942640e+00 1.00000000000000e+00 8.92979511569249e-01 9.02745292950934e-01 1.00000000000000e+00 1.19063934875900e+00 1.50457548825154e+00 1.99999999999397e+00 2.77815847933857e+00
R
Lanczos' approximation is a loose
.
<lang r> stirling <- function(z) sqrt(2*pi/z) * (exp(-1)*z)^z
nemes <- function(z) sqrt(2*pi/z) * (exp(-1)*(z + (12*z - (10*z)^-1)^-1))^z
lanczos <- function(z) {
if(length(z) > 1) { sapply(z, lanczos) } else { g <- 7 p <- c(0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7) z <- as.complex(z) if(Re(z) < 0.5) { pi / (sin(pi*z) * lanczos(1-z)) } else { z <- z - 1 x <- p[1] + sum(p[-1]/seq.int(z+1, z+g+1)) tt <- z + g + 0.5 sqrt(2*pi) * tt^(z+0.5) * exp(-tt) * x } }
}
spouge <- function(z, a=49) {
if(length(z) > 1) { sapply(z, spouge) } else { z <- z-1 k <- seq.int(1, a-1) ck <- rep(c(1, -1), len=a-1) / factorial(k-1) * (a-k)^(k-0.5) * exp(a-k) (z + a)^(z+0.5) * exp(-z - a) * (sqrt(2*pi) + sum(ck/(z+k))) }
}
- Checks
z <- (1:10)/3 all.equal(gamma(z), stirling(z)) # Mean relative difference: 0.07181942 all.equal(gamma(z), nemes(z)) # Mean relative difference: 0.003460549 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)) </lang>
z stirling nemes lanczos spouge builtin 1 0.3333333 2.1569760 2.6290752 2.6789385+0i 2.6789385 2.6789385 2 0.6666667 1.2028507 1.3515736 1.3541179+0i 1.3541179 1.3541179 3 1.0000000 0.9221370 0.9996275 1.0000000+0i 1.0000000 1.0000000 4 1.3333333 0.8397427 0.8928835 0.8929795+0i 0.8929795 0.8929795 5 1.6666667 0.8591902 0.9027098 0.9027453+0i 0.9027453 0.9027453 6 2.0000000 0.9595022 0.9999831 1.0000000+0i 1.0000000 1.0000000 7 2.3333333 1.1491064 1.1906296 1.1906393+0i 1.1906393 1.1906393 8 2.6666667 1.4584904 1.5045690 1.5045755+0i 1.5045755 1.5045755 9 3.0000000 1.9454032 1.9999951 2.0000000+0i 2.0000000 2.0000000 10 3.3333333 2.7097638 2.7781544 2.7781585+0i 2.7781585 2.7781585
RLaB
RLaB through GSL has the following functions related to the Gamma function:
- , the Gamma function;
- , the regularized Gamma function which is also known as the normalized incomplete Gamma function,;
- , which the GSL calls the complementary normalized Gamma function;