// $Id: WCMath.cpp,v 1.2 2011/01/07 04:21:52 samn Exp $ 

#include "WCMath.h"

// a lot of this code from :
/**
 * cvsid $Id: WCMath.cpp,v 1.2 2011/01/07 04:21:52 samn Exp $
 * @(#)Probability.java
 * Taken from Balasubramanian Narasimhan's Applet Demos for JSM 96 Talk
 * http://www-stat.stanford.edu/~naras/jsm/
 * @version 1.0 01/09/2002
 */


// Normal Distribution Cumulative Distribution function.
double normalCDF (double y) 
{
	double r[] = {
		1.253314137315500251207883,
			1.193182964731915311846094,1.137490921203604514832235,
			1.085827027468003637553896,1.037824575853726812300365,
			.9931557904881572182738326,.9515271920712067152786701,
			.9126755670832121676776603,.8763644564536923467278531,
			.8423810914521299997866721,.8105337152790304152036357,
			.7806492378708633711733062,.7525711790634080514554734,
			.7261578617139919103863031,.7012808218544300582109727,
			.6778234075911775329302582,.6556795424187984715438712,
			.6347526319769262711052108,.6149545961509297292775566,
			.5962050108690213064751457,.5784303460476310766336125,
			.5615632879362914156483528,.5455421356582169186076368,
			.5303102630712526699958747,.5158156382179633550265125,
			.5020103936204170322841234,.488850441527573754354743,
			.4762951289605100272316751,.4643069280394421644372609,
			.452851157630626619621641,.4418957328326000054087525,
			.4314109392400032535140663,.4213692292880544732249343,
			.4117450382989767017176231,.4025146181296716932830159,
			.3936558865630571131481576,.3851482907984342415801495,
			.3769726835829618014159496,.3691112106902638389489919,
			.3615472085963400644161054,.354265111329793337375764,
			.3472503655851968568924767,.3404893532870850071346728,
			.3339693208791823593693459,.3276783146905523474475952,
			.3216051217986084063997594,.3157392158694103491188545,
			.3100707075093597582907416,.304590298710103019232964,
			.2992892410108769288187367,.2941592970402895988586268,
			.2891927051332122944730683,.2843821467484925828542051,
			.2797207164400090390048569,.2752018941576061687414437,
			.2708195196759090041951434,.2665677689682234829510997,
			.2624411323600359552832388,.2584343943120382172559107,
			.2545426146965892806142785,.2507611114439652148255265,
			.2470854444460805703298742,.2435114006154562183725053,
			.2400349800063908654015814,.2366523829135604830915503,
			.233359997870698598664295,.2301543904788006381072361,
			.2270322929993801119871592,.2239905946538290863869083,
			.2210263325749769736284363,.2181366833614710122297952,
			.2153189551897363262218578,.2125705804420320545972856,
			.2098891088125368083169621,.2072722008565007589748572,
			.2047176219503302188107064,.2022232366330545215419131,
			.1997870033019862546952077,.1974069692375194490844627,
			.1950812659339918105426953,.1928081047153155616100713,
			.1905857726157402721374016,.1884126285076001815699527,
			.1862870994592909823499507,.1842076773079703381780956,
			.1821729154326492082990465,.1801814257143915475980612,
			.1782318756713315633990563,.1763229857571025870546387,
			.17445352681211268567823,.1726223176578507652332691,
			.170828222825113676267847,.1690701504076941880139726,
			.1673470500336553419899068,.1656579109468773975553577,
			.1640017601920640363366404,.1623776608968673374563811,
			.1607847106452193635505097,.1592220399363673265836744,
			.1576888107244717415286631,.1561842150339760769159709,
			.154707473646271358338463,.1532578348534790212960446,
			.1518345732754411325827937,.1504369887362691412736569,
			.1490644051970330214282687,.1477161697413932577413847,
			.1463916516111827416630501,.1450902412891308370578393,
			.1438113496261050352444228,.1425544070104022432905889,
			.1413188625767789529834851,.1401041834530502056148988,
			.1389098540422202650857274,.1377353753382303588480118,
			.1365802642735279803607799,.1354440530967635155710012,
			.1343262887790271962841785,.1332265324471292504019244,
			.132144358842515398166367,.1310793558044917945593207,
			.1300311237765102200812464,.128999275334337470011875,
			.1279834347349965694864678,.1269832374854368818978314,
			.1259983299299428153278645,.1250283688553503087964205,
			.1240730211131909094124509,.1231319632579322598468361,
			.1222048812005296276989127,.1212914698765461287919247,
			.1203914329281397110711456,.1195044823992530794107805,
			.1186303384433775019338515,.1177687290432979327966364,
			.1169193897422535131538919,.1160820633859823480155247,
			.1152564998751443189754465,.1144424559276431412284366,
			.1136396948503935650514457,.112847986320103044088645,
			.1120671061726592771214108,.1112968362007358601364944,
			.110536963959247939376068,.1097872825783083063597883,
			.1090475905833518904388312,.1083176917221130451719685,
			.1075973947981563778083775,.1068865135106744244241717,
			.1061848663002832119442509,.105492276200556096942253,
			.1048085706950511306815753,.1041335815795982142447371,
			.1034671448296236160489718,.1028091004723000198182652,
			.1021592924633203069380762,.1015175685681027854865852,
			.1008837802472445895049806,.1002577825460485147279854,
			.09963943398795665776104485,.09902859647173190962510087,
			.09842513517223564521296569,.09782891844465686959119527,
			.09723981773205465097029204,.09665770747608198672456013,
			.09608246503076461364872451,.09551397057921561379174917,
			.09495210705316911518666266,.09439676005522442883814663,
			.09384781778369499237091277,.09330517095996170483952101,
			.09276871275823452392594484,.09223833873763036123879576,
			.09171394677647927541850185,.09119543700877473684364846,
			.09068271176268733191388065,.09017567550106469857171747,
			.08967423476384374680079322,.08917829811230432667975505,
			.08868777607509647008527077,.08820258109597615778665339,
			.08772262748318725850402136,.0872478313604298571655505,
			.08677811061935764238292468,.08631338487354936400705558,
			.0858535754139016061424034,.08539860516539225449532534,
			.08494839864516607442904103,.08450288192189570024131742,
			.08406198257637363654606902,.08362562966329130783291586,
			.08319375367416516011749277,.0827662865013691388259681,
			.08234316140323582329192271,.08192431297018950814835363,
			.08150967709187601159489754,.08109919092525534990138904,
			.08069279286362471847049943,.0802904225065404650858022,
			.07989202063060893323638823,.07949752916111719512006792,
			.07910689114447578744239717,.0787200507214466106865758,
			.07833695310113015625477634,.07795754453568718779269906,
			.07758177229577092502292493,.07720958464664666235002043,
			.07684093082497660209183881,.07647576101624849508185813,
			.07611402633282746114038384,.07575567879261111001491475,
			.07540067129826880125610838,.07504895761704657047089306,
			.07470049236111991075842844,.07435523096847723310605399,
			.07401312968431743926073765,.07367414554294562620094298,
			.07333823635015150386570458,.0730053606660556482535154,
			.07267547778840923133747565,.07234854773633336836416045,
			.07202453123448470287828978,.07170338969763431106948403,
			.07138508521564745055865869,.07106958053885214609220417,
			.07075683906378472602241689,.07044682481930169454732137,
			.07013950245304622696078506,.06983483721825943539080115,
			.06953279496092599670031216,.06923334210724437899739519,
			.06893644565141218490800048,.06864207314371744366250278,
			.06835019267892698635104373,.06806077288496332988399061,
			.06777378291186177570382577,.06748919242099969956446575,
			.06720697157459026913544459,.06692709102543307719554012,
			.06664952190691442013000059,.06637423582325018469784634 };

		double f, h;
		int j;
		double  dcphi, x, z, f1, f2, f3, f4, f5;

		x = y;
		if (fabs(x) > 15.) {
			dcphi = 0.;
		} else {
			j = (int) floor(fabs(x) * 16. + .5);
			z = j * .0625;
			h = fabs(x) - z;
			f = r[j];
			f1 = f * z - 1;
			f2 = f + z * f1;
			f3 = f1 * 2. + z * f2;
			f4 = f2 * 3 + z * f3;
			f5 = f3 * 4 + z * f4;
			dcphi = f + h * (f1 * 120. + h * (f2 * 60. + h * (f3 * 20. + h * (f4 * 5. + h * f5)))) / 120.;
			dcphi = dcphi * .3989422804014326779 * exp(x * -.5 * x);
		}
		if (x < 0.) {
			return dcphi;
		} else {
			return (1.0 - dcphi);
		}
}

// From Numerical Recipes, with normal approximation from Appl. Stat. 239

static double EPSILON = 1.0e-14, LARGE_A = 10000.0;
static int ITMAX = 1000;

// Compute gamma cdf by a normal approximation
static double gnorm(double a, double x)
{
	double sx;

	if (x <= 0.0 || a <= 0.0)
		return 0.0;
	else {
		sx = sqrt(a) * 3.0 *
			(pow(x / a, 1.0 / 3.0) + 1.0 / (a * 9.0) - 1.0);
		return normalCDF(sx);
	}
}

// compute gamma cdf by its series representation
double gser(double a, double x, double gln) {
	double p, sum, del, ap;
	int n;
	bool done = false;

	if (x <= 0.0 || a <= 0.0)
		p = 0.0;
	else {
		ap = a;
		del = 1.0 / a;
		sum = del;
		for (n = 1; (! done) && (n < ITMAX); n++) {
			ap += 1.0;
			del *= x / ap;
			sum += del;
			if (fabs(del) < EPSILON)
				done = true;
		}
		p = sum * exp(- x + a * log(x) - gln);
	}
	return p;
}


// compute complementary gamma cdf by its continued fraction expansion
double gcf(double a, double x, double gln) {
	double gold = 0.0, g, fac = 1.0, b1 = 1.0;
	double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
	double p;
	bool done = false;

	a1 = x;
	p = 0.0;
	for (an = 1.0; (! done) && (an <= ITMAX); an += 1.0) {
		ana = an - a;
		a0 = (a1 + a0 * ana) * fac;
		b0 = (b1 + b0 * ana) * fac;
		anf = an * fac;
		a1 = x * a0 + anf * a1;
		b1 = x * b0 + anf * b1;
		if (a1 != 0.0) {
			fac = 1.0 / a1;
			g = b1 * fac;
			if (fabs((g - gold) / g) < EPSILON) {
				p = exp(-x + a * log(x) - gln) * g;
				done = true;
			}
			gold = g;
		}
	}
	return p;
}

double gammaCDF(double a, double x) {
	double gln;

	if (x <= 0.0 || a <= 0.0)
		return 0.0;
	else if (a > LARGE_A)
		return gnorm(a, x);
	else {
		gln = lngamma(a);
		if (x < (a + 1.0))
			return gser(a, x, gln);
		else
			return (1.0 - gcf(a, x, gln));
	}
}

double chisqCDF(double x, double df)
{
	return gammaCDF(0.5 * df, 0.5 * x);
}