ncbi_math.c

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: ncbi_math.c,v $
 * PRODUCTION Revision 1000.2  2004/06/01 18:08:20  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
 * PRODUCTION
 * ===========================================================================
 */

/* $Id: ncbi_math.c,v 1000.2 2004/06/01 18:08:20 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Authors:  Gish, Kans, Ostell, Schuler
 *
 * Version Creation Date:   10/23/91
 *
 * ==========================================================================
 */

/** @file ncbi_math.c
 * Definitions for portable math library (ported from C Toolkit)
 * @todo FIXME doxygen comments and formatting
 */

static char const rcsid[] = 
    "$Id: ncbi_math.c,v 1000.2 2004/06/01 18:08:20 gouriano Exp $";

#define THIS_MODULE g_corelib
#define THIS_FILE _this_file

#include <algo/blast/core/ncbi_math.h>

#if 0
extern char * g_corelib;
static char * _this_file = __FILE__;
#endif


/*
    BLAST_Expm1(x)
    Return values accurate to approx. 16 digits for the quantity exp(x)-1
    for all x.
*/
extern double BLAST_Expm1(double	x)
{
  double	absx;

  if ((absx = ABS(x)) > .33)
    return exp(x) - 1.;

  if (absx < 1.e-16)
    return x;

  return x * (1. + x *
              (0.5 + x * (1./6. + x *
                          (1./24. + x * (1./120. + x *
                                         (1./720. + x * (1./5040. + x *
                                                         (1./40320. + x * (1./362880. + x *
                                                                           (1./3628800. + x * (1./39916800. + x *
                                                                                               (1./479001600. + x/6227020800.)
                                                                                               ))
                                                                           ))
                                                         ))
                                         ))
                          )));
}

#ifndef DBL_EPSILON
#define DBL_EPSILON 2.2204460492503131e-16
#endif

/*
    Log1p(x)
    Return accurate values for the quantity log(x+1) for all x > -1.
*/
extern double BLAST_Log1p(double x)
{
	Int4	i;
	double	sum, y;

	if (ABS(x) >= 0.2)
		return log(x+1.);

	/* Limit the loop to 500 terms. */
	for (i=0, sum=0., y=x; i<500 ; ) {
		sum += y/++i;
		if (ABS(y) < DBL_EPSILON)
			break;
		y *= x;
		sum -= y/++i;
		if (y < DBL_EPSILON)
			break;
		y *= x;
	}
	return sum;
}

static double LogDerivative(Int4 order, double* u) /* nth derivative of ln(u) */
        /* order is order of the derivative */
        /* u is values of u, u', u", etc. */
{
  Int4          i;
  double   y[LOGDERIV_ORDER_MAX+1];
  double  value, tmp;

  if (order < 0 || order > LOGDERIV_ORDER_MAX) {
InvalidOrder:
/*
    ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: unsupported derivative order");
*/
    /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
    return HUGE_VAL;
  }

  if (order > 0 && u[0] == 0.) {
/*
    ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: divide by 0");
*/
    /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
    return HUGE_VAL;
  }
  for (i = 1; i <= order; i++)
    y[i] = u[i] / u[0];

  switch (order) {
  case 0:
    if (u[0] > 0.)
      value = log(u[0]);
    else {
/*
      ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: log(x<=0)");
*/
      /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(x<=0)"));**/
      return HUGE_VAL;
    }
    break;
  case 1:
    value = y[1];
    break;
  case 2:
    value = y[2] - y[1] * y[1];
    break;
  case 3:
    value = y[3] - 3. * y[2] * y[1] + 2. * y[1] * y[1] * y[1];
    break;
  case 4:
    value = y[4] - 4. * y[3] * y[1] - 3. * y[2] * y[2]
      + 12. * y[2] * (tmp = y[1] * y[1]);
    value -= 6. * tmp * tmp;
    break;
  default:
    goto InvalidOrder;
  }
  return value;
}

static double _default_gamma_coef [] = {
         4.694580336184385e+04,
        -1.560605207784446e+05,
         2.065049568014106e+05,
        -1.388934775095388e+05,
         5.031796415085709e+04,
        -9.601592329182778e+03,
         8.785855930895250e+02,
        -3.155153906098611e+01,
         2.908143421162229e-01,
        -2.319827630494973e-04,
         1.251639670050933e-10
         };

static int      xgamma_dim = DIM(_default_gamma_coef);

static double
general_lngamma(double x, Int4 order)      /* nth derivative of ln[gamma(x)] */
        /* x is 10-digit accuracy achieved only for 1 <= x */
        /* order is order of the derivative, 0..POLYGAMMA_ORDER_MAX */
{
        Int4            i;
        double     xx, tx;
        double     y[POLYGAMMA_ORDER_MAX+1];
        double    tmp, value;
	double    *coef;

        xx = x - 1.;  /* normalize from gamma(x + 1) to xx! */
        tx = xx + xgamma_dim;
        for (i = 0; i <= order; ++i) {
                tmp = tx;
                /* sum the least significant terms first */
                coef = &_default_gamma_coef[xgamma_dim];
                if (i == 0) {
                        value = *--coef / tmp;
                        while (coef > _default_gamma_coef)
                                value += *--coef / --tmp;
                }
                else {
                        value = *--coef / BLAST_Powi(tmp, i + 1);
                        while (coef > _default_gamma_coef)
                                value += *--coef / BLAST_Powi(--tmp, i + 1);
                        tmp = BLAST_Factorial(i);
                        value *= (i%2 == 0 ? tmp : -tmp);
                }
                y[i] = value;
        }
        ++y[0];
        value = LogDerivative(order, y);
        tmp = tx + 0.5;
        switch (order) {
        case 0:
                value += ((NCBIMATH_LNPI+NCBIMATH_LN2) / 2.)
                                        + (xx + 0.5) * log(tmp) - tmp;
                break;
        case 1:
                value += log(tmp) - xgamma_dim / tmp;
                break;
        case 2:
                value += (tmp + xgamma_dim) / (tmp * tmp);
                break;
        case 3:
                value -= (1. + 2.*xgamma_dim / tmp) / (tmp * tmp);
                break;
        case 4:
                value += 2. * (1. + 3.*xgamma_dim / tmp) / (tmp * tmp * tmp);
                break;
        default:
                tmp = BLAST_Factorial(order - 2) * BLAST_Powi(tmp, 1 - order)
                                * (1. + (order - 1) * xgamma_dim / tmp);
                if (order % 2 == 0)
                        value += tmp;
                else
                        value -= tmp;
                break;
        }
        return value;
}


static double PolyGamma(double x, Int4 order) /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
	/* x is and derivatives */
	/* order is order of the derivative */
/* order = 0, 1, 2, ...  ln(gamma), digamma, trigamma, ... */
/* CAUTION:  the value of order is one less than the suggested "di" and
"tri" prefixes of digamma, trigamma, etc.  In other words, the value
of order is truly the order of the derivative.  */
{
	Int4		k;
	double	value, tmp;
	double	y[POLYGAMMA_ORDER_MAX+1], sx;

	if (order < 0 || order > POLYGAMMA_ORDER_MAX) {
/*
		ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: unsupported derivative order");
*/
		/**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
		return HUGE_VAL;
	}

	if (x >= 1.)
		return general_lngamma(x, order);

	if (x < 0.) {
		value = general_lngamma(1. - x, order);
		value = ((order - 1) % 2 == 0 ? value : -value);
		if (order == 0) {
			sx = sin(NCBIMATH_PI * x);
			sx = ABS(sx);
			if ( (x < -0.1 && (ceil(x) == x || sx < 2.*DBL_EPSILON))
					|| sx == 0.) {
/*
				ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
*/
				/**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
				return HUGE_VAL;