Gamma function: Difference between revisions
m
syntax highlighting fixup automation
m (Reduced number of digits in coefficients to reflect float64 math.) |
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
||
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 496:
{{trans|BBC Basic}} - Lanczos method.
<
110
120 DEF FNgamma(z) = EXP(FNlngamma(z))
Line 521:
330 NEXT i
340 LET FNlngamma = (LOG(SQR(2*PI)) + LOG(a) - b) + LOG(b) * (z+0.5)
350 END FUNCTION</
=={{header|Arturo}}==
<
1.00000000000000000000 0.57721566490153286061 neg 0.65587807152025388108
neg 0.04200263503409523553 0.16653861138229148950 neg 0.04219773455554433675
Line 556:
pad (to :string v1-v2) 30
]
]</
{{out}}
Line 574:
{{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 655:
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 {
Line 677:
return exp(b*log(a))
}
</syntaxhighlight>
{{out}}
<pre>
Line 707:
{{trans|Phix}} - Lanczos method.
<
print
for i = 1 to 20
Line 735:
next i
return sqr(2.0 * pi) * (t ^ (x + 0.5)) * exp(-t) * a
end function</
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
Uses the Lanczos approximation.
<
INSTALL @lib$+"FNUSING"
Line 762:
a += lz(i%) / (z + i%)
NEXT
= (LNSQR(2*PI) + LN(a) - b) + LN(b) * (z+0.5)</
'''Output:'''
<pre>
Line 792:
This implements [[wp:Stirling's approximation|Stirling's approximation]] and [[wp:Spouge's approximation|Spouge's approximation]].
<
#include <stdlib.h>
#include <math.h>
Line 843:
}
return 0;
}</
{{out}}
<pre> Stirling Spouge GSL libm
Line 860:
=={{header|C sharp}}==
This is just rewritten from the Wikipedia Lanczos article. Works with complex numbers as well as reals.
<
using System.Numerics;
Line 887:
}
}
</syntaxhighlight>
=={{header|C++}}==
<
#include <numbers>
#include <stdio.h>
Line 948:
}
}
</syntaxhighlight>
{{out}}
<pre>
Line 974:
=={{header|Clojure}}==
<
"Returns Gamma(z + 1 = number) using Lanczos approximation."
[number]
Line 988:
(Math/exp (- (+ n 7 0.5)))
(+ (first c)
(apply + (map-indexed #(/ %2 (+ n %1 1)) (next c))))))))</
{{out}}
<
<pre>
0.1 9.5135
Line 1,025:
=={{header|Common Lisp}}==
<
(defconstant tcoeff
'( 1.00000000000000000000 0.57721566490153286061 -0.65587807152025388108
Line 1,050:
(loop for i from 1 to 10
do (
format t "~12,10f~%" (gamma (/ i 3.0))))</
{{out|Produces}}
<pre>
Line 1,068:
====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 1,110:
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 1,144:
=={{header|D}}==
<
real taylorGamma(in real x) pure nothrow @safe @nogc {
Line 1,205:
x.taylorGamma, x.lanczosGamma, x.gamma);
}
}</
{{out}}
<pre>0.333333: 2.6789385347077476335e+00 2.6789385347077470551e+00 2.6789385347077476339e+00
Line 1,221:
{{libheader| System.Math}}
{{Trans|Go}}
<syntaxhighlight lang="delphi">
program Gamma_function;
Line 1,247:
writeln(format('%5.1f %24.16g', [x, lanczos7(x)]));
readln;
end.</
{{out}}
<pre> x Lanczos7
Line 1,263:
=={{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,283:
Enum.each(Enum.map(1..10, &(&1/3)), fn x ->
:io.format "~f ~18.15f~n", [x, Gamma.taylor(x)]
end)</
{{out}}
Line 1,303:
Solved using the Lanczos Coefficients described in Numerical Recipes (Press et al)
<syntaxhighlight lang="f sharp">
open System
Line 1,320:
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,357:
=={{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,379:
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,387:
</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,411:
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,423:
{{works with|Fortran|2008}}
{{works with|Fortran|95 with extensions}}
<
implicit none
Line 1,520:
end function lacz_gamma
end program ComputeGammaInt</
{{out}}
<pre> Simpson Lanczos Builtin
Line 1,536:
=={{header|FreeBASIC}}==
{{trans|Java}}
<
Const pi = 3.1415926535897932
Line 1,582:
Print
Print "Press any key to quit"
Sleep</
{{out}}
Line 1,611:
=={{header|Go}}==
<
import (
Line 1,637:
1.5056327351493116e-7/(z+7)
return math.Sqrt2 * math.SqrtPi * math.Pow(t, z-.5) * math.Exp(-t) * x
}</
{{out}}
<pre>
Line 1,655:
=={{header|Groovy}}==
{{trans|Ada}}
<
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
Line 1,669:
(1..10).each{ printf("% 1.9e\n", gamma(it / 3.0)) }
</syntaxhighlight>
{{out}}
<pre> 2.678938535e+00
Line 1,685:
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,705:
main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]</
Or equivalently, as a point-free applicative expression:
<
cof :: [Double]
Line 1,729:
main :: IO ()
main = mapM_ print $ gammaln <$> [0.1,0.2 .. 1.0]</
{{Out}}
<pre>2.252712651734255
Line 1,746:
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,752:
procedure gamma(z) # Stirling's approximation
return (2*&pi/z)^0.5 * (z/&e)^z
end</
{{Out}}
Line 1,772:
=={{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).
<
<code><:</code> subtracts one from a number. It's sort of like <code>--lvalue</code> in C, except it always accepts an "rvalue" as an argument (which means it does not modify that argument). And <code>!value</code> finds the factorial of value if value is a positive integer. This illustrates the close relationship between the factorial and gamma functions.
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,820:
}
}
}</
{{out}}
<pre>
Line 1,848:
=={{header|JavaScript}}==
Implementation of Lanczos approximation.
<
var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
Line 1,867:
return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}</
=={{header|jq}}==
Line 1,873:
====Taylor Series====
{{trans|Ada}}
<
[
1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108, -0.04200263503409523553,
Line 1,888:
| 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,905:
($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,933:
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,992:
5.5 +5.234278e+01
=!EXPECTEND!=
*/</
{{out}}
Line 2,027:
'''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 2,040:
{{works with|Julia|1.0}}
'''Library function''':
<
gamma(1/2) - sqrt(pi)</
{{out}}
Line 2,047:
=={{header|Kotlin}}==
<
fun gammaStirling(x: Double): Double = Math.sqrt(2.0 * Math.PI / x) * Math.pow(x / Math.E, x)
Line 2,081:
println("%17.15f".format(gammaLanczos(d)))
}
}</
{{out}}
Line 2,114:
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 2,158:
return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;
}
</syntaxhighlight>
{{output}}
Line 2,176:
=={{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 2,186:
else return (z - 1) * gammafunc(z-1)
end
end</
=={{header|M2000 Interpreter}}==
<syntaxhighlight lang="m2000 interpreter">
Module PrepareLambdaFunctions {
Const e = 2.7182818284590452@
Line 2,229:
Print $("")
clipboard doc$
</syntaxhighlight>
χ Stirling Lanczos
0.10 5.697187148977170 9.5135076986687024462927178610
Line 2,254:
=={{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,265:
=={{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,312:
=={{header|Maxima}}==
<
gamma_coeff(n) := block([a: makelist(1, n)],
Line 2,332:
gamma_approx(12.3b0) - gamma(12.3b0);
/* -9.25224705314470500985141176997b-15 */</
=={{header|МК-61/52}}==
Line 2,345:
=={{header|Modula-3}}==
{{trans|Ada}}
<
FROM IO IMPORT Put;
Line 2,381:
Put(Extended(Taylor(FLOAT(i, EXTENDED) / 3.0X0), style := Style.Sci) & "\n");
END;
END Gamma.</
{{out}}
<pre>
Line 2,399:
{{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).
<
const A = [
Line 2,425:
let val1 = gamma(i.toFloat / 3)
let val2 = math.gamma(i.toFloat / 3)
echo &"{val1:18.16f} {val2:18.16f} {val1 - val2:11.4e}"</
{{Out}}
Line 2,442:
=={{header|OCaml}}==
<
let pi = 4. *. atan 1.
let sqrttwopi = sqrt (2. *. pi)
Line 2,504:
(Stirling.gamma z)
(Stirling2.gamma z)
done</
{{out}}
<pre>
Line 2,538:
=={{header|Octave}}==
<
p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \
771.32342877765313, -176.61502916214059, 12.507343278686905, \
Line 2,559:
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,575:
=={{header|Oforth}}==
<
[
Line 2,591:
2 Pi * sqrt *
t z 0.5 + powf *
t neg exp * ;</
{{out}}
Line 2,618:
=={{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}}==
Line 2,636:
Gamma(x)*Gamma(1 - x) = pi/sin(pi*x).
<
program GammaTest;
{$mode objfpc}{$H+}
Line 2,713:
WriteLn( SysUtils.Format( '-(32/945)Sqrt(pi) = %g', [-k*Sqrt(PI)]));
end.
</syntaxhighlight>
{{out}}
<pre>
Line 2,736:
=={{header|Perl}}==
<
use warnings;
use constant pi => 4*atan2(1, 1);
Line 2,805:
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,814:
=={{header|Phix}}==
{{trans|C}}
<!--<
<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>
Line 2,855:
<span style="color: #000000;">sq</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</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;">si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma</span><span style="color: #0000FF;">(</span><span style="color: #000000;">4</span><span style="color: #0000FF;">))</span>
<!--</
{{out}}
<small>(slightly different results under pwa/p2js...)</small>
Line 2,873:
Above translated to mpfr, with higher accuracy and more iterations as per REXX, and compared against the builtin.
{{libheader|Phix/mpfr}}
<!--<
<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: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.0"</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (mpfr_set_default_prec[ision] has been renamed)</span>
Line 2,969:
<span style="color: #000000;">sq</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma2</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;">si</span><span style="color: #0000FF;">(</span><span style="color: #000000;">gamma2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">4</span><span style="color: #0000FF;">))</span>
<!--</
{{out}}
<pre>
Line 2,996:
=={{header|Phixmonti}}==
<
-0.65587807152056 var coeff
-0.042002635033944 var quad
Line 3,036:
for
dup print " = " print gammafunc print nl
endfor</
=={{header|PicoLisp}}==
{{trans|Ada}}
<
(de *A
Line 3,059:
(for A (cdr *A)
(setq Sum (+ A (*/ Sum Y 1.0))) )
(*/ 1.0 1.0 Sum) ) )</
{{out}}
<pre>: (for I (range 1 10)
Line 3,075:
=={{header|PL/I}}==
<
test: procedure options (main);
Line 3,119:
end lanczos_gamma;
end test;</
{{out}}
<pre>
Line 3,137:
=={{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 3,169:
=={{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,
Line 3,206:
X2 is X1 + C*Power,
recip_gamma(X, Power, Cs, X1, X2, Y).
</syntaxhighlight>
{{Out}}
<pre>
Line 3,253:
* 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 3,334:
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 3,343:
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 3,358:
===Procedural===
{{trans|Ada}}
<
-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
Line 3,380:
for i in range(1,11):
print " %20.14e" % gamma(i/3.0)
</syntaxhighlight>
{{out}}
<pre> 2.67893853470775e+00
Line 3,396:
In terms of fold/reduce:
{{Works with|Python|3.7}}
<
from functools import reduce
Line 3,477:
# MAIN ---
if __name__ == '__main__':
main()</
{{Out}}
<pre> i -> gamma(i/3):
Line 3,495:
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,544:
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,559:
=={{header|Racket}}==
<
(define (gamma number)
(if (> 1/2 number)
Line 3,587:
; 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,612:
}
say Γ($_) for 1/3, 2/3 ... 10/3;</
{{out}}
<pre>2.67893853470775
Line 3,632:
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,699:
sum=sum*xm + #.k
end /*k*/
return 1/sum</
{{out|output|text= when using the input of: <tt> 0.5 </tt>}}
<pre>
Line 3,723:
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,791:
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,807:
=={{header|Ring}}==
<
decimals(3)
gamma = 0.577
Line 3,826:
but fabs(z) <= 0.5 return 1 / recigamma(z)
else return (z - 1) * gammafunc(z-1) ok
</syntaxhighlight>
=={{header|RLaB}}==
Line 3,841:
====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,857:
end
(1..10).each {|i| puts format("%.14e", gamma(i/3.0))}</
{{out}}
<pre>2.67893853470775e+00
Line 3,870:
2.77815847933857e+00</pre>
====Built in====
<
{{out}}
<pre>2.678938534707748
Line 3,885:
=={{header|Scala}}==
<
object Gamma {
Line 3,910:
println("%.1f -> %.16f %.16f".formatLocal(ENGLISH, x, stGamma(x), laGamma(x)))
}
}</
{{out}}
<pre>Gamma Stirling Lanczos
Line 3,939:
{{trans|Ruby}} for Taylor
<
(import (scheme base)
(scheme inexact)
Line 4,002:
(number->string (gamma-taylor i))
"\n")))
</syntaxhighlight>
{{out}}
Line 4,089:
=={{header|Scilab}}==
<syntaxhighlight lang="text">function x=gammal(z) // Lanczos'
lz=[ 1.000000000190015 ..
76.18009172947146 -86.50532032941677 24.01409824083091 ..
Line 4,110:
x=i/10
printf("%4.1f %9f %9f\n",x,gamma(x),gammal(x))
end</
{{out}}
<pre style="height:20ex">
Line 4,148:
=={{header|Seed7}}==
{{trans|Ada}}
<
include "float.s7i";
Line 4,189:
writeln((gamma(flt(I) / 3.0)) digits 15);
end for;
end func;</
{{out}}
<pre>2.678937911987305
Line 4,204:
=={{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 4,222:
for i in 1..10 {
say ("%.14e" % gamma(i/3))
}</
{{out}}
<pre>
Line 4,238:
Lanczos approximation:
<
var epsilon = 0.0000001
func withinepsilon(x) {
Line 4,274:
for i in 1..10 {
say ("%.14e" % gamma(i/3))
}</
{{out}}
Line 4,291:
A simpler implementation:
<
define τ = Num.tau
Line 4,301:
for i in (1..10) {
say ("%.14e" % Γ(i/3))
}</
{{out}}
Line 4,321:
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 4,381:
v=map(&gamma_(),x)
max(abs((v-u):/u))
end</
'''Output'''
Line 4,389:
Here follows the Python program to compute coefficients.
<
mp.dps = 50
Line 4,419:
for x in gc:
print(mp.nstr(x, 25))</
=={{header|Tcl}}==
Line 4,425:
{{tcllib|math}}
{{tcllib|math::calculus}}
<
package require math::calculus
Line 4,459:
set f2 [expr $f2 * $x]
}
}</
{{out}}
<pre> 1.0 1.000000 1.000000 1.000000
Line 4,486:
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,529:
===Lanczos' approximation===
<
I/2→X
e^(ln((1.0
Line 4,542:
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,581:
=={{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,642:
ENDDEFINE
</syntaxhighlight>
{{out}}
<pre>
Line 4,670:
=={{header|Vlang}}==
{{trans|go}}
<
fn main() {
println(" x math.Gamma Lanczos7")
Line 4,690:
1.5056327351493116e-7/(z+7)
return math.sqrt2 * math.sqrt_pi * math.pow(t, z-.5) * math.exp(-t) * x
}</
{{out}}
<pre> x math.Gamma Lanczos7
Line 4,709:
{{libheader|Wren-math}}
The ''gamma'' method in the Math class is based on the Lanczos approximation.
<
import "/math" for Math
Line 4,720:
System.write("%(Fmt.f(16, stirling.call(d), 14))\t")
System.print("%(Fmt.f(16, Math.gamma(d), 14))")
}</
{{out}}
Line 4,750:
=={{header|Yabasic}}==
{{trans|Phix}}
<
sub gamma(z)
Line 4,779:
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,805:
return(1.0 / sm);
}</
<
// Coefficients used by the GNU Scientific Library.
// http://en.wikipedia.org/wiki/Lanczos_approximation
Line 4,833:
return((2.0 * PI).sqrt() * t.pow(z + 0.5) * exp(-t) * x);
}
}</
{{out}}
<pre>
|