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

/* etopoissonp -- convert an Expect value to a Poisson P-value */

double
etopoissonp(cnt, x)
	unsigned    cnt;
	register double  x;
{
	int		j;
	register unsigned	i;
	register unsigned	imin;
	register double	sum, prod, xi, tail;

	if (cnt < 2) {
		if (cnt == 1)
			return -fct_expm1(-x);
		return 1.;
	}

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

	sum = prod = fct_powi(x, cnt) / fct_factorial(cnt);
	/* 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/i);
		sum += prod;
	}

	sum = exp(-x) * (sum + tail);
	return MIN(sum, 1.);
}
#if 0 /* the old old way... */
#if 0 /* the old way... */
{
	register long    i; 
	register double  sum, prod;
	register double	p;

	sum = prod = 1.0;
	for (i = 1; i < cnt; ++i) {
		prod *= (x/i);
		sum += prod;  
	}
	if (x > 0.1*cnt)   /* Is limit approximation not applicable? */
		p = 1.0 - exp(-x)*sum;
	else
		p = prod * exp(-x) * x * (1. + x/(i+1.)) / i;

	return p;
}
#else /* the old new and improved way */
{
	register long	i;
	register double	sum, prod;
	register long	imax;

	/* Watch out when x is much greater than cnt--  SLOW!! */
	imax = x;
	imax = MAX(cnt + 10, imax + 15);

	sum = prod = fct_powi(x, cnt) / fct_factorial(cnt);
	for (i=cnt+1; i<imax; ++i) {
		prod *= (x/i);
		sum += prod;
	}

	return exp(-x) * sum;
}
#endif
#endif