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

/* ehighs() --  return the expected high score */

Score_t
ehighs(Lambda, K, H, dblen, qlen)
	double	Lambda, K, H; /* parameters for Karlin statistics */
	unsigned long	dblen;	/* length of database */
	unsigned long	qlen;	/* length of query sequence */
{
	Score_t	highs, s;
	double	C, f, g, xqlen = qlen, xdblen = dblen;


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

	/* Calculate the expected high score */
	highs = log(C)/Lambda;
#ifdef REFINED_STATS
	do {
		f = MAX(xdblen - stolen(highs, Lambda, H), 1.) / xdblen;
		g = MAX(xqlen - stolen(highs, Lambda, H), 1.) / xqlen;
		s = log(C * f * g)/Lambda;
	} while (s < highs && (highs = s) > 0);
#endif

	return highs;
}