/* Philip T.L.C. Clausen Jan 2017 plan@dtu.dk */

/*
 * Copyright (c) 2017, Philip Clausen, Technical University of Denmark
 * All rights reserved.
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *		http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
*/

#include <math.h>
#include "stdstat.h"

int (*cmp)(int, int) = &cmp_or;

int cmp_or(int t, int q) {
	return (t || q);
}

int cmp_and(int t, int q) {
	return (t && q);
}

double fastp(long double q) {
	/* P-value from quantile in a chi-square distribution */
	double p = 1.0;
	if(q > 114.5242) {
		p = 1e-26;
	} else if(q > 109.9604) {
		p = 1e-25;
	} else if(q > 105.3969) {
		p = 1e-24;
	} else if(q > 100.8337) {
		p = 1e-23;
	} else if(q > 96.27476) {
		p = 1e-22;
	} else if(q > 91.71701) {
		p = 1e-21;
	} else if(q > 87.16164) {
		p = 1e-20;
	} else if(q > 82.60901) {
		p = 1e-19;
	} else if(q > 78.05917) {
		p = 1e-18;
	} else if(q > 73.51245) {
		p = 1e-17;
	} else if(q > 68.96954) {
		p = 1e-16;
	} else if(q > 64.43048) {
		p = 1e-15;
	} else if(q > 59.89615) {
		p = 1e-14;
	} else if(q > 55.36699) {
		p = 1e-13;
	} else if(q > 50.84417) {
		p = 1e-12;
	} else if(q > 46.32844) {
		p = 1e-11;
	} else if(q > 41.82144) {
		p = 1e-10;
	} else if(q > 37.32489) {
		p = 1e-9;
	} else if(q > 32.84127) {
		p = 1e-8;
	} else if(q > 28.37395) {
		p = 1e-7;
	} else if(q > 23.92814) {
		p = 1e-6;
	} else if(q > 19.51139) {
		p = 1e-5;
	} else if(q > 15.13671) {
		p = 1e-4;
	} else if(q > 10.82759) {
		p = 1e-3;
	} else if(q > 6.634897) {
		p = 0.01;
	} else if(q > 3.841443) {
		p = 0.05;
	} else if(q > 2.705532) {
		p = 0.1;
	} else if(q > 2.072251) {
		p = 0.15;
	} else if(q > 1.642374) {
		p = 0.2;
	} else if(q > 1.323304) {
		p = 0.25;
	} else if(q > 1.074194) {
		p = 0.3;
	} else if(q > 0.8734571) {
		p = 0.35;
	} else if(q > 0.7083263) {
		p = 0.4;
	} else if(q > 0.5706519) {
		p = 0.45;
	} else if(q > 0.4549364) {
		p = 0.5;
	} else if(q > 0.3573172) {
		p = 0.55;
	} else if(q > 0.2749959) {
		p = 0.6;
	} else if(q > 0.2059001) {
		p = 0.65;
	} else if(q > 0.1484719) {
		p = 0.7;
	} else if(q > 0.1015310) {
		p = 0.75;
	} else if(q > 0.06418475) {
		p = 0.8;
	} else if(q > 0.03576578) {
		p = 0.85;
	} else if(q > 0.01579077) {
		p = 0.9;
	} else if(q > 0.00393214) {
		p = 0.95;
	} else if(q >= 0.0) {
		p = 1.0;
	} else {
		p = 1.00 - fastp(-1 * q);
	}
	return p;
}

double p_chisqr(long double q) {
	
	if(q < 0) {
		/* Handle negative overflow */
		return 1e-26;
	} else if(q > 49) {
		/* Get p-val from table, to avoid overflow */
		return fastp(q);
	}
	/* claculate p-value */
	return 1 - 1.772453850 * erf(sqrt(0.5 * q)) / tgamma(0.5);
}

double power(double x, unsigned n) {
	
	double y;
	
	if(n) {
		y = power(x, n >> 1);
		return (n & 1) ? y * y * x : y * y;
	}
	
	return 1.0;
}

double binP(int n, int k, double p) {
	
	int i, j, nk;
	double P, q, pq;
	
	/*
		P = n! / (k! (n-k)!) * p^k * q^(n - k)
	*/
	
	q = 1 - p;
	if(k == 0) {
		P = power(q, n);
		return P != 0.0 ? P : 1.0e-308;
	} else if(n == k) {
		P = power(p, n);
		return P != 0.0 ? P : 1.0e-308;
	} else if(p == 0 || q == 0) {
		return 0.0;
	}
	
	P = 1.0;
	nk = n - k;
	pq = p * q;
	i = n + 1;
	if(k < nk) {
		j = k + 1;
	} else {
		j = nk + 1;
	}
	
	while(--j) {
		P *= (--i * pq / j);
	}
	
	if(nk < k) {
		P *= power(p, k - nk);
	} else if(k < nk) {
		P *= power(q, nk - k);
	}
	
	return P != 0.0 ? P : 1.0e-308;
}

unsigned minimum(unsigned *src, unsigned n) {
	
	unsigned min;
	
	min = *src;
	while(--n) {
		if(*++src < min) {
			min = *src;
		}
	}
	
	return min;
}