ncbi_math.c

ncbi源码

			}
			value += NCBIMATH_LNPI - log(sx);
		}
		else {
			y[0] = sin(x *= NCBIMATH_PI);
			tmp = 1.;
			for (k = 1; k <= order; k++) {
				tmp *= NCBIMATH_PI;
				y[k] = tmp * sin(x += (NCBIMATH_PI/2.));
			}
			value -= LogDerivative(order, y);
		}
	}
	else {
		value = general_lngamma(1. + x, order);
		if (order == 0) {
			if (x == 0.) {
/*
				ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
*/
				/**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
				return HUGE_VAL;
			}
			value -= log(x);
		}
		else {
			tmp = BLAST_Factorial(order - 1) * BLAST_Powi(x,  -order);
			value += (order % 2 == 0 ? tmp : - tmp);
		}
	}
	return value;
}

static double LnGamma(double x)               /* ln(ABS[gamma(x)]) - 10 dig
its of accuracy */
{
        return PolyGamma(x, 0);
}

#define FACTORIAL_PRECOMPUTED   36

extern double BLAST_Factorial(Int4 n)
{
        static double      precomputed[FACTORIAL_PRECOMPUTED]
                = { 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880., 3628800.};
        static Int4     nlim = 10;
        Int4   m;
        double    x;

        if (n >= 0) {
                if (n <= nlim)
                        return precomputed[n];
                if (n < DIM(precomputed)) {
                        for (x = precomputed[m = nlim]; m < n; ) {
                                ++m;
                                precomputed[m] = (x *= m);
                        }
                        nlim = m;
                        return x;
                }
                return exp(LnGamma((double)(n+1)));
        }
        return 0.0; /* Undefined! */
}

/* LnGammaInt(n) -- return log(Gamma(n)) for integral n */
extern double BLAST_LnGammaInt(Int4 n)
{
	static double	precomputed[FACTORIAL_PRECOMPUTED];
	static Int4	nlim = 1; /* first two entries are 0 */
	Int4 m;

	if (n >= 0) {
		if (n <= nlim)
			return precomputed[n];
		if (n < DIM(precomputed)) {
			for (m = nlim; m < n; ++m) {
				precomputed[m+1] = log(BLAST_Factorial(m));
			}
			return precomputed[nlim = m];
		}
	}
	return LnGamma((double)n);
}


/*
	Romberg numerical integrator

	Author:  Dr. John Spouge, NCBI
	Received:  July 17, 1992
	Reference:
		Francis Scheid (1968)
		Schaum's Outline Series
		Numerical Analysis, p. 115
		McGraw-Hill Book Company, New York
*/
#define F(x)  ((*f)((x), fargs))
#define ROMBERG_ITMAX 20

extern double BLAST_RombergIntegrate(double (*f) (double,void*), void* fargs, double p, double q, double eps, Int4 epsit, Int4 itmin)

{
	double	romb[ROMBERG_ITMAX];	/* present list of Romberg values */
	double	h;	/* mesh-size */
	Int4		i, j, k, npts;
	long	n;	/* 4^(error order in romb[i]) */
	Int4		epsit_cnt = 0, epsck;
	double	y;
	double	x;
	double	sum;

	/* itmin = min. no. of iterations to perform */
	itmin = MAX(1, itmin);
	itmin = MIN(itmin, ROMBERG_ITMAX-1);

	/* epsit = min. no. of consecutive iterations that must satisfy epsilon */
	epsit = MAX(epsit, 1); /* default = 1 */
	epsit = MIN(epsit, 3); /* if > 3, the problem needs more prior analysis */

	epsck = itmin - epsit;

	npts = 1;
	h = q - p;
	x = F(p);
	if (ABS(x) == HUGE_VAL)
		return x;
	y = F(q);
	if (ABS(y) == HUGE_VAL)
		return y;
	romb[0] = 0.5 * h * (x + y);	/* trapezoidal rule */
	for (i = 1; i < ROMBERG_ITMAX; ++i, npts *= 2, h *= 0.5) {
		sum = 0.;	/* sum of ordinates for */
		/* x = p+0.5*h, p+1.5*h, ..., q-0.5*h */
		for (k = 0, x = p+0.5*h; k < npts; k++, x += h) {
			y = F(x);
			if (ABS(y) == HUGE_VAL)
				return y;
			sum += y;
		}
		romb[i] = 0.5 * (romb[i-1] + h*sum); /* new trapezoidal estimate */

		/* Update Romberg array with new column */
		for (n = 4, j = i-1; j >= 0; n *= 4, --j)
			romb[j] = (n*romb[j+1] - romb[j]) / (n-1);

		if (i > epsck) {
			if (ABS(romb[1] - romb[0]) > eps * ABS(romb[0])) {
				epsit_cnt = 0;
				continue;
			}
			++epsit_cnt;
			if (i >= itmin && epsit_cnt >= epsit)
				return romb[0];
		}
	}
/*
	ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"RombergIntegrate: iterations > ROMBERG_ITMAX");
*/
	/**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > ROMBERG_ITMAX"));**/
	return HUGE_VAL;
}

/*
	Gcd(a, b)

	Return the greatest common divisor of a and b.

	Adapted 8-15-90 by WRG from code by S. Altschul.
*/
long BLAST_Gcd(long a, long b)
{
	long	c;

	b = ABS(b);
	if (b > a)
		c=a, a=b, b=c;

	while (b != 0) {
		c = a%b;
		a = b;
		b = c;
	}
	return a;
}

/* Round a floating point number to the nearest integer */
long BLAST_Nint(double x)	/* argument */
{
	x += (x >= 0. ? 0.5 : -0.5);
	return (long)x;
}

/*
integer power function

Original submission by John Spouge, 6/25/90
Added to shared library by WRG
*/
extern double BLAST_Powi(double x, Int4 n)	/* power */
{
	double	y;

	if (n == 0)
		return 1.;

	if (x == 0.) {
		if (n < 0) {
/*
			ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"Powi: divide by 0");
*/
			/**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
			return HUGE_VAL;
		}
		return 0.;
	}

	if (n < 0) {
		x = 1./x;
		n = -n;
	}

	while (n > 1) {
		if (n&1) {
			y = x;
			goto Loop2;
		}
		n /= 2;
		x *= x;
	}
	return x;

Loop2:
	n /= 2;
	x *= x;
	while (n > 1) {
		if (n&1)
			y *= x;
		n /= 2;
		x *= x;
	}
	return y * x;
}

extern double BLAST_LnFactorial (double x) {

    if(x<=0.0)
        return 0.0;
    else
        return LnGamma(x+1.0);
        
}

/*
 * ===========================================================================
 *
 * $Log: ncbi_math.c,v $
 * Revision 1000.2  2004/06/01 18:08:20  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
 *
 * Revision 1.8  2004/05/19 14:52:03  camacho
 * 1. Added doxygen tags to enable doxygen processing of algo/blast/core
 * 2. Standardized copyright, CVS $Id string, $Log and rcsid formatting and i
 *    location
 * 3. Added use of @todo doxygen keyword
 *
 * Revision 1.7  2003/12/05 16:03:57  camacho
 * Remove compiler warnings
 *
 * Revision 1.6  2003/09/26 20:39:32  dondosha
 * Rearranged code so it compiles
 *
 * Revision 1.5  2003/09/26 19:01:59  madden
 * Prefix ncbimath functions with BLAST_
 *
 * Revision 1.4  2003/09/10 21:36:29  dondosha
 * Removed Nlm_ prefix from math functions definitions
 *
 * Revision 1.3  2003/08/25 22:32:51  dondosha
 * Added #ifndef for definition of DBL_EPSILON
 *
 * Revision 1.2  2003/08/11 15:02:00  dondosha
 * Added algo/blast/core to all #included headers
 *
 * Revision 1.1  2003/08/02 16:31:48  camacho
 * Moved ncbimath.c -> ncbi_math.c
 *
 * Revision 1.1  2003/08/01 21:03:46  madden
 * Cleaned up version of file for C++ toolkit
 *
 * ===========================================================================
 */