开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
ncbi_math.c
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：16k
源码类别：

生物技术

开发平台：

C/C++

ncbi_math.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: ncbi_math.c,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:08:20  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: ncbi_math.c,v 1000.2 2004/06/01 18:08:20 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's official duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Authors:  Gish, Kans, Ostell, Schuler
  36.  *
  37.  * Version Creation Date:   10/23/91
  38.  *
  39.  * ==========================================================================
  40.  */
  42. /** @file ncbi_math.c
  43.  * Definitions for portable math library (ported from C Toolkit)
  44.  * @todo FIXME doxygen comments and formatting
  45.  */
  47. static char const rcsid[] = 
  48.     "$Id: ncbi_math.c,v 1000.2 2004/06/01 18:08:20 gouriano Exp $";
  50. #define THIS_MODULE g_corelib
  51. #define THIS_FILE _this_file
  53. #include <algo/blast/core/ncbi_math.h>
  55. #if 0
  56. extern char * g_corelib;
  57. static char * _this_file = __FILE__;
  58. #endif
  61. /*
  62.     BLAST_Expm1(x)
  63.     Return values accurate to approx. 16 digits for the quantity exp(x)-1
  64.     for all x.
  65. */
  66. extern double BLAST_Expm1(double x)
  67. {
  68.   double absx;
  70.   if ((absx = ABS(x)) > .33)
  71.     return exp(x) - 1.;
  73.   if (absx < 1.e-16)
  74.     return x;
  76.   return x * (1. + x *
  77.               (0.5 + x * (1./6. + x *
  78.                           (1./24. + x * (1./120. + x *
  79.                                          (1./720. + x * (1./5040. + x *
  80.                                                          (1./40320. + x * (1./362880. + x *
  81.                                                                            (1./3628800. + x * (1./39916800. + x *
  82.                                                                                                (1./479001600. + x/6227020800.)
  83.                                                                                                ))
  84.                                                                            ))
  85.                                                          ))
  86.                                          ))
  87.                           )));
  88. }
  90. #ifndef DBL_EPSILON
  91. #define DBL_EPSILON 2.2204460492503131e-16
  92. #endif
  94. /*
  95.     Log1p(x)
  96.     Return accurate values for the quantity log(x+1) for all x > -1.
  97. */
  98. extern double BLAST_Log1p(double x)
  99. {
  100. Int4 i;
  101. double sum, y;
  103. if (ABS(x) >= 0.2)
  104. return log(x+1.);
  106. /* Limit the loop to 500 terms. */
  107. for (i=0, sum=0., y=x; i<500 ; ) {
  108. sum += y/++i;
  109. if (ABS(y) < DBL_EPSILON)
  110. break;
  111. y *= x;
  112. sum -= y/++i;
  113. if (y < DBL_EPSILON)
  114. break;
  115. y *= x;
  116. }
  117. return sum;
  118. }
  120. static double LogDerivative(Int4 order, double* u) /* nth derivative of ln(u) */
  121.         /* order is order of the derivative */
  122.         /* u is values of u, u', u", etc. */
  123. {
  124.   Int4          i;
  125.   double   y[LOGDERIV_ORDER_MAX+1];
  126.   double  value, tmp;
  128.   if (order < 0 || order > LOGDERIV_ORDER_MAX) {
  129. InvalidOrder:
  130. /*
  131.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: unsupported derivative order");
  132. */
  133.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  134.     return HUGE_VAL;
  135.   }
  137.   if (order > 0 && u[0] == 0.) {
  138. /*
  139.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: divide by 0");
  140. */
  141.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  142.     return HUGE_VAL;
  143.   }
  144.   for (i = 1; i <= order; i++)
  145.     y[i] = u[i] / u[0];
  147.   switch (order) {
  148.   case 0:
  149.     if (u[0] > 0.)
  150.       value = log(u[0]);
  151.     else {
  152. /*
  153.       ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: log(x<=0)");
  154. */
  155.       /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(x<=0)"));**/
  156.       return HUGE_VAL;
  157.     }
  158.     break;
  159.   case 1:
  160.     value = y[1];
  161.     break;
  162.   case 2:
  163.     value = y[2] - y[1] * y[1];
  164.     break;
  165.   case 3:
  166.     value = y[3] - 3. * y[2] * y[1] + 2. * y[1] * y[1] * y[1];
  167.     break;
  168.   case 4:
  169.     value = y[4] - 4. * y[3] * y[1] - 3. * y[2] * y[2]
  170.       + 12. * y[2] * (tmp = y[1] * y[1]);
  171.     value -= 6. * tmp * tmp;
  172.     break;
  173.   default:
  174.     goto InvalidOrder;
  175.   }
  176.   return value;
  177. }
  179. static double _default_gamma_coef [] = {
  180.          4.694580336184385e+04,
  181.         -1.560605207784446e+05,
  182.          2.065049568014106e+05,
  183.         -1.388934775095388e+05,
  184.          5.031796415085709e+04,
  185.         -9.601592329182778e+03,
  186.          8.785855930895250e+02,
  187.         -3.155153906098611e+01,
  188.          2.908143421162229e-01,
  189.         -2.319827630494973e-04,
  190.          1.251639670050933e-10
  191.          };
  193. static int      xgamma_dim = DIM(_default_gamma_coef);
  195. static double
  196. general_lngamma(double x, Int4 order)      /* nth derivative of ln[gamma(x)] */
  197.         /* x is 10-digit accuracy achieved only for 1 <= x */
  198.         /* order is order of the derivative, 0..POLYGAMMA_ORDER_MAX */
  199. {
  200.         Int4            i;
  201.         double     xx, tx;
  202.         double     y[POLYGAMMA_ORDER_MAX+1];
  203.         double    tmp, value;
  204. double    *coef;
  206.         xx = x - 1.;  /* normalize from gamma(x + 1) to xx! */
  207.         tx = xx + xgamma_dim;
  208.         for (i = 0; i <= order; ++i) {
  209.                 tmp = tx;
  210.                 /* sum the least significant terms first */
  211.                 coef = &_default_gamma_coef[xgamma_dim];
  212.                 if (i == 0) {
  213.                         value = *--coef / tmp;
  214.                         while (coef > _default_gamma_coef)
  215.                                 value += *--coef / --tmp;
  216.                 }
  217.                 else {
  218.                         value = *--coef / BLAST_Powi(tmp, i + 1);
  219.                         while (coef > _default_gamma_coef)
  220.                                 value += *--coef / BLAST_Powi(--tmp, i + 1);
  221.                         tmp = BLAST_Factorial(i);
  222.                         value *= (i%2 == 0 ? tmp : -tmp);
  223.                 }
  224.                 y[i] = value;
  225.         }
  226.         ++y[0];
  227.         value = LogDerivative(order, y);
  228.         tmp = tx + 0.5;
  229.         switch (order) {
  230.         case 0:
  231.                 value += ((NCBIMATH_LNPI+NCBIMATH_LN2) / 2.)
  232.                                         + (xx + 0.5) * log(tmp) - tmp;
  233.                 break;
  234.         case 1:
  235.                 value += log(tmp) - xgamma_dim / tmp;
  236.                 break;
  237.         case 2:
  238.                 value += (tmp + xgamma_dim) / (tmp * tmp);
  239.                 break;
  240.         case 3:
  241.                 value -= (1. + 2.*xgamma_dim / tmp) / (tmp * tmp);
  242.                 break;
  243.         case 4:
  244.                 value += 2. * (1. + 3.*xgamma_dim / tmp) / (tmp * tmp * tmp);
  245.                 break;
  246.         default:
  247.                 tmp = BLAST_Factorial(order - 2) * BLAST_Powi(tmp, 1 - order)
  248.                                 * (1. + (order - 1) * xgamma_dim / tmp);
  249.                 if (order % 2 == 0)
  250.                         value += tmp;
  251.                 else
  252.                         value -= tmp;
  253.                 break;
  254.         }
  255.         return value;
  256. }
  259. static double PolyGamma(double x, Int4 order) /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
  260. /* x is and derivatives */
  261. /* order is order of the derivative */
  262. /* order = 0, 1, 2, ...  ln(gamma), digamma, trigamma, ... */
  263. /* CAUTION:  the value of order is one less than the suggested "di" and
  264. "tri" prefixes of digamma, trigamma, etc.  In other words, the value
  265. of order is truly the order of the derivative.  */
  266. {
  267. Int4 k;
  268. double value, tmp;
  269. double y[POLYGAMMA_ORDER_MAX+1], sx;
  271. if (order < 0 || order > POLYGAMMA_ORDER_MAX) {
  272. /*
  273. ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: unsupported derivative order");
  274. */
  275. /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  276. return HUGE_VAL;
  277. }
  279. if (x >= 1.)
  280. return general_lngamma(x, order);
  282. if (x < 0.) {
  283. value = general_lngamma(1. - x, order);
  284. value = ((order - 1) % 2 == 0 ? value : -value);
  285. if (order == 0) {
  286. sx = sin(NCBIMATH_PI * x);
  287. sx = ABS(sx);
  288. if ( (x < -0.1 && (ceil(x) == x || sx < 2.*DBL_EPSILON))
  289. || sx == 0.) {
  290. /*
  291. ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  292. */
  293. /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  294. return HUGE_VAL;
  295. }
  296. value += NCBIMATH_LNPI - log(sx);
  297. }
  298. else {
  299. y[0] = sin(x *= NCBIMATH_PI);
  300. tmp = 1.;
  301. for (k = 1; k <= order; k++) {
  302. tmp *= NCBIMATH_PI;
  303. y[k] = tmp * sin(x += (NCBIMATH_PI/2.));
  304. }
  305. value -= LogDerivative(order, y);
  306. }
  307. }
  308. else {
  309. value = general_lngamma(1. + x, order);
  310. if (order == 0) {
  311. if (x == 0.) {
  312. /*
  313. ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  314. */
  315. /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  316. return HUGE_VAL;
  317. }
  318. value -= log(x);
  319. }
  320. else {
  321. tmp = BLAST_Factorial(order - 1) * BLAST_Powi(x,  -order);
  322. value += (order % 2 == 0 ? tmp : - tmp);
  323. }
  324. }
  325. return value;
  326. }
  328. static double LnGamma(double x)               /* ln(ABS[gamma(x)]) - 10 dig
  329. its of accuracy */
  330. {
  331.         return PolyGamma(x, 0);
  332. }
  334. #define FACTORIAL_PRECOMPUTED   36
  336. extern double BLAST_Factorial(Int4 n)
  337. {
  338.         static double      precomputed[FACTORIAL_PRECOMPUTED]
  339.                 = { 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880., 3628800.};
  340.         static Int4     nlim = 10;
  341.         Int4   m;
  342.         double    x;
  344.         if (n >= 0) {
  345.                 if (n <= nlim)
  346.                         return precomputed[n];
  347.                 if (n < DIM(precomputed)) {
  348.                         for (x = precomputed[m = nlim]; m < n; ) {
  349.                                 ++m;
  350.                                 precomputed[m] = (x *= m);
  351.                         }
  352.                         nlim = m;
  353.                         return x;
  354.                 }
  355.                 return exp(LnGamma((double)(n+1)));
  356.         }
  357.         return 0.0; /* Undefined! */
  358. }
  360. /* LnGammaInt(n) -- return log(Gamma(n)) for integral n */
  361. extern double BLAST_LnGammaInt(Int4 n)
  362. {
  363. static double precomputed[FACTORIAL_PRECOMPUTED];
  364. static Int4 nlim = 1; /* first two entries are 0 */
  365. Int4 m;
  367. if (n >= 0) {
  368. if (n <= nlim)
  369. return precomputed[n];
  370. if (n < DIM(precomputed)) {
  371. for (m = nlim; m < n; ++m) {
  372. precomputed[m+1] = log(BLAST_Factorial(m));
  373. }
  374. return precomputed[nlim = m];
  375. }
  376. }
  377. return LnGamma((double)n);
  378. }
  381. /*
  382. Romberg numerical integrator
  384. Author:  Dr. John Spouge, NCBI
  385. Received:  July 17, 1992
  386. Reference:
  387. Francis Scheid (1968)
  388. Schaum's Outline Series
  389. Numerical Analysis, p. 115
  390. McGraw-Hill Book Company, New York
  391. */
  392. #define F(x)  ((*f)((x), fargs))
  393. #define ROMBERG_ITMAX 20
  395. extern double BLAST_RombergIntegrate(double (*f) (double,void*), void* fargs, double p, double q, double eps, Int4 epsit, Int4 itmin)
  397. {
  398. double romb[ROMBERG_ITMAX]; /* present list of Romberg values */
  399. double h; /* mesh-size */
  400. Int4 i, j, k, npts;
  401. long n; /* 4^(error order in romb[i]) */
  402. Int4 epsit_cnt = 0, epsck;
  403. double y;
  404. double x;
  405. double sum;
  407. /* itmin = min. no. of iterations to perform */
  408. itmin = MAX(1, itmin);
  409. itmin = MIN(itmin, ROMBERG_ITMAX-1);
  411. /* epsit = min. no. of consecutive iterations that must satisfy epsilon */
  412. epsit = MAX(epsit, 1); /* default = 1 */
  413. epsit = MIN(epsit, 3); /* if > 3, the problem needs more prior analysis */
  415. epsck = itmin - epsit;
  417. npts = 1;
  418. h = q - p;
  419. x = F(p);
  420. if (ABS(x) == HUGE_VAL)
  421. return x;
  422. y = F(q);
  423. if (ABS(y) == HUGE_VAL)
  424. return y;
  425. romb[0] = 0.5 * h * (x + y); /* trapezoidal rule */
  426. for (i = 1; i < ROMBERG_ITMAX; ++i, npts *= 2, h *= 0.5) {
  427. sum = 0.; /* sum of ordinates for */
  428. /* x = p+0.5*h, p+1.5*h, ..., q-0.5*h */
  429. for (k = 0, x = p+0.5*h; k < npts; k++, x += h) {
  430. y = F(x);
  431. if (ABS(y) == HUGE_VAL)
  432. return y;
  433. sum += y;
  434. }
  435. romb[i] = 0.5 * (romb[i-1] + h*sum); /* new trapezoidal estimate */
  437. /* Update Romberg array with new column */
  438. for (n = 4, j = i-1; j >= 0; n *= 4, --j)
  439. romb[j] = (n*romb[j+1] - romb[j]) / (n-1);
  441. if (i > epsck) {
  442. if (ABS(romb[1] - romb[0]) > eps * ABS(romb[0])) {
  443. epsit_cnt = 0;
  444. continue;
  445. }
  446. ++epsit_cnt;
  447. if (i >= itmin && epsit_cnt >= epsit)
  448. return romb[0];
  449. }
  450. }
  451. /*
  452. ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"RombergIntegrate: iterations > ROMBERG_ITMAX");
  453. */
  454. /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > ROMBERG_ITMAX"));**/
  455. return HUGE_VAL;
  456. }
  458. /*
  459. Gcd(a, b)
  461. Return the greatest common divisor of a and b.
  463. Adapted 8-15-90 by WRG from code by S. Altschul.
  464. */
  465. long BLAST_Gcd(long a, long b)
  466. {
  467. long c;
  469. b = ABS(b);
  470. if (b > a)
  471. c=a, a=b, b=c;
  473. while (b != 0) {
  474. c = a%b;
  475. a = b;
  476. b = c;
  477. }
  478. return a;
  479. }
  481. /* Round a floating point number to the nearest integer */
  482. long BLAST_Nint(double x) /* argument */
  483. {
  484. x += (x >= 0. ? 0.5 : -0.5);
  485. return (long)x;
  486. }
  488. /*
  489. integer power function
  491. Original submission by John Spouge, 6/25/90
  492. Added to shared library by WRG
  493. */
  494. extern double BLAST_Powi(double x, Int4 n) /* power */
  495. {
  496. double y;
  498. if (n == 0)
  499. return 1.;
  501. if (x == 0.) {
  502. if (n < 0) {
  503. /*
  504. ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"Powi: divide by 0");
  505. */
  506. /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  507. return HUGE_VAL;
  508. }
  509. return 0.;
  510. }
  512. if (n < 0) {
  513. x = 1./x;
  514. n = -n;
  515. }
  517. while (n > 1) {
  518. if (n&1) {
  519. y = x;
  520. goto Loop2;
  521. }
  522. n /= 2;
  523. x *= x;
  524. }
  525. return x;
  527. Loop2:
  528. n /= 2;
  529. x *= x;
  530. while (n > 1) {
  531. if (n&1)
  532. y *= x;
  533. n /= 2;
  534. x *= x;
  535. }
  536. return y * x;
  537. }
  539. extern double BLAST_LnFactorial (double x) {
  541.     if(x<=0.0)
  542.         return 0.0;
  543.     else
  544.         return LnGamma(x+1.0);
  545.         
  546. }
  548. /*
  549.  * ===========================================================================
  550.  *
  551.  * $Log: ncbi_math.c,v $
  552.  * Revision 1000.2  2004/06/01 18:08:20  gouriano
  553.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.8
  554.  *
  555.  * Revision 1.8  2004/05/19 14:52:03  camacho
  556.  * 1. Added doxygen tags to enable doxygen processing of algo/blast/core
  557.  * 2. Standardized copyright, CVS $Id string, $Log and rcsid formatting and i
  558.  *    location
  559.  * 3. Added use of @todo doxygen keyword
  560.  *
  561.  * Revision 1.7  2003/12/05 16:03:57  camacho
  562.  * Remove compiler warnings
  563.  *
  564.  * Revision 1.6  2003/09/26 20:39:32  dondosha
  565.  * Rearranged code so it compiles
  566.  *
  567.  * Revision 1.5  2003/09/26 19:01:59  madden
  568.  * Prefix ncbimath functions with BLAST_
  569.  *
  570.  * Revision 1.4  2003/09/10 21:36:29  dondosha
  571.  * Removed Nlm_ prefix from math functions definitions
  572.  *
  573.  * Revision 1.3  2003/08/25 22:32:51  dondosha
  574.  * Added #ifndef for definition of DBL_EPSILON
  575.  *
  576.  * Revision 1.2  2003/08/11 15:02:00  dondosha
  577.  * Added algo/blast/core to all #included headers
  578.  *
  579.  * Revision 1.1  2003/08/02 16:31:48  camacho
  580.  * Moved ncbimath.c -> ncbi_math.c
  581.  *
  582.  * Revision 1.1  2003/08/01 21:03:46  madden
  583.  * Cleaned up version of file for C++ toolkit
  584.  *
  585.  * ===========================================================================
  586.  */