/*
 *  myGammaFunctions.cpp
 *  fp
 *
 *  Created by Robert Delaney on 4/24/09.
 *  Copyright 2009 __Bob Delaney's Science Software__. All rights reserved.
 *
 */

#include "fp.h"
#include "theBernoulliNumbers.h"

extern INT32 blockPrec;
// Gamma(x+1) is given by
// pow(x+a, (x+.5))*exp(-(x+a)) times
// 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(fp& z, const fp& 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 fp GammaSum(const fp& x, INT32 a)
{
	fp		z;
	
	GammaSum(z, x, a);
	
	return z;
	
}/* GammaSum, */



// uses Spouge formula for Gamma(x+1)
void GammaSpouge(fp& z, const fp& x)
{
	fp				xt;
	INT32			a, P;
	
	xt = x - 1.; // so we get Gamma(x)
	P = getDecPrec();
	
	a = (P-log(P)/log(10))*659/526 - .5; // change log to Lg2
	
	z = GammaSum(xt, a);
	
	z = z*pow(xt+a, (xt+.5))*exp(-(xt+a));
	
}/* GammaSpouge */


fp GammaSpouge(const fp& x)
{
	fp		z;
	
	GammaSpouge(z, x);
	
	return z;
}


// Euler's Constant calculation uses Gamma, so never calculate real arg Gamma using Euler's Constant!

// use reflection formula for negative arguments?
// for x+n>1 use Gamma(x+n) = x(x+1)(x+2)...(x+n-1)Gamma(x) to reduce to x<1
// if x is negative find integer n for which x+n is first positive
// then Gamma(x) = Gamma(x+n)/[x(x+1)(x+2)...(x+n-1)]
// use smallIncompleteGamma and largeIncompleteGamma to calculate the Gamma(x) or Gamma(x+n)

// I've reversed the meanings of x and a from previous work
// now x is the argument of the Gamma function and a is used in the integral definitions of the small and large Gamma functions
// the order of the arguments is the same as in the Handbook, but their a is my x and their x is my a
// uses first formula of 6.5.12 in "Handbook of Mathematical Functions"
void smallIncompleteGamma(fp& z, const fp& x, const fp& a)
{
	static fp		s, s1, t, t1, xinc;
	static bool		initGood=false;
	INT32			i;
	INT32			p;
	
	if(!initGood)
	{
		init(s);
		init(s1);
		init(t);
		init(t1);
		init(xinc);
		initGood = true;
	}
	
	if(a.i.n<=0 || x.i.n<=0)
	{
		z = 0;
		return;
	}
	
	p = getBlockPrec();
	setBlockPrec(p+1);
	
	equate(xinc, x);
	s = 0;
	s1 = 0;
	div(t, 1, x);
	
  	for (i=1; ; i++)
  	{
  		add(s1, s, t);
  		if(!compare(s, s1))
			break;
		equate(s, s1);
  		mul(t1, t, a);
		add(xinc, xinc, 1);
  		div(t, t1, xinc);
  	}
  	
  	//res = s*exp(-a + x*log(a));
	log(s1, a);
	mul(s1, s1, x);
	sub(s1, s1, a);
	exp(s1, s1);
	mul(z, s, s1);
	
	setBlockPrec(p);
  	
}/* smallIncompleteGamma */


// Essentially this uses the repeated fraction formula 6.5.31 in "Handbook of Mathematical Functions"
// The recursion relation in this routine carries out the repeated fraction.
void largeIncompleteGamma(fp& z, const fp& x, const fp& a)
{
	static fp		a0, a1, b0, b1, s, s1, s2;
	static bool		initGood=false;
	INT32			n;
	INT32			p;
	
	if(!initGood)
	{
		init(a0);
		init(a1);
		init(b0);
		init(b1);
		init(s);
		init(s1);
		init(s2);
		initGood = true;
	}
	
	if(a.i.n<=0)
	{
		z = 0;
		return;
	}
	
	p = getBlockPrec();
	setBlockPrec(2*p);
	
	a0 = 0;
	equate(b0, 1.);
	equate(a1, 1.);
	equate(b1, a);
	div(s, a1, b1);
	
	for(n=1;;n++)
	{
		//a0 = a1 + (n - x)*a0;
		sub (s1, n, x);
		mul(s1, s1, a0);
		add(a0, s1, a1);
		
		//b0 = b1 + (n - x)*b0;
		sub (s1, n, x);
		mul(s1, s1, b0);
		add(b0, s1, b1);
		
		//a1 = a*a0 + n*a1;
		mul(s1, n, a1);
		mul(s2, a, a0);
		add(a1, s1, s2);
		
		//b1 = a*b0 + n*b1;
		mul(s1, n, b1);
		mul(s2, a, b0);
		add(b1, s1, s2);
		
		//s1 = a1/b1;
		div(s1, a1, b1);
		setBlockPrec((3*p)/2);
		if (!compare(s, s1))
			break;
		setBlockPrec(2*p);
  		equate(s, s1);
	}
	
	//res = s*exp(-a + x*log(a));
	log(s1, a);
	mul(s1, s1, x);
	sub(s1, s1, a);
	exp(s1, s1);
	mul(z, s, s1);
	
	setBlockPrec(p);
	
}/* largeIncompleteGamma */


// only for x>=1
static void GammaSpecial(fp& z, const fp& x)
{
	fp			n, xt, a, gS, gL, one=1;
	INT32		nInt, i, p, p1;
	
	if(x<one)
		return;
	
	if(equate(nInt, x))
	{
		// x is the integer xInt so use factorial
		z = 1;
		for(i=1;i<nInt;++i)
			z = i*z;
		return;
	}
	
	// now x is not an integer, but it is > 1 and nInt>=1 is the truncated x
	
	xt = x - nInt + 1; // 1<xt<2
	// Gamma(x) = Gamma(xt+nInt-1) = xt(xt+1)(xt+2)...(xt+nInt-2)Gamma(xt)
	
	// we should consider increasing a with precision so that only smallIncompleteGamma is needed
	// integrand = t^(x-1)*e^(-t) ; for t = 300 and x=1.9 this is 8.7 x 10^(-129) so a good estimate
	// of the ratio of largeG to smallG (which is always about one for the x's given to this routine)
	// so for decPrec<=100 we only need smallG
	// we want the abs value of the exponent of 10 to be somewhat more than 8*blockPrec, say 10*blockPrec
	// try e^(-a) = 10^(-10*blockPrec)
	// or a = 23*blockPrec
	//a = 300.; // this appears to minimize time taken
	a = 20*getBlockPrec();  // was 23, 18 starts to go wrong, so use 20
	if(a<=300)
	{
		smallIncompleteGamma(z, xt, a);
	}
	else
	{
		a = 300;
		smallIncompleteGamma(gS, xt, a);
		p = getBlockPrec();
		p1 = p - 15; // largeGamma is about 10^(-128) * smallGamma so we can reduce the precision
		if(p1<2)
			p1 = 2;
		setBlockPrec(p1);
		largeIncompleteGamma(gL, xt, a);
		setBlockPrec(p);
		z = gS + gL;
	}
	
	for(i=0;i<nInt-1;++i)
		z = (xt + i)*z;
	
}/* GammaSpecial */


// need to use optimum a for given argument given to small and large incomplete functions
// for x+n>1 use Gamma(x+n) = x(x+1)(x+2)...(x+n-1)Gamma(x) to reduce to x<1
// might use a = xt-1 wich maximizes t^(xt-1)*e(-t) as a function of t, so want 1<xt<2 ( Gamma(xt) = (xt-1)!so argument of factorial is positive )
// xMax avoids overflow in pow function - dropped
//86100000 gives 6.529519196446844e645811020
//86200000 gives 7.064028851869211e-646309441
void Gamma(fp& z, const fp& x)
{
	fp			n, xt, a, gS, gL, one=1, xCrit=20000, xMax=86100000;
	INT32		nInt, i;
	INT32		p;
	
	if(x>xMax)
	{
		z = 0;
		return;
	}
	
	z = 0;  // z = 1 fouls up hyperF for x > .9 ! Why?
	
	if(x>=one)
	{
		if(x>xCrit)
			GammaSpouge(z, x);
		else
			GammaSpecial(z, x);
		return;
	}
	
	// now x<1 so bring argument up to >1 and use GammaSpecial
	if(equate(nInt, x))
		return; // have x = 0, -1, -2, ... where Gamma is unbounded.
	
	p = getBlockPrec();
	setBlockPrec(p+1);
	
	// Gamma(x) = Gamma(x+n)/[x(x+1)(x+2)...(x+n-1)]
	nInt = 2 - nInt; // nInt is one of 2, 3, ...
	xt = x + nInt; // xt>1
	GammaSpecial(z, xt);
	for(i=0;i<nInt;++i)
		z = z/(x+i);
	
	setBlockPrec(p);
	
}/* Gamma */


fp Gamma(const fp& x)
{
	fp		z;
	
	Gamma(z, x);
	
	return z;
	
}/* Gamma */


void Gamma(fp& z, double x)
{
	fp		xt;
	
	xt = x;
	
	Gamma(z, xt);
	
}/* Gamma */


fp Gamma(double x)
{
	fp		z;
	
	Gamma(z, x);
	
	return z;
	
}/* Gamma */


bool LnFactorial(fp& z, const fp& x)
{
	fp		ten, temp, temp1, lasttemp1;
	INT32	maxBernoulli=500;
	fp		thisTerm, zPower;
	fp		numRR, denRR;
	INT32	longProduct;
	INT32	n;
	bool	calcGood;
	INT32	saveNumBlocks;
	
	// temp = x
	if(!x.i.n || x==1.)
	{
		z = 0;
		return true;
	}
	
	ten = 10.;
	
	if(x<0)
	{
		z = 0;
		return false;
	}
	
	if(x>9*pow(ten, 10000000.))  // use Ln2
	{
		z = 0;
		return false;
	}
	
	if(x<=10000000)
	{
		Gamma(z, x+1);
		z = log(z);
		return true;
	}
	
	// now x>1e7 and we use the asymptotic series 6.1.40 on p. 257 of "Handbook of Mathematical Functions"
	// no longer check for too large a precision
	
	calcGood = true;
	saveNumBlocks = getBlockPrec();
	setBlockPrec(saveNumBlocks+2);
	
	temp = x + 1;
	temp1 = (temp-.5)*log(temp) - temp;
	
	temp1 = temp1 + .5*log(2*Pi());
	
	// if temp>10^10000 that's all we need
	
	static fp num[501], den[501];
	static bool	isGood=false;
	
	if(!isGood)
	{
		loadNumbers(num, den);
		isGood = true;
	}
	
	if(temp<=pow(ten, 10000.))
	{
		zPower = temp;  // z
		temp = 1/temp; // try to avoid overflow errors
		temp = temp*temp; // 1/z^2
		
		lasttemp1 = temp1;
		calcGood = false;
		for(n=1;n<=maxBernoulli;n++)
		{
			zPower = zPower*temp;
			longProduct = 2*n*(2*n-1);
			mul(thisTerm, zPower, num[n]);
			div(thisTerm, thisTerm, den[n]);
			thisTerm = thisTerm/longProduct;
			temp1 = temp1 + thisTerm;
			if(lasttemp1==temp1)
			{
				calcGood = true;
				break;
			}
			lasttemp1 = temp1;
		}
		// cout << "n = " << n << endl;
	}
	z = temp1;
	setBlockPrec(saveNumBlocks);
	
	return calcGood;
	
}/* LnFactorial */


fp LnFactorial(const fp& x)
{
	fp	z;
	
	LnFactorial(z, x);
	
	return z;
	
}/* LnFactorial */


char* TenToLogTenFactorial(const fp& x, bool& calcGood)
{
	char	*outString;
	fp		z;
	
	if(x<0)
	{
		outString = (char*)malloc(2*sizeof(char));
		outString[0] = '0';
		outString[1] = 0;
		return outString;
	}
	
	if(!x.i.n || x==1)
	{
		outString = (char*)malloc(2*sizeof(char));
		outString[0] = '1';
		outString[1] = 0;
		return outString;
	}
	
	calcGood = LnFactorial(z, x);
	
	outString = GetNumberFromLog(z);
	
	return outString;

}/* TenToLogTenFactorial */

