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 19574 days, 10 hours, 42 minutes ago.