#include <ncbi.h>
#include <gishlib.h>
#include "blastapp.h"

static double R PROTO((unsigned i, unsigned j));

/* etoconsistp -- returns "consistent Poisson P-values" */
double
etoconsistp(cnt, x)
	unsigned	cnt; /* Number of consistent events */
	double	x;	/* Expected frequency for one event */
{
	register double	sum, prod, r, xi, xfact;
	register double	tail;
	Boolean		onezero = FALSE;
	unsigned	i, imin;
	int		j;

	if (cnt == 0)
		return 1.;
	if (cnt == 1)
		return -fct_expm1(-x);

	j = ceil(x);
	imin = MAX(j, 1) - 1;

	xfact = fct_factorial(cnt);
	sum = prod = fct_powi(x, cnt) / xfact;
	r = R(cnt, cnt);
	if (r != 0.)
		sum *= r / xfact;
	else
		onezero = TRUE;
	/* Using tail approximation suggested by Phil Green */
	for (xi = i = cnt+1; i <= imin || (tail = prod/(1.-x/xi)) > 0.01*sum; ++i, ++xi) {
		prod *= x / xi;
		if (!onezero) {
			r = R(i, cnt);
			if (r > 0.) {
				xfact *= xi;
				sum += prod * (r / xfact);
				continue;
			}
			onezero = TRUE;
		}
		sum += prod;
	}
	sum = exp(-x) * (sum + tail);
	return MIN(sum, 1.);
}


static double
R(i, j)
	unsigned	i, j;
{
	/*
	Precomputed values for R(i,4), for i=4,5,6,...

                            i
    R(i,4) = i! - 2(i!)^2  Sum ( (2k)! (3k^2 + 2k + 1 - 2ik)
                           k=0  --------------------------
                              k!^4 (i-k)!^2 (k+1)^2 (k+2)(i-k+1)
	*/
	static double	Ri4[] = {
			1.,
			17.,
			207.,
			2278.,
			24553.,
			268521.,
			3042210.,
			36153510.,
			454208895.,
			6059942223.,
			86030083110.,
			1299647574882.,
			20865826165777.,
			355277740280849.,
			6399391841784282.,
			121623163346687170.,
			2432739049821422100.,
			5.1089720946192155e+19,
			1.1239915020483751e+21,
			2.5851946529035947e+22,
			6.2044786379762287e+23,
			1.5511205895988435e+25,
			4.0329142896669799e+26,
			1.0888869199700814e+28,
			3.0488834264736675e+29,
			8.8417619782775417e+30,
			2.6525285968995213e+32,
			8.2228386532075972e+33,
			2.631308369259616e+35,
			8.6833176187500536e+36,
			2.9523279903910797e+38
		};
	register double	x, y, x2, x3;

	switch (j) {
	case 1:
		x = fct_factorial(i);
		if (x == HUGE_VAL)
			return 0.;
		return x;
	case 2:
		x = fct_factorial(i);
		if (x == HUGE_VAL)
			return 0.;
		return x - 1.;
	case 3:
		y = fct_factorial(2*i);
		if (y == HUGE_VAL)
			return 0.;
		x = fct_factorial(i);
		return x - y / (x * x * (i+1));
	case 4:
		if (i - 4 < DIM(Ri4))
			return Ri4[i - 4];
	default:
		break;
	}

	switch (i - j) {
	case 0:
		return 1.;
	case 1:
		x = j;
		return x * x + 1.;
	case 2:
		x = j;
		x2 = x * x;
		return (x2*x2 + x2) / 2. + x2*x + x + 3.;
	case 3:
		x = j;
		x2 = x * x;
		x3 = x2 * x;
		return (x3*x3 + 10.*x2*x2 + 4.*x3 + 19.*x2 + 62.*x) / 6. + x3*x2 + 11.;
	default:
		break;
	}

	return 0.;
}

#if 0
/* program to print out the R(i,4) terms, i >= 4 */
#include <ncbi.h>
#include <gishlib.h>

double	term(int i, int k);

main()
{
	int		i, k;
	double	x, val;

	for (i = 4; i < 35; ++i) {
		val = 0.;
		for (k = 0; k <= i; ++k)
			val += term(i, k);
		x = fct_factorial(i);
		val = x * (1. - 2. * x * val);
		printf("%d %.19lg\n", i, val);
	}
	exit(0);
}

double term(int i, int k)
{
	double	x, y;

	x = fct_factorial(2*k) * ((3.*k)*k + 2*k + 1 - i - (2.*i)*k);
	y = fct_factorial(i-k);
	x /= fct_powi(fct_factorial(k), 4) * y * y * (k+2) * (i-k+1);
	x /= (k+1);
	x /= (k+1);
	return x;
}
#endif