/* the functions int this file allow to calculate 
   the gamma and related functions */


#include <gamma.h>

double gamma_function(double z)
{
     double pi = 3.141592654;

     if(z > 0.)
        return exp(gammln(z));
     if(z < 0.)
        return pi / ( sin(pi * z) * exp(gammln(1. - z)) );   
}

/* the function gammln(xx) returns the natural logarithm of the gamma function 
   for positive argument xx */

double gammln(double xx)
{
     double x, y, tmp, ser;
     static double cof[6] = {76.18009172947146, -86.50532032941677,
			     24.01409824083091, -1.231739572450155,
			     0.1208650973866179E-2, -0.5395239384953E-5};
     int j;

     y = x = xx;
     tmp = x + 5.5;
     tmp -= (x + .5) * log(tmp);
     ser = 1.000000000190015;
	 for (j = 0; j<= 5; j++) ser += cof[j]/++y;

     return -tmp + log(2.5066282746310005 * ser / x);
}

/* the function factrl(n) returns the factorial of an integer n*/

double factrl(int n)
{
     static int ntop = 4;
     static double a[33] = {1.0,1.0,2.0,6.0,24.0};
     int j;

     if (n > 32) return exp(gammln( n + 1.0));
     else
     {
	while (ntop < n)
	{
	   j = ntop++;
	   a[ntop] = a[j] * ntop;
	}
     return a[n];
     }

}

/* the function bico(m,k) returns the binomial coefficient */
/* with m as upper entry and k as lower entry */

double bico(int m, int k)
{
     return floor(.5 + exp(log(factrl(m)) - log(factrl(k)) - log(factrl(m - k)) ));

}

/* the function beta(z, w) returns the beta function of z and w */

double beta(double z, double w)
{
     return exp(gammln(z) + gammln(w) - gammln(z + w));
}


/* the function art_gamma_1plusy() gives the argument of gamma(1 + i y) */

double arg_gamma_1plusy(double y)
{
     double sum = 0.0;
     int n;
     
     for(n = 0; n <= 100; n++)
     { 
       sum += y / (1. + n) - atan(y / (1. + n));
       
     }
     
     return sum - y * EULER;
     
}

dcomplex Cgamma_function(dcomplex z)
{
     double pi = 3.141592654;

     if(z.re > 0.)
        return Cexp(Cgammln(z));
     if(z.re < 0.)
        return CRmul(pi,Cinv(Cmul(Csin(CRmul(pi,z)),
		Cexp(Cgammln(Complex(1.-z.re,-z.im))))));   
}

/* the function Cgammln(xx) returns the natural logarithm of the gamma function 
   for positive argument xx */

dcomplex Cgammln(dcomplex xx)
{
     dcomplex x, y, tmp, ser, ser2;
     static double cof[6] = {76.18009172947146, -86.50532032941677,
			     24.01409824083091, -1.231739572450155,
			     0.1208650973866179E-2, -0.5395239384953E-5};
     int j;

     y = x = xx;
     tmp = Cadd(x,Complex(5.5,0.));
     tmp = Csub(tmp,Cmul(Cadd(x,Complex(.5,0.)),Clog(tmp)));
     ser = Complex(1.000000000190015,0.);
     for (j = 0; j<= 5; j++) 
	 {
		 y=Cadd(y,Complex(1.,0.));
		 ser2=CRmul(cof[j],Cinv(y));
		 ser=Cadd(ser,ser2);
	 }	 

     return Cadd(CRmul(-1.,tmp),
		 Clog(Cmul(CRmul(2.5066282746310005,ser),Cinv(x))));
	
}