#include <ncbi.h>
#include "karlin.h"


/* etos() -- given an Expect value, return the associated cutoff score */
Score_t
etos(E, Lambda, K, H, dblen, qlen)
	double	E, Lambda, K, H; /* parameters for Karlin statistics */
	unsigned long	dblen;	/* length of database */
	unsigned long	qlen;	/* length of query sequence */
{
	register Score_t	S;
	double	C, f, g, x, prod1, prod2, xqlen = qlen, xdblen = dblen;
	double	elen;
	register Score_t	lo, hi;

	C = K * xqlen * xdblen; /* Frequently used constant */

	S = ceil( log(C/E) / Lambda );

#ifdef REFINED_STATS
	/* Binary search to find appropriate cutoff */
	hi = S;
	lo = S = 1;
	while (lo < hi) {
		S = lo + (hi - lo)/2;
		x = K * exp(-Lambda * S);
		elen = stolen(S, Lambda, H);
		f = MAX(xqlen - elen, 1.);
		g = MAX(xdblen - elen, 1.);
		prod1 = x*f*g;
		x = K * exp(-Lambda * (S+1));
		elen = stolen(S+1, Lambda, H);
		f = MAX(xqlen - elen, 1.);
		g = MAX(xdblen - elen, 1.);
		prod2 = x*f*g;
		if (prod1 <= E) {
			if (prod1 == E)
				break;
			hi = S;
		}
		else {
			lo = S+1;
			if (prod2 <= E) {
				S = lo;
				break;
			}
		}
	}
#endif /* REFINED_STATS */

	return S;
}