Gamma.frink

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, digits, debug] arbitraryExp[-k+a, digits, debug]
   }

   // 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["positiveGamma:  z is $z, a is $a"]
      // TODO:  Use cached sqrt of 2 pi
      return series[z-1, a, digits, debug] 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[max[50,ceil[(digits+3)*1.8]]]

         retval = 0
         if Re[z] > 0
            retval = positiveGamma[z, max[50, ceil[(digits+3)*1.8]], debug]
         else
         {
            pi = Pi.getPi[]
            // TODO:  Doh!  Implement arbitrary sin for complex numbers!
            retval = pi/(arbitrarySin[pi z, digits] 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 19944 days, 22 hours, 15 minutes ago.