Download or view Gamma.frink in plain text format
/* This calculates the Gamma function using Spouge's formula.
See:
http://en.literateprograms.org/Gamma_function_with_Spouge%27s_formula_%28Mathematica%29
and
http://en.literateprograms.org/Special:Downloadcode/Gamma_function_with_Spouge%27s_formula_%28Mathematica%29
Spouge's paper:
http://dx.doi.org/10.1137/0731050
For integers,
gamma[x] = (x-1)!
or
x! = gamma[x+1]
*/
use ArbitraryPrecision.frink
/**
This class has methods for calculating the Gamma function.
TODO: Calculate and store roots of pi.
*/
class Gamma
{
// Function to calculate the number of iterations required for convergence.
class gammaA[digits] := ceil[digits ln[10]/ln[2 pi]]
class coeff[k, a, digits, debug=false] :=
{
if debug
println["Entering coeff $k $a"]
if k == 0
return 1
else
return (-1)^(k-1)/(k-1)! arbitraryPow[-k+a, k-1/2] arbitraryExp[-k+a]
}
// TODO: Cache coefficients? They depend on a, though.
class series[z, a, digits, debug=false] :=
{
sum = 0
for k = a-1 to 1 step -1
{
term = coeff[k,a,digits,debug]/(z+k)
if debug
println["term $k is $term"]
sum = sum + term
}
if debug
println["About to return from series"]
return 1 + sum arbitraryPow[Pi.get2Pi[digits], -1/2, digits, debug]
}
// This is only valid for Re[z] > 0
class positiveGamma[z, digits, debug=false] :=
{
a = gammaA[digits]
if debug
println["a is $a"]
// TODO: Use cached sqrt of 2 pi
return series[z-1, a, digits] arbitraryPow[z-1+a, z-1/2, digits] arbitraryExp[-z+1-a, digits, debug] sqrt[Pi.get2Pi[digits], digits]
}
// Calculate gamma[z] to the specified number of digits.
class gamma[z, digits = getPrecision[], debug = false] :=
{
// Return the much faster factorial calculation if the number is an
// integer.
// THINK ABOUT: This may return a huge exact integer. Should
// this be comverted to a floating-point value first?
if isInteger[z] && z >= 1
return (z-1)!
origPrec = getPrecision[]
try
{
// TODO: This is too many digits, but needed for convergence. Find out why it's not
// converging with fewer digits.
setPrecision[ceil[(digits+3)*1.8]]
retval = 0
if Re[z] > 0
retval = positiveGamma[z, digits+10, debug]
else
{
pi = Pi.getPi[]
// TODO: Doh! Implement arbitrary sin for complex numbers!
retval = pi/(arbitrarySin[pi z] positiveGamma[1-z, digits, debug])
}
setPrecision[digits]
return 1. retval
}
finally
setPrecision[origPrec]
}
/** The beta function is closely tied to Gamma. */
class beta[a, b, digits=getPrecision[]] :=
{
origPrecision = getPrecision[]
try
{
workingPrec = digits + 4
setPrecision[workingPrec]
beta = gamma[a, workingPrec] gamma[b, workingPrec] / gamma[a+b, workingPrec]
setPrecision[digits]
return 1. * beta
}
finally
setPrecision[origPrecision]
}
}
/*n = 1000
digits = 20
println["\nResult:"]
println[" " + Gamma.gamma[n, digits, true]]
println[n]
origPrecision = getPrecision[]
try
{
setPrecision[digits]
if isInteger[n] and n>=0
println["Factorial: " + (n-1)! * 1.]
}
finally
setPrecision[origPrecision]
*/
/*
for n = 1 to 100
println[Gamma.gamma[n,12]]
for d=1 to 40
println[Gamma.gamma[1999.001, d]]
*/
/*
for d=1 to 40
println[Gamma.gamma[100.001, d]]
*/
Download or view Gamma.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19764 days, 4 hours, 7 minutes ago.