Gamma function: Difference between revisions

* [[wp:Stirling's approximation|Stirling's approximation]]

<br><br>

 

=={{header|360 Assembly}}==

For maximum compatibility, this program uses only the basic instruction set.

USING GAMMAT,R13

SAVEAR B STM-SAVEAR(R15)

EQUREGS

EQUREGS REGS=FPR

{{out}}

<pre style="height:20ex">

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,

end loop;

return 1.0 / Sum;

Test program:

with Gamma;

 

Put_Line (Long_Float'Image (Gamma (Long_Float (I) / 3.0)));

end loop;

{{Out}}

<pre>

{{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,

FOR x FROM 10 BY 10 TO 70 DO

printf((sample exp fmt, x, sample(x)))

{{out}}

<pre>

gamma( 70)= 1.711224524e98, 1.711224524281e98, 1.711224524281e98, 7.57303907062e-29, 1.709188578191e98

</pre>

 

 

{{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.

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

 

{{works with|BBC BASIC for Windows}}

Uses the Lanczos approximation.

INSTALL @lib$+"FNUSING"

a += lz(i%) / (z + i%)

NEXT

'''Output:'''

<pre>

This implements [[wp:Stirling's approximation|Stirling's approximation]] and [[wp:Spouge's approximation|Spouge's approximation]].

 

#include <stdlib.h>

#include <math.h>

}

return 0;

{{out}}

<pre> Stirling Spouge GSL libm

2.70976382 2.77815848 2.77815848 2.77815848

</pre>

=={{header|C sharp}}==

This is just rewritten from the Wikipedia Lanczos article. Works with complex numbers as well as reals.

using System.Numerics;

 

}

}

 

=={{header|Clojure}}==

"Returns Gamma(z + 1 = number) using Lanczos approximation."

[number]

(Math/exp (- (+ n 7 0.5)))

(+ (first c)

{{out}}

<pre>

0.1 9.5135

 

=={{header|Common Lisp}}==

(defconstant tcoeff

'( 1.00000000000000000000 0.57721566490153286061 -0.65587807152025388108

(loop for i from 1 to 10

do (

{{out|Produces}}

<pre>

====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,

end

 

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)] }

{{out}}

<pre>

 

=={{header|D}}==

 

real taylorGamma(in real x) pure nothrow @safe @nogc {

x.taylorGamma, x.lanczosGamma, x.gamma);

}

{{out}}

<pre>0.333333: 2.6789385347077476335e+00 2.6789385347077470551e+00 2.6789385347077476339e+00

3.000000: 1.9999999999992207405e+00 2.0000000000000015575e+00 2.0000000000000000000e+00

3.333333: 2.7781584802531739378e+00 2.7781584804376666336e+00 2.7781584804376642124e+00</pre>

 

=={{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,

Enum.each(Enum.map(1..10, &(&1/3)), fn x ->

:io.format "~f ~18.15f~n", [x, Gamma.taylor(x)]

 

{{out}}

3.333333 2.778158479338573

</pre>

 

=={{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)

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

</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)

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

{{works with|Fortran|2008}}

{{works with|Fortran|95 with extensions}}

 

implicit none

end function lacz_gamma

 

{{out}}

<pre> Simpson Lanczos Builtin

=={{header|FreeBASIC}}==

{{trans|Java}}

 

Const pi = 3.1415926535897932

Print

Print "Press any key to quit"

 

{{out}}

1.90 0.920842721894229 0.961765831907388

2.00 0.959502175744492 1.000000000000000

</pre>

 

=={{header|Go}}==

 

import (

1.5056327351493116e-7/(z+7)

return math.Sqrt2 * math.SqrtPi * math.Pow(t, z-.5) * math.Exp(-t) * x

{{out}}

<pre>

=={{header|Groovy}}==

{{trans|Ada}}

-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,

-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,

 

(1..10).each{ printf("% 1.9e\n", gamma(it / 3.0)) }

{{out}}

<pre> 2.678938535e+00

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

 

main :: IO ()

 

Or equivalently, as a point-free applicative expression:

 

cof :: [Double]

 

main :: IO ()

{{Out}}

<pre>2.252712651734255

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

procedure gamma(z) # Stirling's approximation

return (2*&pi/z)^0.5 * (z/&e)^z

 

{{Out}}

=={{header|J}}==

This code shows the built-in method, which works for any value (positive, negative and complex numbers -- but note that negative integer arguments give infinite results).

 

The following direct coding of the task comes from the [[J:Essays/Stirling's%20Approximation|Stirling's approximation essay]] on the J wiki:

scorr =: 1 1r12 1r288 _139r51840 _571r2488320&p.@%

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

(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){

}

}

{{out}}

<pre>

=={{header|JavaScript}}==

Implementation of Lanczos approximation.

var p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,

771.32342877765313, -176.61502916214059, 12.507343278686905,

 

return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;

 

=={{header|jq}}==

====Taylor Series====

{{trans|Ada}}

[

1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108, -0.04200263503409523553,

| reduce range(2; 1+$n) as $an

($a[$n-1]; (. * $y) + $a[$n - $an])

====Lanczos Approximation====

def gamma_by_lanczos:

def pow(x): if x == 0 then 1 elif x == 1 then . else x * log | exp end;

($p[0]; . + ($p[$i]/($x + $i) ))

* ((2*$pi) | sqrt) * ($t | pow($x+0.5)) * ((-$t)|exp)

====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

====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

{{Out}}

i: gamma lanczos tgamma

0.3333333333333333: 2.6789385347077483 : 2.6789385347077483 : 2.678938534707748

2.6666666666666665: 1.5045754882515399 : 1.5045754882515576 : 1.5045754882515558

3 : 1.9999999999939684 : 2.0000000000000013 : 2

 

=={{header|Jsish}}==

{{trans|Javascript}}

/* Gamma function, in Jsish, using the Lanczos approximation */

function gamma(x) {

5.5 +5.234278e+01

=!EXPECTEND!=

 

{{out}}

 

'''Built-in function''':

 

'''By adaptive Gauss-Kronrod integration''':

gammaquad(t::Float64) = first(quadgk(x -> x ^ (t - 1) * exp(-x), zero(t), Inf, reltol = 100eps(t)))

 

{{out}}

{{works with|Julia|1.0}}

'''Library function''':

 

{{out}}

<pre>2.220446049250313e-16</pre>

 

=={{header|Kotlin}}==

 

fun gammaStirling(x: Double): Double = Math.sqrt(2.0 * Math.PI / x) * Math.pow(x / Math.E, x)

println("%17.15f".format(gammaLanczos(d)))

}

 

{{out}}

2.00 0.959502175744492 1.000000000000000

</pre>

 

=={{header|Limbo}}==

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;

return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;

}

 

{{output}}

=={{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

else return (z - 1) * gammafunc(z-1)

end

 

=={{header|M2000 Interpreter}}==

Module PrepareLambdaFunctions {

Const e = 2.7182818284590452@

Print $("")

clipboard doc$

χ Stirling Lanczos

0.10 5.697187148977170 9.5135076986687024462927178610

=={{header|Maple}}==

Built-in method that accepts any value.

GAMMA(7*I);

M := Matrix(2, 3, 'fill' = -3.6);

{{Out|Output}}

<pre>2027025*sqrt(Pi)*(1/256)

GAMMA(7*I)

Matrix(2, 3, [[.2468571430, .2468571430, .2468571430], [.2468571430, .2468571430, .2468571430]])</pre>

This code shows the built-in method, which works for any value (positive, negative and complex numbers).

Output integers and half-integers (a space is multiplication in Mathematica):

<pre>

 

=={{header|Maxima}}==

 

gamma_coeff(n) := block([a: makelist(1, n)],

 

gamma_approx(12.3b0) - gamma(12.3b0);

 

=={{header|Modula-3}}==

{{trans|Ada}}

 

FROM IO IMPORT Put;

Put(Extended(Taylor(FLOAT(i, EXTENDED) / 3.0X0), style := Style.Sci) & "\n");

END;

{{out}}

<pre>

1.9999999999939684e+000

2.7781584793385790e+000

</pre>

 

=={{header|Nim}}==

{{trans|Ada}}

1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108,

-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,

 

proc gamma(x: float): float =

 

for i in 1..10:

{{Out}}

 

=={{header|OCaml}}==

let pi = 4. *. atan 1.

let sqrttwopi = sqrt (2. *. pi)

(Stirling.gamma z)

(Stirling2.gamma z)

{{out}}

<pre>

 

=={{header|Octave}}==

p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \

771.32342877765313, -176.61502916214059, 12.507343278686905, \

for i = 1:10

printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));

{{out}}

<pre>2.678939 2.678939

=={{header|Oforth}}==

 

 

[

2 Pi * sqrt *

t z 0.5 + powf *

 

{{out}}

=={{header|PARI/GP}}==

===Built-in===

 

===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|Perl}}==

use warnings;

use constant pi => 4*atan2(1, 1);

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

stirling: 2.678938532866 1.354117938504 0.999999999306 0.892979510955 0.902745292336 0.999999999306 1.190639347940 1.504575487227 1.999999998611 2.778158478527</pre>

 

{{out}}

 

 

 

 

 

 

=={{header|PicoLisp}}==

{{trans|Ada}}

 

(de *A

(for A (cdr *A)

(setq Sum (+ A (*/ Sum Y 1.0))) )

{{out}}

<pre>: (for I (range 1 10)

 

=={{header|PL/I}}==

test: procedure options (main);

 

end lanczos_gamma;

 

{{out}}

<pre>

=={{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/

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)}

 

{{Out}}

0.961765831907388

1

</pre>

 

* 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)

Procedure Factorial(x) ; AKJ 01-May-10

ProcedureReturn Gamma(x+1)

;Examples

For i = 12 To 15

Debug StrD(i/3.0, 3)+" "+StrD(Gamma(i/3.0))

Debug "Ln(Gamma(5.0)) = "+StrD(GammLn(5.0), 16) ; Ln(24)

Debug ""

{{out}}

<pre>[Debug] Gamma():

===Procedural===

{{trans|Ada}}

-0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,

-0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,

for i in range(1,11):

print " %20.14e" % gamma(i/3.0)

{{out}}

<pre> 2.67893853470775e+00

In terms of fold/reduce:

{{Works with|Python|3.7}}

 

from functools import reduce

# MAIN ---

if __name__ == '__main__':

{{Out}}

<pre> i -> gamma(i/3):

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

{

z <- z - 1

tt <- z + g + 0.5

sqrt(2*pi) * tt^(z+0.5) * exp(-tt) * x

all.equal(as.complex(gamma(z)), lanczos(z)) # TRUE

all.equal(gamma(z), spouge(z)) # TRUE

{{out}}

z stirling nemes lanczos spouge builtin

 

=={{header|Racket}}==

(define (gamma number)

(if (> 1/2 number)

; 1.1642297137253037

; 1.068628702119319

 

=={{header|REXX}}==

 

Note: &nbsp; The Taylor series isn't much good above values of &nbsp; <big>'''6½'''</big>.

/*The GAMMA function symbol is the Greek capital letter: Γ */

numeric digits 90 /*be able to handle extended precision.*/

sum=sum*xm + #.k

end /*k*/

{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 0.5 </tt>}}

<pre>

 

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*/

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*/

{{out|output|text=&nbsp; when using the default input:}}

<pre>

 

=={{header|Ring}}==

decimals(3)

gamma = 0.577

but fabs(z) <= 0.5 return 1 / recigamma(z)

else return (z - 1) * gammafunc(z-1) ok

 

=={{header|RLaB}}==

:<math>RecGamma(a) = \frac{1}{\Gamma(a)}</math>;

:<math>Pochhammer(a,x) = \frac{\Gamma(a+x)}{\Gamma(x)}</math>.

 

=={{header|Ruby}}==

====Taylor series====

{{trans|Ada}}

-0.04200_26350_34095_23553, 0.16653_86113_82291_48950, -0.04219_77345_55544_33675,

-0.00962_19715_27876_97356, 0.00721_89432_46663_09954, -0.00116_51675_91859_06511,

end

 

{{out}}

<pre>2.67893853470775e+00

2.77815847933857e+00</pre>

====Built in====

{{out}}

<pre>2.678938534707748

2.0

2.7781584804376647

</pre>

 

=={{header|Scala}}==

 

object Gamma {

println("%.1f -> %.16f %.16f".formatLocal(ENGLISH, x, stGamma(x), laGamma(x)))

}

{{out}}

<pre>Gamma Stirling Lanczos

{{trans|Ruby}} for Taylor

 

(import (scheme base)

(scheme inexact)

(number->string (gamma-taylor i))

"\n")))

 

{{out}}

 

=={{header|Scilab}}==

lz=[ 1.000000000190015 ..

76.18009172947146 -86.50532032941677 24.01409824083091 ..

x=i/10

printf("%4.1f %9f %9f\n",x,gamma(x),gammal(x))

{{out}}

<pre style="height:20ex">

=={{header|Seed7}}==

{{trans|Ada}}

include "float.s7i";

 

writeln((gamma(flt(I) / 3.0)) digits 15);

end for;

{{out}}

<pre>2.678937911987305

=={{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,

for i in 1..10 {

say ("%.14e" % gamma(i/3))

{{out}}

<pre>

 

Lanczos approximation:

var epsilon = 0.0000001

func withinepsilon(x) {

for i in 1..10 {

say ("%.14e" % gamma(i/3))

 

{{out}}

 

A simpler implementation:

define τ = Num.tau

 

for i in (1..10) {

say ("%.14e" % Γ(i/3))

 

{{out}}

 

=={{header|Stata}}==

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,

 

function gamma_(x_) {

}

function map(f,a) {

}

 

u=map(&gamma(),x)

v=map(&gamma_(),x)

 

'''Output'''

 

 

=={{header|Tcl}}==

{{tcllib|math}}

{{tcllib|math::calculus}}

package require math::calculus

 

set f2 [expr $f2 * $x]

}

{{out}}

<pre> 1.0 1.000000 1.000000 1.000000

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

{{out}}

The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) .

 

===Lanczos' approximation===

I/2→X

e^(ln((1.0

Disp X,(X-1)!,Y

Pause

{{out}}

The output display for x=0.5 to 5 by 0.5 : x, (x-1)!, Y(x) .

24

</pre>

 

=={{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

 

ENDDEFINE

{{out}}

<pre>

1.9 0.961765831907391

2.0 1.000000000000000

</pre>

 

=={{header|Yabasic}}==

{{trans|Phix}}

sub gamma(z)

for i = 0.1 to 2.1 step .1

print i, " = "; : si(gamma(i))

 

=={{header|zkl}}==

{{trans|D}} but without a built in gamma function.

var table = T(

0x1p+0, 0x1.2788cfc6fb618f4cp-1,

 

return(1.0 / sm);

// Coefficients used by the GNU Scientific Library.

// http://en.wikipedia.org/wiki/Lanczos_approximation

return((2.0 * PI).sqrt() * t.pow(z + 0.5) * exp(-t) * x);

}

{{out}}

<pre>