Gamma function: Difference between revisions

Content deleted Content added
→‎{{header|R}}: Issue found in previous R implementation step using 'seq.int' function: imaginary parts discarded in coercion; The edit resolved the issue and retain the imaginary part to make sure the precision. The edited lanczos(z) was benchmarked with gammaz(z) and complex_gamma(z) from 'pracma' and 'hypergeo' libraries to cross validate the solution accuracy
Zeddicus (talk | contribs)
 
(11 intermediate revisions by 8 users not shown)
Line 492:
gamma( 70)= 1.711224524e98, 1.711224524281e98, 1.711224524281e98, 7.57303907062e-29, 1.709188578191e98
</pre>
 
=={{header|ANSI Standard BASIC}}==
 
{{trans|BBC Basic}} - Lanczos method.
<syntaxhighlight lang="ansi standard 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>
 
=={{header|Arturo}}==
Line 703 ⟶ 673:
</pre>
 
=={{header|BASIC256BASIC}}==
==={{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.
 
Line 737 ⟶ 761:
end function</syntaxhighlight>
 
==={{header|BBC BASIC}}===
{{works with|BBC BASIC for Windows}}
Uses the Lanczos approximation.
Line 1,260 ⟶ 1,284:
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}}==
Line 2,279 ⟶ 2,321:
<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}}==
<syntaxhighlight lang="scala">// version 1.0.6
Line 4,029 ⟶ 4,153:
 
Note: &nbsp; The Taylor series isn't much good above values of &nbsp; <big>'''6½'''</big>.
Already on modest values of x (say > 3) you loose precision. See below for a solution.
<syntaxhighlight lang="rexx">/*REXX program calculates GAMMA using the Taylor series coefficients; ≈80 decimal digits*/
/*The GAMMA function symbol is the Greek capital letter: Γ */
Line 4,201 ⟶ 4,326:
3.3233509704478425511840640312646472177454 3.1415926535897932384626433832795028841972
6
</pre>
This program also has some minor issues. The number of steps, as well as e and pi, are hardcoded. This limits the precision to about 80 digits. See below for a solution.
 
===Libraries===
REXX does not have an 'include' facility nor 'arbitrary precision' mathematical functions out of the
box. In previous REXX entries it was custom the have all needed routines or procedures copied into
the program. Newer entries redirect to<br>
[[Mathematics.rex|Libraries]] and<br>
[[Compiler/Simple file inclusion pre processor#rexx|Pre processor and Include clause]].<br>
At the end of each entry you'll find several Include clauses, showing the libraries that the program
needs (cf '#include', 'import', 'uses' or '::requires').
===A generic Gamma function===
Below program calculates the Gamma function for any (real) element x (except integers < 1) in the specified precision (but over 100 digits it quickly becomes slow).<br>
Improvements over above programs:
 
Closed formulas added for all half integers and positive integers (much faster than series). Argument reduction (mapping to interval 0.5...1.5) added to the Lanczos solution. Results are now also for larger values of x accurate to about 60 digits. Dynamic determination of digits and iterations added to the Spouge solution. It's now slightly faster on low x values and gives correct results using > 87 digits. Depending on the parameters, the program selects the optimal method for calculating Gamma.
 
<syntaxhighlight lang="rexx">
arg n; if n = '' then n = 100; numeric digits n
parse version version; say version; glob. = ''
say 'Gamma function to' n 'decimal places'
say
say '(Half)integers formulas'
w = '-99.5 -10.5 -5.5 -2.5 -1.5 -0.5 0.5 1 1.5 2 2.5 5 5.5 10 10.5 99 99.5'
numeric digits n
do i = 1 to words(w)
x = word(w,i); call time('r'); r = Gamma(x); e = format(time('e'),,3)
say 'Formulas' format(x,4,1) r '('e 'seconds)'
end
say
say 'Lanczos (max 60 decimals) vs Spouge (no limit) vs Stirling (no limit) approximation'
w = '-12.8 -6.4 -3.2 -1.6 -0.8 -0.4 -0.2 -0.1 0.1 0.2 0.4 0.8 1.6 3.2 6.4 12.8'
do i = 1 to words(w)
x = word(w,i)
numeric digits min(60,n)
call time('r'); r = Gamma(x); e = format(time('e'),,3)
say 'Lanczos ' format(x,4,1) r '('e 'seconds)'
numeric digits n
call time('r'); r = Gamma(x); e = format(time('e'),,3)
say 'Spouge ' format(x,4,1) r '('e 'seconds)'
if x > 0 then do
call time('r'); r = Stirling(x); e = format(time('e'),,3)
say 'Stirling' format(x,4,1) r '('e 'seconds)'
end
end
say
say 'Same for a bigger number'
w = '-99.9 99.9'
do i = 1 to words(w)
x = word(w,i)
numeric digits min(60,n)
call time('r'); r = Gamma(x); e = format(time('e'),,3)
say 'Lanczos ' format(x,4,1) r '('e 'seconds)'
numeric digits n
call time('r'); r = Gamma(x); e = format(time('e'),,3)
say 'Spouge ' format(x,4,1) r '('e 'seconds)'
if x > 0 then do
call time('r'); r = Stirling(x); e = format(time('e'),,3)
say 'Stirling' format(x,4,1) r '('e 'seconds)'
end
end
exit
 
Gamma:
/* Gamma */
procedure expose glob.
arg x
/* Validity */
if x < 1 & Whole(x) then
return 'X'
/* Formulas for negative and positive (half)integers */
if x < 0 then do
if Half(x) then do
numeric digits Digits()+2
i = Abs(Floor(x)); y = (-1)**i*2**(2*i)*Fact(i)*Sqrt(Pi())/Fact(2*i)
numeric digits Digits()-2
return y+0
end
end
if x > 0 then do
if Whole(x) then
return Fact(x-1)
if Half(x) then do
numeric digits Digits()+2
i = Floor(x); y = Fact(2*i)*Sqrt(Pi())/(2**(2*i)*Fact(i))
numeric digits Digits()-2
return y+0
end
end
p = Digits()
if p < 61 then do
/* Lanczos with predefined coefficients */
/* Map negative x to positive x */
if x < 0 then
return Pi()/(Gamma(1-x)*Sin(Pi()*x))
/* Argument reduction to interval (0.5,1.5) */
numeric digits p+2
c = Trunc(x); x = x-c
if x < 0.5 then do
x = x+1; c = c-1
end
/* Series coefficients 1/Gamma(x) in 80 digits Fransen & Wrigge */
c.1 = 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000
c.2 = 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467
c.3 = -0.65587807152025388107701951514539048127976638047858434729236244568387083835372210
c.4 = -0.04200263503409523552900393487542981871139450040110609352206581297618009687597599
c.5 = 0.16653861138229148950170079510210523571778150224717434057046890317899386605647425
c.6 = -0.04219773455554433674820830128918739130165268418982248637691887327545901118558900
c.7 = -0.00962197152787697356211492167234819897536294225211300210513886262731167351446074
c.8 = 0.00721894324666309954239501034044657270990480088023831800109478117362259497415854
c.9 = -0.00116516759185906511211397108401838866680933379538405744340750527562002584816653
c.10 = -0.00021524167411495097281572996305364780647824192337833875035026748908563946371678
c.11 = 0.00012805028238811618615319862632816432339489209969367721490054583804120355204347
c.12 = -0.00002013485478078823865568939142102181838229483329797911526116267090822918618897
c.13 = -0.00000125049348214267065734535947383309224232265562115395981534992315749121245561
c.14 = 0.00000113302723198169588237412962033074494332400483862107565429550539546040842730
c.15 = -0.00000020563384169776071034501541300205728365125790262933794534683172533245680371
c.16 = 0.00000000611609510448141581786249868285534286727586571971232086732402927723507435
c.17 = 0.00000000500200764446922293005566504805999130304461274249448171895337887737472132
c.18 = -0.00000000118127457048702014458812656543650557773875950493258759096189263169643391
c.19 = 0.00000000010434267116911005104915403323122501914007098231258121210871073927347588
c.20 = 0.00000000000778226343990507125404993731136077722606808618139293881943550732692987
c.21 = -0.00000000000369680561864220570818781587808576623657096345136099513648454655443000
c.22 = 0.00000000000051003702874544759790154813228632318027268860697076321173501048565735
c.23 = -0.00000000000002058326053566506783222429544855237419746091080810147188058196444349
c.24 = -0.00000000000000534812253942301798237001731872793994898971547812068211168095493211
c.25 = 0.00000000000000122677862823826079015889384662242242816545575045632136601135999606
c.26 = -0.00000000000000011812593016974587695137645868422978312115572918048478798375081233
c.27 = 0.00000000000000000118669225475160033257977724292867407108849407966482711074006109
c.28 = 0.00000000000000000141238065531803178155580394756670903708635075033452562564122263
c.29 = -0.00000000000000000022987456844353702065924785806336992602845059314190367014889830
c.30 = 0.00000000000000000001714406321927337433383963370267257066812656062517433174649858
c.31 = 0.00000000000000000000013373517304936931148647813951222680228750594717618947898583
c.32 = -0.00000000000000000000020542335517666727893250253513557337960820379352387364127301
c.33 = 0.00000000000000000000002736030048607999844831509904330982014865311695836363370165
c.34 = -0.00000000000000000000000173235644591051663905742845156477979906974910879499841377
c.35 = -0.00000000000000000000000002360619024499287287343450735427531007926413552145370486
c.36 = 0.00000000000000000000000001864982941717294430718413161878666898945868429073668232
c.37 = -0.00000000000000000000000000221809562420719720439971691362686037973177950067567580
c.38 = 0.00000000000000000000000000012977819749479936688244144863305941656194998646391332
c.39 = 0.00000000000000000000000000000118069747496652840622274541550997151855968463784158
c.40 = -0.00000000000000000000000000000112458434927708809029365467426143951211941179558301
c.41 = 0.00000000000000000000000000000012770851751408662039902066777511246477487720656005
c.42 = -0.00000000000000000000000000000000739145116961514082346128933010855282371056899245
c.43 = 0.00000000000000000000000000000000001134750257554215760954165259469306393008612196
c.44 = 0.00000000000000000000000000000000004639134641058722029944804907952228463057968680
c.45 = -0.00000000000000000000000000000000000534733681843919887507741819670989332090488591
c.46 = 0.00000000000000000000000000000000000032079959236133526228612372790827943910901464
c.47 = -0.00000000000000000000000000000000000000444582973655075688210159035212464363740144
c.48 = -0.00000000000000000000000000000000000000131117451888198871290105849438992219023663
c.49 = 0.00000000000000000000000000000000000000016470333525438138868182593279063941453996
c.50 = -0.00000000000000000000000000000000000000001056233178503581218600561071538285049997
c.51 = 0.00000000000000000000000000000000000000000026784429826430494783549630718908519485
c.52 = 0.00000000000000000000000000000000000000000002424715494851782689673032938370921241
/* Series expansion */
x = x-1; s = 0
do k = 52 by -1 to 1
s = s*x+c.k
end
y = 1/s
/* Undo reduction */
if c = -1 then
y = y/x
else do
do i = 1 to c
y = (x+i)*y
end
end
end
else do
x = x-1
/* Spouge */
/* Estimate digits and iterations */
q = Floor(p*1.5); a = Floor(p*1.3)
numeric digits q
/* Series */
s = 0
do k = 1 to a-1
s = s+((-1)**(k-1)*Power(a-k,k-0.5)*Exp(a-k))/(Fact(k-1)*(x+k))
end
s = s+Sqrt(2*Pi()); y = Power(x+a,x+0.5)*Exp(-a-x)*s
end
/* Normalize */
numeric digits p
return y+0
 
Stirling:
/* Sterling */
procedure expose glob.
arg x
return sqrt(2*pi()/x) * power(x/e(),x)
 
include Constants
include Functions
include Numbers
include Sequences
</syntaxhighlight>
{{out|Output:}}
<pre style="font-size:70%">
REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024
Gamma function in arbitrary precision
 
(Half)integers formulas (100 digits)
Formulas -99.5 3.370459273906717035419140191178165368212858243169804822385940657239537725368041400224226278392502618E-157 (0.003 seconds)
Formulas -10.5 -2.640121820547716316246385325311240439682468432522587656059168154777653141232089673307782521306919765E-7 (0.000 seconds)
Formulas -5.5 0.01091265478190986298673234429377905644050439299584730891829569009625650045347052040480970101310888490 (0.001 seconds)
Formulas -2.5 -0.9453087204829418812256893244486107641586930432652731350473641545882193517818838300666403502605571549 (0.000 seconds)
Formulas -1.5 2.363271801207354703064223311121526910396732608163182837618410386470548379454709575166600875651392887 (0.000 seconds)
Formulas -0.5 -3.544907701811032054596334966682290365595098912244774256427615579705822569182064362749901313477089331 (0.000 seconds)
Formulas 0.5 1.772453850905516027298167483341145182797549456122387128213807789852911284591032181374950656738544665 (0.000 seconds)
Formulas 1.0 1 (0.000 seconds)
Formulas 1.5 0.8862269254527580136490837416705725913987747280611935641069038949264556422955160906874753283692723327 (0.001 seconds)
Formulas 2.0 1 (0.000 seconds)
Formulas 2.5 1.329340388179137020473625612505858887098162092091790346160355842389683463443274136031212992553908499 (0.000 seconds)
Formulas 5.0 24 (0.000 seconds)
Formulas 5.5 52.34277778455352018114900849241819367949013237611424488006401129409378637307891910622901158181014715 (0.000 seconds)
Formulas 10.0 362880 (0.000 seconds)
Formulas 10.5 1133278.388948785567334574165588892475560298308275159776608723414529483390056004153717630538727607291 (0.000 seconds)
Formulas 99.0 9.426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729113E+153 (0.001 seconds)
Formulas 99.5 9.367802114655996591305637999137598430056780942876744878408158310966911153528806387157595583300106562E+154 (0.001 seconds)
 
Lanczos (60 digits) vs Spouge (61 digits) vs Stirling (61 digits) approximations
Lanczos -12.8 -1.44241809348812392704183496823417171228485786984071442421769E-9 (0.001 seconds)
Spouge -12.8 -1.442418093488123927041834968234171712284857869840714424217876E-9 (0.061 seconds)
Lanczos -6.4 -0.00214311842968855616971089012135323947575979728834356196087007 (0.001 seconds)
Spouge -6.4 -0.002143118429688556169710890121353239475759797288343561960870050 (0.062 seconds)
Lanczos -3.2 0.689056412005979742919224040168358601528149721171394266681422 (0.001 seconds)
Spouge -3.2 0.6890564120059797429192240401683586015281497211713942666814006 (0.061 seconds)
Lanczos -1.6 2.31058285808092523235318127282049942788600677267861414816871 (0.001 seconds)
Spouge -1.6 2.310582858080925232353181272820499427886006772678614148168727 (0.063 seconds)
Lanczos -0.8 -5.73855463999850381650594784491144000429263749786675377223601 (0.001 seconds)
Spouge -0.8 -5.738554639998503816505947844911440004292637497866753772236066 (0.063 seconds)
Lanczos -0.4 -3.72298062203204275598583347080335570330149759689981183834669 (0.001 seconds)
Spouge -0.4 -3.722980622032042755985833470803355703301497596899811838346699 (0.062 seconds)
Lanczos -0.2 -5.82114856862651686818160469134229346570980884445593876492447 (0.001 seconds)
Spouge -0.2 -5.821148568626516868181604691342293465709808844455938764924472 (0.062 seconds)
Lanczos -0.1 -10.6862870211931935489730533569448077816983878506097317904937 (0.001 seconds)
Spouge -0.1 -10.68628702119319354897305335694480778169838785060973179049371 (0.066 seconds)
Lanczos 0.1 9.51350769866873183629248717726540219255057862608837734305000 (0.000 seconds)
Spouge 0.1 9.513507698668731836292487177265402192550578626088377343050001 (0.063 seconds)
Stirling 0.1 5.697187148977169278607259516105973271247103273209125657453591 (0.003 seconds)
Lanczos 0.2 4.59084371199880305320475827592915200343410999829340301778885 (0.000 seconds)
Spouge 0.2 4.590843711998803053204758275929152003434109998293403017788853 (0.062 seconds)
Stirling 0.2 3.325998424022392556831252338044446631141936051050822874561565 (0.003 seconds)
Lanczos 0.4 2.21815954375768822305905402190767945077056650177146958224198 (0.000 seconds)
Spouge 0.4 2.218159543757688223059054021907679450770566501771469582241978 (0.063 seconds)
Stirling 0.4 1.841476335936235407774504215721792671787348602355998395147546 (0.003 seconds)
Lanczos 0.8 1.16422971372530337363632093826845869314196176889118775298489 (0.000 seconds)
Spouge 0.8 1.164229713725303373636320938268458693141961768891187752984894 (0.059 seconds)
Stirling 0.8 1.053370968425608533997094555812234627247042644411690511716880 (0.003 seconds)
Lanczos 1.6 0.893515349287690261436600032992805368792359423255954841203208 (0.000 seconds)
Spouge 1.6 0.8935153492876902614366000329928053687923594232559548412032077 (0.064 seconds)
Stirling 1.6 0.8486932421525737351870671884548739160413838607490704857265386 (0.003 seconds)
Lanczos 3.2 2.42396547993536801209211236969059225781321007909891679339251 (0.000 seconds)
Spouge 3.2 2.423965479935368012092112369690592257813210079098916793392514 (0.069 seconds)
Stirling 3.2 2.361851203186240411720280147420856619690317833030951977155556 (0.003 seconds)
Lanczos 6.4 240.833779983445938793193921422181728753410453111118999895588 (0.000 seconds)
Spouge 6.4 240.8337799834459387931939214221817287534104531111189998955875 (0.060 seconds)
Stirling 6.4 237.7207525902404787072183649425603774456984237805478743551168 (0.002 seconds)
Lanczos 12.8 289487660.334241583580309266192984495005040185968050409506654 (0.000 seconds)
Spouge 12.8 289487660.3342415835803092661929844950050401859680504095066538 (0.059 seconds)
Stirling 12.8 287609477.0875058075487675889481380625378498188858381901124308 (0.003 seconds)
 
Same for a bigger number
Lanczos -99.9 1.72726520939328005075554286550876635176440408500391962571907E-157 (0.000 seconds)
Spouge -99.9 1.727265209393280050755542865508766351764404085003870332156550733790282578981487471194267631617799447E-157 (0.274 seconds)
Lanczos 99.9 5.89173215164436165675854187059394609460949390669692020855752E+155 (0.000 seconds)
Spouge 99.9 5.891732151644361656758541870593946094609493906696920208557519945524581875362460658552892944423816715E+155 (0.266 seconds)
Stirling 99.9 5.886819525828908144161960074977041612792659661747365945060453973672990129231211036627906746201731112E+155 (0.010 seconds)
</pre>
 
Line 4,234 ⟶ 4,628:
:<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}}==
Line 4,279 ⟶ 4,703:
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>
 
Line 5,235 ⟶ 5,699:
{{libheader|Wren-math}}
The ''gamma'' method in the Math class is based on the Lanczos approximation.
<syntaxhighlight lang="ecmascriptwren">import "./fmt" for Fmt
import "./math" for Math
 
var stirling = Fn.new { |x| (2 * Num.pi / x).sqrt * (x / Math.e).pow(x) }
Line 5,243 ⟶ 5,707:
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))
System.write("%(Fmt.f(4, d, 2))\t")
System.write("%(Fmt.f(16, stirling.call(d), 14))\t")
System.print("%(Fmt.f(16, Math.gamma(d), 14))")
}</syntaxhighlight>
 
Line 5,254 ⟶ 5,716:
0.10 5.69718714897717 9.51350769866875
0.20 3.32599842402239 4.59084371199881
0.30 2.36253003626962 2.9915689876876099156898768759
0.40 1.84147633593624 2.21815954375769
0.50 1.52034690106628 1.77245385090552
Line 5,272 ⟶ 5,734:
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>