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