/*
 *  GammaFunctionsC.cpp
 *  fp
 *
 *  Created by Robert Delaney on 11/13/09.
 *  Copyright 2009 Bob Delaney's Science Software. All rights reserved.
 *
 */

#include "fpComplex.h"

// Spouge formula sum from Wikipedia
// Sqrt(2*Pi) + Sum(k=1 to N) c(k)/(x+k)
// where N = a - 1 and
// c(0) = sqrt(2*Pi)
// c(k) = (-1)^(k-1) (a-k)^(k-.5) e^(a-k) / (k-1)!  k = 1. 2. ... , a-1
static void GammaSum(fpComplex& z, const fpComplex& x, INT32 a)
{
	INT32		k, p;
	fp			ck, den, temp, myE;
	
	p = getBlockPrec();
	setBlockPrec(1.4*p);
	
	exp(myE, 1.);
	
	// c(1)
	temp = a - 1;
	den = exp(-temp);
	z = sqrt(2*Pi());
	z = z + sqrt(temp)/((x+1.)*den);
	
	for(k=2;k<a;k++)
	{
		den = den*(k-1)*myE;
		temp = a - k;
		ck = pow(temp, k-.5)/den;
		if(k==2*(k/2))
			ck = -ck;
		z = z + ck/(x+k);
	}
	
	setBlockPrec(p);
	
}/* GammaSum */


static fpComplex GammaSum(const fpComplex& x, INT32 a)
{
	fpComplex		z;
	
	GammaSum(z, x, a);
	
	return z;
	
}/* GammaSum, */



// uses Spouge formula for Gamma(x+1)
static void GammaCalc(fpComplex& z, const fpComplex& x)
{
	fpComplex		xt;
	INT32			a, P, p;
	
	p = getBlockPrec();
	setBlockPrec(p+1);
	
	xt = x - 1.; // so we get Gamma(x)
	P = getDecPrec();
	
	a = (P-log(P)/log(10))*659/526 - .5; // change log to Lg2
	a = a+1;
	
	//cout << "a = " << a << endl; // a = 37 for decPrec = 32
	
	z = GammaSum(xt, a);
	
	z = z*pow(xt+a, (xt+.5))*exp(-(xt+a));
	
	setBlockPrec(p);
	
}/* GammaCalc */


void Gamma(fpComplex& z, const fpComplex& x)
{
	
	if(!x.im.i.n)
	{
		z.re = Gamma(x.re);
		z.im = 0;
		return;
	}
	
	if(x.re.i.n>=0)
	{
		GammaCalc(z, x);
		return;
	}
	
	// now use reflection formula
	fpComplex		xt, zt;
	xt = 1. - x;
	GammaCalc(zt, xt);
	
	z = Pi()/(zt*sin(Pi()*xt));
	
}/* Gamma */


fpComplex Gamma(const fpComplex& x)
{
	fpComplex		z;
	
	Gamma(z, x);
	
	return z;

}/* Gamma */