#version 3.1;
global_settings {assumed_gamma 1.0 }

// Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.

// Note that the functions Gamma and LogGamma are mutually dependent.

#macro Gamma (X) 

	// We require x > 0

	#if (X <= 0.0)
		#warning concat ("Invalid input argument ", str (X, 0, 4), ". Argument must be positive. \n")
		#break
	#end
	
	// Split the function domain into three intervals:
	// (0, 0.001), [0.001, 12), and (12, infinity)
	
	///////////////////////////////////////////////////////////////////////////
	// First interval: (0, 0.001)
	//
	// For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
	// So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
	// The relative error over this interval is less than 6e-7.
	
	#local _gamma = 0.577215664901532860606512090; // Euler's gamma constant
	
	#if (X < 0.001)
		#local Result = 1.0/(X*(1.0 + _gamma*X));
	#end
	
	///////////////////////////////////////////////////////////////////////////
	// Second interval: [0.001, 12)
	
	#if (X < 12.0)

		// The algorithm directly approximates gamma over (1,2) and uses
		// reduction identities to reduce other arguments to this interval.

		#local Y = X;
		#local n = 0;
		#local arg_was_less_than_one = (Y < 1.0);

		// Add or subtract integers as necessary to bring y into (1,2)
		// Will correct for this below
		#if (arg_was_less_than_one)
			#local Y = Y + 1.0;
		#else
			#local n = floor(Y) - 1;  // will use n later
			#local Y = Y - n;
		#end

		// numerator coefficients for approximation over the interval (1,2)
		#local p = array [8]
		{
			-1.71618513886549492533811E+0,
			 2.47656508055759199108314E+1,
			-3.79804256470945635097577E+2,
			 6.29331155312818442661052E+2,
			 8.66966202790413211295064E+2,
			-3.14512729688483675254357E+4,
			-3.61444134186911729807069E+4,
			 6.64561438202405440627855E+4
		};

		// denominator coefficients for approximation over the interval (1,2)
		#local q = array [8]
		{
			-3.08402300119738975254353E+1,
			 3.15350626979604161529144E+2,
			-1.01515636749021914166146E+3,
			-3.10777167157231109440444E+3,
			 2.25381184209801510330112E+4,
			 4.75584627752788110767815E+3,
			-1.34659959864969306392456E+5,
			-1.15132259675553483497211E+5
		};

		#local num = 0.0;
		#local den = 1.0;
		#local Z = Y - 1;
		
		#for (i, 0, 7)
			#local num = (num + p[i])*Z;
			#local den = den*Z + q[i];
		#end
		
		#local Result = num/den + 1.0;
		
		// Apply correction if argument was not initially in (1,2)
		#if (arg_was_less_than_one)
			// Use identity gamma(z) = gamma(z+1)/z
			// The variable "result" now holds gamma of the original y + 1
			// Thus we use y-1 to get back the orginal y.
			#local Result = Result / (Y-1.0);
		#else

			// Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
			#for (i, 0, n)
				// result *= y++;
				#local Result = Result * (Y+1);
			#end
		#end


		///////////////////////////////////////////////////////////////////////////
		// Third interval: [12, infinity)
		
		#if (X > 171.624)

		// Correct answer too large to display. Force +infinity.
			#local temp = 1e6;
			#local Result = temp*2.0;
			
			#local Result = exp(LogGamma(X));
		#end
		
		
		
		Result
		
#end // end macro Gamma
#end



#macro LogGamma (X) // x must be positive

	#if (X <= 0.0)
		#error concat ("Invalid input argument ", str (X, 0, 4), ". Argument must be positive. \n")
	#end

	#if (X < 12.0)
		#local Result = log( abs( Gamma (X)));
	#end
	
	// Abramowitz and Stegun 6.1.41
	// Asymptotic series should be good to at least 11 or 12 figures
	// For error analysis, see Whittiker and Watson
	// A Course in Modern Analysis (1927), page 252
	
	#local c = array [8]
	{
		 1.0/12.0,
		-1.0/360.0,
		 1.0/1260.0,
		-1.0/1680.0,
		 1.0/1188.0,
		-691.0/360360.0,
		 1.0/156.0,
		-3617.0/122400.0
	};
	
	#local Z = 1.0/(X*X);
	#local Sum = c[7];
	
	#for (i, 6, 0, -1)
		#local Sum = Sum * Z;
		#local Sum = Sum + c[i];
	#end
	
	#local series = Sum/X;

	#local halfLogTwoPi = 0.91893853320467274178032973640562;
	#local logGamma = (X - 0.5)*log(X) - X + halfLogTwoPi + series;    
	#local Result = logGamma;
	
	Result

#end // end LogGamma

#local Test = Gamma (0.5)

#debug concat( "Test = ", str(Test, 0, 16), "\n")

// Can delete these functions and the #include for <iostream> if not testing.

#macro TestGamma ()

	#local test = array [11][2] 
	{
		// Test near branches in code for (0, 0.001), [0.001, 12), (12, infinity)
		{1e-20, 1e+20},
		{2.19824158876e-16, 4.5490905327e+15},	// 0.99*DBL_EPSILON
		{2.24265050974e-16, 4.45900953205e+15}, // 1.01*DBL_EPSILON
		{0.00099, 1009.52477271},
		{0.00100, 999.423772485},
		{0.00101, 989.522792258},
		{6.1, 142.451944066},
		{11.999, 39819417.4793},
		{12, 39916800.0},
		{12.001, 40014424.1571},
		{15.2, 149037380723.0}
	};

	#local numTests = dimension_size (test, 1)-1;

	#local worst_absolute_error = 0.0;
	#local worst_relative_error = 0.0;
	#local worst_absolute_error_case = 0;
	#local worst_relative_error_case = 0;

	#for (T, 0, numTests)
		#local computed = Gamma (test[T][0])
		#local absolute_error = abs (computed - test[T][1]);
		#local relative_error = absolute_error / test[T][1];

		#if (absolute_error > worst_absolute_error)
			#local worst_absolute_error = absolute_error;
			#local worst_absolute_error_case = T;
		#end

		#if (relative_error > worst_relative_error)
			#local worst_relative_error = absolute_error;
			#local worst_relative_error_case = T;
		#end
	#end

	#local T = worst_absolute_error_case;
	#local X = test[T][0];
	#local Y = test[T][1];
	
	#debug concat (
		"Worst absolute error: ",
		str (abs(Gamma(X) - Y), 0, 16),
		"\nGamma ( ",
		str (X, 0, 4),
		") computed as " ,
		str (Gamma(X), 0, 4), 
		" but exact value is ",
		str (Y, 0, 4),
		"\n"
	)

	#local T = worst_relative_error_case;
	#local X = test[T][0];
	#local Y = test[T][1];
	
	#debug concat (
		"Worst relative error: ",
		str (abs(Gamma(X) - Y)/Y, 0, 16),
		"\nGamma ( ",
		str (X, 0, 4),
		") computed as " ,
		str (Gamma(X), 0, 4), 
		" but exact value is ",
		str (Y, 0, 4),
		"\n"
	)
#end // end macro 



TestGamma()


/*
void TestLogGamma()
{
	struct TestCase
	{
		double input;
		double expected;
	};

	TestCase test[] = 
	{
		{1e-12, 27.6310211159},
		{0.9999, 5.77297915613e-05},
		{1.0001, -5.77133422205e-05},
		{3.1, 0.787375083274},
		{6.3, 5.30734288962},
		{11.9999, 17.5020635801},
		{12, 17.5023078459},
		{12.0001, 17.5025521125},
		{27.4, 62.5755868211}
	};

	size_t numTests = sizeof(test) / sizeof(TestCase);

	double worst_absolute_error = 0.0;
	double worst_relative_error = 0.0;
	size_t worst_absolute_error_case = 0;
	size_t worst_relative_error_case = 0;

	for (size_t t = 0; t < numTests; t++)
	{
		double computed = LogGamma( test[t].input );
		double absolute_error = fabs(computed - test[t].expected);
		double relative_error = absolute_error / test[t].expected;

		if (absolute_error > worst_absolute_error)
		{
			worst_absolute_error = absolute_error;
			worst_absolute_error_case = t;
		}

		if (relative_error > worst_relative_error)
		{
			worst_relative_error = absolute_error;
			worst_relative_error_case = t;
		}
	}

	size_t t = worst_absolute_error_case;
	double x = test[t].input;
	double y = test[t].expected;
	std::cout << "Worst absolute error: "
			  << fabs(LogGamma(x) - y)
			  << "\nGamma( "
		      <<  x
			  << ") computed as " 
			  << LogGamma(x) 
			  << " but exact value is "
			  << y
			  << "\n";

	t = worst_relative_error_case;
	x = test[t].input;
	y = test[t].expected;
	std::cout << "Worst relative error: "
			  << (LogGamma(x) - y) / y
		      << "\nGamma( "
		      <<  x
			  << ") computed as " 
			  << LogGamma(x) 
			  << " but exact value is "
			  << y
			  << "\n";
}
*/



//TestLogGamma()