home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / corelib / ncbimath.c < prev    next >
Encoding:
Text File  |  1995-12-17  |  17.4 KB  |  751 lines  |  [TEXT/R*ch]
  1. /*   ncbimath.c
  2. * ===========================================================================
  3. *
  4. *                            PUBLIC DOMAIN NOTICE
  5. *               National Center for Biotechnology Information
  6. *
  7. *  This software/database is a "United States Government Work" under the
  8. *  terms of the United States Copyright Act.  It was written as part of
  9. *  the author's official duties as a United States Government employee and
  10. *  thus cannot be copyrighted.  This software/database is freely available
  11. *  to the public for use. The National Library of Medicine and the U.S.
  12. *  Government have not placed any restriction on its use or reproduction.
  13. *
  14. *  Although all reasonable efforts have been taken to ensure the accuracy
  15. *  and reliability of the software and data, the NLM and the U.S.
  16. *  Government do not and cannot warrant the performance or results that
  17. *  may be obtained by using this software or data. The NLM and the U.S.
  18. *  Government disclaim all warranties, express or implied, including
  19. *  warranties of performance, merchantability or fitness for any particular
  20. *  purpose.
  21. *
  22. *  Please cite the author in any work or product based on this material.
  23. *
  24. * ===========================================================================
  25. *
  26. * File Name:  ncbimath.c
  27. *
  28. * Author:  Gish, Kans, Ostell, Schuler
  29. *
  30. * Version Creation Date:   10/23/91
  31. *
  32. * $Revision: 2.7 $
  33. *
  34. * File Description:
  35. *       portable math functions
  36. *
  37. * Modifications:
  38. * --------------------------------------------------------------------------
  39. * Date     Name        Description of modification
  40. * -------  ----------  -----------------------------------------------------
  41. * 04-15-93 Schuler     Changed _cdecl to LIBCALL
  42. * 12-22-93 Schuler     Converted ERRPOST((...)) to ErrPostEx(...)
  43. *
  44. * ==========================================================================
  45. */
  46.  
  47. #undef  THIS_MODULE
  48. #define THIS_MODULE g_corelib
  49. #define THI_FILE _this_file
  50.  
  51. #include <ncbi.h>
  52. #include <ncbiwin.h>
  53.  
  54. extern char * g_corelib;
  55. static char * _this_file = __FILE__;
  56.  
  57. /*
  58.     Nlm_Expm1(x)
  59.     Return values accurate to approx. 16 digits for the quantity exp(x)-1
  60.     for all x.
  61. */
  62. double LIBCALL
  63. Nlm_Expm1(x)
  64.     register double    x;
  65. {
  66.     register double    absx;
  67.  
  68.     if ((absx = ABS(x)) > .33)
  69.         return exp(x) - 1.;
  70.  
  71.     if (absx < 1.e-16)
  72.         return x;
  73.  
  74.     return x * (1. + x *
  75.                     (0.5 + x * (1./6. + x *
  76.                         (1./24. + x * (1./120. + x *
  77.                             (1./720. + x * (1./5040. + x *
  78.                                 (1./40320. + x * (1./362880. + x *
  79.                                     (1./3628800. + x * (1./39916800. + x *
  80.                                         (1./479001600. + x/6227020800.)
  81.                                     ))
  82.                                 ))
  83.                             ))
  84.                         ))
  85.                     )));
  86. }
  87.  
  88. /*
  89.     Nlm_Log1p(x)
  90.     Return accurate values for the quantity log(x+1) for all x > -1.
  91. */
  92. double LIBCALL
  93. Nlm_Log1p(x)
  94.     register double    x;
  95. {
  96.     register int    i;
  97.     register double    sum, y;
  98.  
  99.     if (ABS(x) >= 0.2)
  100.         return log(x+1.);
  101.  
  102.     for (i=0, sum=0., y=x; ; ) {
  103.         sum += y/++i;
  104.         if (y < DBL_EPSILON)
  105.             break;
  106.         y *= x;
  107.         sum -= y/++i;
  108.         if (y < DBL_EPSILON)
  109.             break;
  110.         y *= x;
  111.     }
  112.     return sum;
  113. }
  114.  
  115.  
  116. /*
  117. ** Special thanks to Dr. John ``Gammas Galore'' Spouge for deriving the
  118. ** method for calculating the gamma coefficients and their use.
  119. ** (See the #ifdef-ed program included at the end of this file).
  120. **/
  121.  
  122. /*
  123.     For discussion of the Gamma function, see "Numerical Recipes in C",
  124.     Press et al. (1988) pages 167-169.
  125. */
  126.  
  127. static double NEAR general_lngamma PROTO((double x,int n));
  128.  
  129. static double _default_gamma_coef [] = {
  130.      4.694580336184385e+04,
  131.     -1.560605207784446e+05,
  132.      2.065049568014106e+05,
  133.     -1.388934775095388e+05,
  134.      5.031796415085709e+04,
  135.     -9.601592329182778e+03,
  136.      8.785855930895250e+02,
  137.     -3.155153906098611e+01,
  138.      2.908143421162229e-01,
  139.     -2.319827630494973e-04,
  140.      1.251639670050933e-10
  141.      };
  142. static double    PNTR gamma_coef = _default_gamma_coef;
  143. static unsigned    gamma_dim = DIM(_default_gamma_coef);
  144. static double    xgamma_dim = DIM(_default_gamma_coef);
  145.  
  146. void LIBCALL
  147. Nlm_GammaCoeffSet(double PNTR cof, unsigned dim) /* changes gamma coeffs */
  148. {
  149.     if (dim < 3 || dim > 100) /* sanity check */
  150.         return;
  151.     gamma_coef = cof;
  152.     xgamma_dim = gamma_dim = dim;
  153. }
  154.  
  155.  
  156. static double NEAR
  157. general_lngamma(x, order)    /* nth derivative of ln[gamma(x)] */
  158.     double    x;        /* 10-digit accuracy achieved only for 1 <= x */
  159.     int        order;    /* order of the derivative, 0..POLYGAMMA_ORDER_MAX */
  160. {
  161.     int        i;
  162.     double    xx, tx;
  163.     double    y[POLYGAMMA_ORDER_MAX+1];
  164.     register double    tmp, value, PNTR coef;
  165.  
  166.     xx = x - 1.;  /* normalize from gamma(x + 1) to xx! */
  167.     tx = xx + xgamma_dim;
  168.     for (i = 0; i <= order; ++i) {
  169.         tmp = tx;
  170.         /* sum the least significant terms first */
  171.         coef = &gamma_coef[gamma_dim];
  172.         if (i == 0) {
  173.             value = *--coef / tmp;
  174.             while (coef > gamma_coef)
  175.                 value += *--coef / --tmp;
  176.         }
  177.         else {
  178.             value = *--coef / Nlm_Powi(tmp, i + 1);
  179.             while (coef > gamma_coef)
  180.                 value += *--coef / Nlm_Powi(--tmp, i + 1);
  181.             tmp = Nlm_Factorial(i);
  182.             value *= (i%2 == 0 ? tmp : -tmp);
  183.         }
  184.         y[i] = value;
  185.     }
  186.     ++y[0];
  187.     value = Nlm_LogDerivative(order, y);
  188.     tmp = tx + 0.5;
  189.     switch (order) {
  190.     case 0:
  191.         value += ((NCBIMATH_LNPI+NCBIMATH_LN2) / 2.)
  192.                     + (xx + 0.5) * log(tmp) - tmp;
  193.         break;
  194.     case 1:
  195.         value += log(tmp) - xgamma_dim / tmp;
  196.         break;
  197.     case 2:
  198.         value += (tmp + xgamma_dim) / (tmp * tmp);
  199.         break;
  200.     case 3:
  201.         value -= (1. + 2.*xgamma_dim / tmp) / (tmp * tmp);
  202.         break;
  203.     case 4:
  204.         value += 2. * (1. + 3.*xgamma_dim / tmp) / (tmp * tmp * tmp);
  205.         break;
  206.     default:
  207.         tmp = Nlm_Factorial(order - 2) * Nlm_Powi(tmp, 1 - order)
  208.                 * (1. + (order - 1) * xgamma_dim / tmp);
  209.         if (order % 2 == 0)
  210.             value += tmp;
  211.         else
  212.             value -= tmp;
  213.         break;
  214.     }
  215.     return value;
  216. }
  217.  
  218.  
  219. double LIBCALL
  220. Nlm_PolyGamma(x, order) /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
  221.     double    x;      /* and derivatives */
  222.     int        order;    /* order of the derivative */
  223. /* order = 0, 1, 2, ...  ln(gamma), digamma, trigamma, ... */
  224. /* CAUTION:  the value of order is one less than the suggested "di" and
  225. "tri" prefixes of digamma, trigamma, etc.  In other words, the value
  226. of order is truly the order of the derivative.  */
  227. {
  228.     int        k;
  229.     register double    value, tmp;
  230.     double    y[POLYGAMMA_ORDER_MAX+1], sx;
  231.  
  232.     if (order < 0 || order > POLYGAMMA_ORDER_MAX) {
  233.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: unsupported derivative order");
  234.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  235.         return HUGE_VAL;
  236.     }
  237.  
  238.     if (x >= 1.)
  239.         return general_lngamma(x, order);
  240.  
  241.     if (x < 0.) {
  242.         value = general_lngamma(1. - x, order);
  243.         value = ((order - 1) % 2 == 0 ? value : -value);
  244.         if (order == 0) {
  245.             sx = sin(NCBIMATH_PI * x);
  246.             sx = ABS(sx);
  247.             if ( (x < -0.1 && (ceil(x) == x || sx < 2.*DBL_EPSILON))
  248.                     || sx == 0.) {
  249.                 ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  250.                 /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  251.                 return HUGE_VAL;
  252.             }
  253.             value += NCBIMATH_LNPI - log(sx);
  254.         }
  255.         else {
  256.             y[0] = sin(x *= NCBIMATH_PI);
  257.             tmp = 1.;
  258.             for (k = 1; k <= order; k++) {
  259.                 tmp *= NCBIMATH_PI;
  260.                 y[k] = tmp * sin(x += (NCBIMATH_PI/2.));
  261.             }
  262.             value -= Nlm_LogDerivative(order, y);
  263.         }
  264.     }
  265.     else {
  266.         value = general_lngamma(1. + x, order);
  267.         if (order == 0) {
  268.             if (x == 0.) {
  269.                 ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"PolyGamma: log(0)");
  270.                 /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(0)"));**/
  271.                 return HUGE_VAL;
  272.             }
  273.             value -= log(x);
  274.         }
  275.         else {
  276.             tmp = Nlm_Factorial(order - 1) * Nlm_Powi(x,  -order);
  277.             value += (order % 2 == 0 ? tmp : - tmp);
  278.         }
  279.     }
  280.     return value;
  281. }
  282.  
  283. double LIBCALL
  284. Nlm_LogDerivative(order, u) /* nth derivative of ln(u) */
  285.     int        order;        /* order of the derivative */
  286.     double    PNTR u;        /* values of u, u', u", etc. */
  287. {
  288.     int        i;
  289.     double    y[LOGDERIV_ORDER_MAX+1];
  290.     register double    value, tmp;
  291.  
  292.     if (order < 0 || order > LOGDERIV_ORDER_MAX) {
  293. InvalidOrder:
  294.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: unsupported derivative order");
  295.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "unsupported derivative order"));**/
  296.         return HUGE_VAL;
  297.     }
  298.  
  299.     if (order > 0 && u[0] == 0.) {
  300.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: divide by 0");
  301.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  302.         return HUGE_VAL;
  303.     }
  304.     for (i = 1; i <= order; i++)
  305.         y[i] = u[i] / u[0];
  306.  
  307.     switch (order) {
  308.     case 0:
  309.         if (u[0] > 0.)
  310.             value = log(u[0]);
  311.         else {
  312.             ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"LogDerivative: log(x<=0)");
  313.             /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "log(x<=0)"));**/
  314.             return HUGE_VAL;
  315.         }
  316.         break;
  317.     case 1:
  318.         value = y[1];
  319.         break;
  320.     case 2:
  321.         value = y[2] - y[1] * y[1];
  322.         break;
  323.     case 3:
  324.         value = y[3] - 3. * y[2] * y[1] + 2. * y[1] * y[1] * y[1];
  325.         break;
  326.     case 4:
  327.         value = y[4] - 4. * y[3] * y[1] - 3. * y[2] * y[2]
  328.                           + 12. * y[2] * (tmp = y[1] * y[1]);
  329.         value -= 6. * tmp * tmp;
  330.         break;
  331.     default:
  332.         goto InvalidOrder;
  333.     }
  334.     return value;
  335. }
  336.  
  337.  
  338.  
  339. double LIBCALL
  340. Nlm_Gamma(x)        /* ABS[gamma(x)] - 10 digits of accuracy */
  341.     double     x;
  342. {
  343.     return exp(Nlm_PolyGamma(x, 0));
  344. }
  345.  
  346. double LIBCALL
  347. Nlm_LnGamma(x)        /* ln(ABS[gamma(x)]) - 10 digits of accuracy */
  348.     double     x;
  349. {
  350.     return Nlm_PolyGamma(x, 0);
  351. }
  352.  
  353. double LIBCALL
  354. Nlm_DiGamma(x)    /* digamma, 1st order derivative of log(gamma) */
  355.     double    x;
  356. {
  357.     return Nlm_PolyGamma(x, 1);
  358. }
  359.  
  360. double LIBCALL
  361. Nlm_TriGamma(x)    /* trigamma, 2nd order derivative of log(gamma) */
  362.     double    x;
  363. {
  364.     return Nlm_PolyGamma(x, 2);
  365. }
  366.  
  367.  
  368. #ifdef foo
  369. /* A program to calculate the gamma coefficients based on a method
  370.    by John Spouge.
  371.    Cut this program out, save it in a separate file, and compile.
  372.    Be sure to link with a math library.
  373. */
  374. /*
  375.     a[n] = ((gamma+0.5-n)^(n-0.5)) * exp(gamma+0.5-n) *
  376.         ((-1)^(n-1) / (n-1)!) * (1/sqrt(2*Pi))
  377. */
  378. #include <stdio.h>
  379. #include <math.h>
  380.  
  381. main(ac, av)
  382.     int ac;
  383.     char **av;
  384.  
  385. {
  386.     int    i, j, cnt;
  387.     double    a, x, y, z, ifact = 1.;
  388.  
  389.     if (ac != 2 || sscanf(av[1], "%d", &cnt) != 1)
  390.         exit(1);
  391.  
  392.     for (i=0; i<cnt; ++i) {
  393.         x = cnt - (i + 0.5);
  394.         y = log(x) * (i + 0.5) + x;
  395.         y = exp(y);
  396.         if (i > 1)
  397.             ifact *= i;
  398.         y /= ifact;
  399.         if (i%2 == 1)
  400.             y = -y;
  401.         printf("\t\t\t%.17lg", y);
  402.         if (i < cnt-1)
  403.             putchar(',');
  404.         putchar('\n');
  405.     }
  406. }
  407. #endif /* foo */
  408.  
  409.  
  410. #define FACTORIAL_PRECOMPUTED    36
  411.  
  412. double LIBCALL
  413. Nlm_Factorial(int n)
  414.  
  415. {
  416.     static double    precomputed[FACTORIAL_PRECOMPUTED]
  417.                 = { 1., 1., 2., 6., 24., 120., 720.};
  418.     static int    nlim = 6;
  419.  
  420.     if (n < 0)
  421.         return 0.0; /* Undefined! */
  422.  
  423.     if (n >= DIM(precomputed))
  424.         return exp(Nlm_LnGamma((double)(n+1)));
  425.  
  426.     while (nlim < n) {
  427.         ++nlim;
  428.         precomputed[nlim] = precomputed[nlim-1] * nlim;
  429.     }
  430.     return precomputed[n];
  431. }
  432.  
  433. /* Nlm_LnGammaInt(n) -- return log(Gamma(n)) for integral n */
  434. double LIBCALL
  435. Nlm_LnGammaInt(int n)
  436. {
  437.     static double    precomputed[FACTORIAL_PRECOMPUTED];
  438.     static int    nlim = 1; /* first two entries are 0 */
  439.  
  440.     if (n >= 0 && n < DIM(precomputed)) {
  441.         while (nlim < n) {
  442.             ++nlim;
  443.             precomputed[nlim] = log(Nlm_Factorial(nlim-1));
  444.         }
  445.         return precomputed[n];
  446.     }
  447.     return Nlm_LnGamma((double)n);
  448. }
  449.  
  450.  
  451.  
  452. /*
  453.     Combined Newton-Raphson and Bisection root-finder
  454.  
  455.     Original Function Name:  Inv_Xnrbis()
  456.     Author:  Dr. John Spouge
  457.     Location:  NCBI
  458.     Received:  July 16, 1991
  459. */
  460. #define F(x)  ((*f)(x)-y)
  461. #define DF(x) ((*df)(x))
  462. #define NRBIS_ITMAX 100
  463.  
  464. double LIBCALL
  465. Nlm_NRBis(double y, double (LIBCALL *f )PROTO ((double )), double (LIBCALL *df )PROTO ((double )), double p, double x, double q, double tol)        /* tolerance */
  466.  
  467. {
  468.     double    temp;    /* for swapping end-points if necessary */
  469.     double    dx;        /* present interval length */
  470.     double    dxold;    /* old interval length */
  471.     double    fx;        /* f(x)-y */
  472.     double    dfx;    /* Df(x) */
  473.     int    j;                       /* loop index */
  474.     double    fp, fq;
  475.  
  476. /* Preliminary checks for improper bracketing and end-point root. */
  477.     if ((fp = F(p)) == 0.)
  478.         return p;
  479.     if ((fq = F(q)) == 0.)
  480.         return q;
  481.     if ((fp > 0. && fq > 0.) || (fp < 0. && fq < 0.)) {
  482.         ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_INVAL,"NRBis: root not bracketed");
  483.         /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_INVAL, "root not bracketed"));**/
  484.         return HUGE_VAL;
  485.     }
  486. /* Swaps end-points if necessary to make F(p)<0.<F(q) */
  487.     if (fp > 0.) {
  488.         temp = p;
  489.         p = q;
  490.         q = temp;
  491.     }
  492. /* Set up the Bisection & Newton-Raphson iteration. */
  493.     if ((x-p) * (x-q) > 0.)
  494.         x = 0.5*(p+q);
  495.     dxold = dx = p-q;
  496.     for (j = 1; j <= NRBIS_ITMAX; ++j) {
  497.         fx = F(x);
  498.         if (fx == 0.) /* Check for termination. */
  499.             return x;
  500.         if (fx < 0.)
  501.             p = x;
  502.         else
  503.             q = x;
  504.         dfx = DF(x);
  505. /* Check: root out of bounds or bisection faster than Newton-Raphson? */
  506.         if ((dfx*(x-p)-fx)*(dfx*(x-q)-fx) >= 0. || 2.*ABS(fx) > ABS(dfx*dx)) {
  507.             dx = dxold;               /* Bisect */
  508.             dxold = 0.5*(p-q);
  509.             x = 0.5*(p+q);
  510.             if (ABS(dxold) <= tol)
  511.                 return x;
  512.         } else {
  513.             dx = dxold;               /* Newton-Raphson */
  514.             dxold = fx/dfx;
  515.             x -= dxold;
  516.             if (ABS(dxold) < tol && F(x-SIGN(dxold)*tol)*fx < 0.)
  517.                 return x;
  518.         }
  519.     }
  520.     
  521.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"NRBis: iterations > NRBIS_ITMAX");
  522.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > NRBIS_ITMAX"));**/
  523.     return HUGE_VAL;
  524. }
  525. #undef F /* clean-up */
  526. #undef DF /* clean-up */
  527.  
  528.  
  529.  
  530. /*
  531.     Romberg numerical integrator
  532.  
  533.     Author:  Dr. John Spouge, NCBI
  534.     Received:  July 17, 1992
  535.     Reference:
  536.         Francis Scheid (1968)
  537.         Schaum's Outline Series
  538.         Numerical Analysis, p. 115
  539.         McGraw-Hill Book Company, New York
  540. */
  541. #define F(x)  ((*f)((x), fargs))
  542. #define ROMBERG_ITMAX 20
  543.  
  544. double LIBCALL
  545. Nlm_RombergIntegrate(double (LIBCALL *f) (double,Nlm_VoidPtr), Nlm_VoidPtr fargs, double p, double q, double eps, int epsit, int itmin)
  546.  
  547. {
  548.     double    romb[ROMBERG_ITMAX];    /* present list of Romberg values */
  549.     double    h;    /* mesh-size */
  550.     int        i, j, k, npts;
  551.     long    n;    /* 4^(error order in romb[i]) */
  552.     int        epsit_cnt = 0, epsck;
  553.     register double    y;
  554.     register double    x;
  555.     register double    sum;
  556.  
  557.     /* itmin = min. no. of iterations to perform */
  558.     itmin = MAX(1, itmin);
  559.     itmin = MIN(itmin, ROMBERG_ITMAX-1);
  560.  
  561.     /* epsit = min. no. of consecutive iterations that must satisfy epsilon */
  562.     epsit = MAX(epsit, 1); /* default = 1 */
  563.     epsit = MIN(epsit, 3); /* if > 3, the problem needs more prior analysis */
  564.  
  565.     epsck = itmin - epsit;
  566.  
  567.     npts = 1;
  568.     h = q - p;
  569.     x = F(p);
  570.     if (ABS(x) == HUGE_VAL)
  571.         return x;
  572.     y = F(q);
  573.     if (ABS(y) == HUGE_VAL)
  574.         return y;
  575.     romb[0] = 0.5 * h * (x + y);    /* trapezoidal rule */
  576.     for (i = 1; i < ROMBERG_ITMAX; ++i, npts *= 2, h *= 0.5) {
  577.         sum = 0.;    /* sum of ordinates for */
  578.         /* x = p+0.5*h, p+1.5*h, ..., q-0.5*h */
  579.         for (k = 0, x = p+0.5*h; k < npts; k++, x += h) {
  580.             y = F(x);
  581.             if (ABS(y) == HUGE_VAL)
  582.                 return y;
  583.             sum += y;
  584.         }
  585.         romb[i] = 0.5 * (romb[i-1] + h*sum); /* new trapezoidal estimate */
  586.  
  587.         /* Update Romberg array with new column */
  588.         for (n = 4, j = i-1; j >= 0; n *= 4, --j)
  589.             romb[j] = (n*romb[j+1] - romb[j]) / (n-1);
  590.  
  591.         if (i > epsck) {
  592.             if (ABS(romb[1] - romb[0]) > eps * ABS(romb[0])) {
  593.                 epsit_cnt = 0;
  594.                 continue;
  595.             }
  596.             ++epsit_cnt;
  597.             if (i >= itmin && epsit_cnt >= epsit)
  598.                 return romb[0];
  599.         }
  600.     }
  601.     ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"RombergIntegrate: iterations > ROMBERG_ITMAX");
  602.     /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > ROMBERG_ITMAX"));**/
  603.     return HUGE_VAL;
  604. }
  605.  
  606. /*
  607.     Nlm_Gcd(a, b)
  608.  
  609.     Return the greatest common divisor of a and b.
  610.  
  611.     Adapted 8-15-90 by WRG from code by S. Altschul.
  612. */
  613. long LIBCALL
  614. Nlm_Gcd(register long a, register long b)
  615. {
  616.     register long    c;
  617.  
  618.     b = ABS(b);
  619.     if (b > a)
  620.         c=a, a=b, b=c;
  621.  
  622.     while (b != 0) {
  623.         c = a%b;
  624.         a = b;
  625.         b = c;
  626.     }
  627.     return a;
  628. }
  629.  
  630. /* Round a floating point number to the nearest integer */
  631. long LIBCALL
  632. Nlm_Nint(register double x)    /* argument */
  633. {
  634.     x += (x >= 0. ? 0.5 : -0.5);
  635.     return (long)x;
  636. }
  637.  
  638. /*
  639. integer power function
  640.  
  641. Original submission by John Spouge, 6/25/90
  642. Added to shared library by WRG
  643. */
  644. double LIBCALL
  645. Nlm_Powi(register double x, register int n)    /* power */
  646. {
  647.     register double    y;
  648.  
  649.     if (n == 0)
  650.         return 1.;
  651.  
  652.     if (x == 0.) {
  653.         if (n < 0) {
  654.             ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_DOMAIN,"Powi: divide by 0");
  655.             /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_DOMAIN, "divide by 0"));**/
  656.             return HUGE_VAL;
  657.         }
  658.         return 0.;
  659.     }
  660.  
  661.     if (n < 0) {
  662.         x = 1./x;
  663.         n = -n;
  664.     }
  665.  
  666.     while (n > 1) {
  667.         if (n&1) {
  668.             y = x;
  669.             goto Loop2;
  670.         }
  671.         n /= 2;
  672.         x *= x;
  673.     }
  674.     return x;
  675.  
  676. Loop2:
  677.     n /= 2;
  678.     x *= x;
  679.     while (n > 1) {
  680.         if (n&1)
  681.             y *= x;
  682.         n /= 2;
  683.         x *= x;
  684.     }
  685.     return y * x;
  686. }
  687.  
  688. /*
  689.     Additive random number generator
  690.  
  691.     Modelled after "Algorithm A" in
  692.     Knuth, D. E. (1981). The art of computer programming, volume 2, page 27.
  693.  
  694.     7/26/90 WRG
  695. */
  696.  
  697. static long    state[33] = {
  698.     (long)0xd53f1852,  (long)0xdfc78b83,  (long)0x4f256096,  (long)0xe643df7,
  699.     (long)0x82c359bf,  (long)0xc7794dfa,  (long)0xd5e9ffaa,  (long)0x2c8cb64a,
  700.     (long)0x2f07b334,  (long)0xad5a7eb5,  (long)0x96dc0cde,  (long)0x6fc24589,
  701.     (long)0xa5853646,  (long)0xe71576e2,  (long)0xdae30df,  (long)0xb09ce711,
  702.     (long)0x5e56ef87,  (long)0x4b4b0082,  (long)0x6f4f340e,  (long)0xc5bb17e8,
  703.     (long)0xd788d765,  (long)0x67498087,  (long)0x9d7aba26,  (long)0x261351d4,
  704.     (long)0x411ee7ea,  (long)0x393a263,  (long)0x2c5a5835,  (long)0xc115fcd8,
  705.     (long)0x25e9132c,  (long)0xd0c6e906,  (long)0xc2bc5b2d,  (long)0x6c065c98,
  706.     (long)0x6e37bd55 };
  707.  
  708. #define r_off    12
  709.  
  710. static long    *rJ = &state[r_off],
  711.             *rK = &state[DIM(state)-1];
  712.  
  713. void LIBCALL
  714. Nlm_RandomSeed(long x)
  715. {
  716.     register int    i;
  717.  
  718.     state[0] = x;
  719.     /* linear congruential initializer */
  720.     for (i=1; i<DIM(state); ++i) {
  721.         state[i] = 1103515245*state[i-1] + 12345;
  722.     }
  723.  
  724.     rJ = &state[r_off];
  725.     rK = &state[DIM(state)-1];
  726.  
  727.     for (i=0; i<10*DIM(state); ++i)
  728.         (void) Nlm_RandomNum();
  729. }
  730.  
  731.  
  732. /*
  733.     Nlm_RandomNum --  return value in the range 0 <= x <= 2**31 - 1
  734. */
  735. long LIBCALL
  736. Nlm_RandomNum(void)
  737. {
  738.     register long    r;
  739.  
  740.     r = *rK;
  741.     r += *rJ--;
  742.     *rK-- = r;
  743.  
  744.     if (rK < state)
  745.         rK = &state[DIM(state)-1];
  746.     else
  747.         if (rJ < state)
  748.             rJ = &state[DIM(state)-1];
  749.     return (r>>1)&0x7fffffff; /* discard the least-random bit */
  750. }
  751.