/*
   GAMMA.I
   Gamma function, beta function, and binomial coefficients.

   $Id$
 */
/*    Copyright (c) 1994.  The Regents of the University of California.
                    All rights reserved.  */

/* Algorithms from Numerical Recipes, Press, et. al., section 6.1,
   Cambridge University Press (1988).
   Lanczos complex gamma approximation (lngamma) accuracy improved according
   to pgodfrey@intersil.com, http://winnie.fit.edu/~gabdo/gamma.txt
 */

func ln_gamma (z)
/* DOCUMENT ln_gamma(z)
     returns natural log of the gamma function.
     Error is less than 1e-13 for real part of z>=1.
     Use lngamma if you know that all z>=1, or if you don't care much about
     accuracy for z<1.
   SEE ALSO: lngamma, bico, beta
 */
{
  if (structof(z)==complex) big= (z.re>=1.0);
  else big= (z>=1.0);
  list= where(big);
  if (numberof(list)) lg1= lngamma(z(list));
  list= where(!big);
  if (numberof(list)) {
    z= z(list);
    lg2= log(pi/sin(pi*z)) - lngamma(1.0-z)
  }
  return merge(lg1,lg2,big);
}

func lngamma (x)
/* DOCUMENT lngamma(x)
     returns natural log of the gamma function.
     Error is less than 1e-13 for real part of x>=1.
     Use ln_gamma if some x<1.
   SEE ALSO: ln_gamma, bico, beta
 */
{
  /*
  ser = 1.0 + 76.18009173/x - 86.50532033/(x+1.) + 24.01409822/(x+2.) -
    1.231739516/(x+3.) + 0.120858003e-2/(x+4.) - 0.536382e-5/(x+5.);
  tmp = x+4.5;
  */
  ser = 1.000000000000000174663 +
    5716.400188274341379136/x +
    -14815.30426768413909044/(x+1.) +
    14291.49277657478554025/(x+2.) +
    -6348.160217641458813289/(x+3.) +
    1301.608286058321874105/(x+4.) +
    -108.1767053514369634679/(x+5.) +
    2.605696505611755827729/(x+6.) +
    -0.7423452510201416151527e-2/(x+7.) +
    0.5384136432509564062961e-7/(x+8.) +
    -0.4023533141268236372067e-8/(x+9.);
  tmp = x+8.5;
  return log(2.506628274631000502416*ser) + (x-0.5)*log(tmp) - tmp;
}

func bico (n,k)
/* DOCUMENT bico(n,k)
     returns the binomial coefficient n!/(k!(n-k)!) as a double.
   SEE ALSO: ln_gamma, beta
 */
{
  return floor(0.5+exp(lngamma(n+1.)-lngamma(k+1.)-lngamma(n-k+1.)));
}

func beta (z,w)
/* DOCUMENT beta(z,w)
     returns the beta function gamma(z)gamma(w)/gamma(z+w)
   SEE ALSO: ln_gamma, bico
 */
{
  return exp(lngamma(z)+lngamma(w)-lngamma(z+w));
}

/* see dawson.i for more accurate version, be careful not to clobber */
func erfc_nr (x)
/* DOCUMENT erfc_nr(x)
     returns the complementary error function 1-erf with fractional
     error less than 1.2e-7 everywhere.
 */
{
  z= abs(x);
  t= 1.0/(1.0+0.5*z);
  ans= t*exp(-z*z +
             poly(t, -1.26551223, 1.00002368, 0.37409196, 0.09678418,
                  -0.18628806, 0.27886807, -1.13520398, 1.48851587,
                  -0.82215223, 0.17087277));
  z= sign(x);
  return ans*z + (1.-z);
}
if (!is_func(erfc)) erfc = erfc_nr;