/* ncbimath.c * =========================================================================== * * 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. * * =========================================================================== * * File Name: ncbimath.c * * Author: Gish, Kans, Ostell, Schuler * * Version Creation Date: 10/23/91 * * $Revision: 5.1 $ * * File Description: * portable math functions * * Modifications: * -------------------------------------------------------------------------- * Date Name Description of modification * ------- ---------- ----------------------------------------------------- * 04-15-93 Schuler Changed _cdecl to LIBCALL * 12-22-93 Schuler Converted ERRPOST((...)) to ErrPostEx(...) * * $Log: ncbimath.c,v $ * Revision 5.1 1996/06/20 14:08:00 madden * Changed int to Int4, double to Nlm_FloatHi * * Revision 5.0 1996/05/28 13:18:57 ostell * Set to revision 5.0 * * Revision 4.1 1996/03/06 19:47:15 epstein * fix problem observed by Epstein & fixed by Spouge in log calculation * * Revision 4.0 1995/07/26 13:46:50 ostell * force revision to 4.0 * * Revision 2.11 1995/05/15 18:45:58 ostell * added Log line * * * * ========================================================================== */ #define THIS_MODULE g_corelib #define THIS_FILE _this_file #include extern char * g_corelib; static char * _this_file = __FILE__; /* Nlm_Expm1(x) Return values accurate to approx. 16 digits for the quantity exp(x)-1 for all x. */ Nlm_FloatHi LIBCALL Nlm_Expm1(register Nlm_FloatHi x) { register Nlm_FloatHi 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.) )) )) )) )) ))); } /* Nlm_Log1p(x) Return accurate values for the quantity log(x+1) for all x > -1. */ Nlm_FloatHi LIBCALL Nlm_Log1p(register Nlm_FloatHi x) { register Int4 i; register Nlm_FloatHi sum, y; if (ABS(x) >= 0.2) return log(x+1.); for (i=0, sum=0., y=x; ; ) { sum += y/++i; if (ABS(y) < DBL_EPSILON) break; y *= x; sum -= y/++i; if (y < DBL_EPSILON) break; y *= x; } return sum; } /* ** Special thanks to Dr. John ``Gammas Galore'' Spouge for deriving the ** method for calculating the gamma coefficients and their use. ** (See the #ifdef-ed program included at the end of this file). **/ /* For discussion of the Gamma function, see "Numerical Recipes in C", Press et al. (1988) pages 167-169. */ static Nlm_FloatHi NEAR general_lngamma PROTO((Nlm_FloatHi x,Int4 n)); static Nlm_FloatHi _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 Nlm_FloatHi PNTR gamma_coef = _default_gamma_coef; static unsigned gamma_dim = DIM(_default_gamma_coef); static Nlm_FloatHi xgamma_dim = DIM(_default_gamma_coef); void LIBCALL Nlm_GammaCoeffSet(Nlm_FloatHi PNTR cof, unsigned dim) /* changes gamma coeffs */ { if (dim < 3 || dim > 100) /* sanity check */ return; gamma_coef = cof; xgamma_dim = gamma_dim = dim; } static Nlm_FloatHi NEAR general_lngamma(Nlm_FloatHi 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; Nlm_FloatHi xx, tx; Nlm_FloatHi y[POLYGAMMA_ORDER_MAX+1]; register Nlm_FloatHi tmp, value, PNTR 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 = &gamma_coef[gamma_dim]; if (i == 0) { value = *--coef / tmp; while (coef > gamma_coef) value += *--coef / --tmp; } else { value = *--coef / Nlm_Powi(tmp, i + 1); while (coef > gamma_coef) value += *--coef / Nlm_Powi(--tmp, i + 1); tmp = Nlm_Factorial(i); value *= (i%2 == 0 ? tmp : -tmp); } y[i] = value; } ++y[0]; value = Nlm_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 = Nlm_Factorial(order - 2) * Nlm_Powi(tmp, 1 - order) * (1. + (order - 1) * xgamma_dim / tmp); if (order % 2 == 0) value += tmp; else value -= tmp; break; } return value; } Nlm_FloatHi LIBCALL Nlm_PolyGamma(Nlm_FloatHi 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; register Nlm_FloatHi value, tmp; Nlm_FloatHi 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; } 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 -= Nlm_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 = Nlm_Factorial(order - 1) * Nlm_Powi(x, -order); value += (order % 2 == 0 ? tmp : - tmp); } } return value; } Nlm_FloatHi LIBCALL Nlm_LogDerivative(Int4 order, Nlm_FloatHi PNTR u) /* nth derivative of ln(u) */ /* order is order of the derivative */ /* u is values of u, u', u", etc. */ { Int4 i; Nlm_FloatHi y[LOGDERIV_ORDER_MAX+1]; register Nlm_FloatHi 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; } Nlm_FloatHi LIBCALL Nlm_Gamma(Nlm_FloatHi x) /* ABS[gamma(x)] - 10 digits of accuracy */ { return exp(Nlm_PolyGamma(x, 0)); } Nlm_FloatHi LIBCALL Nlm_LnGamma(Nlm_FloatHi x) /* ln(ABS[gamma(x)]) - 10 digits of accuracy */ { return Nlm_PolyGamma(x, 0); } Nlm_FloatHi LIBCALL Nlm_DiGamma(Nlm_FloatHi x) /* digamma, 1st order derivative of log(gamma) */ { return Nlm_PolyGamma(x, 1); } Nlm_FloatHi LIBCALL Nlm_TriGamma(Nlm_FloatHi x) /* trigamma, 2nd order derivative of log(gamma) */ { return Nlm_PolyGamma(x, 2); } #ifdef foo /* A program to calculate the gamma coefficients based on a method by John Spouge. Cut this program out, save it in a separate file, and compile. Be sure to link with a math library. */ /* a[n] = ((gamma+0.5-n)^(n-0.5)) * exp(gamma+0.5-n) * ((-1)^(n-1) / (n-1)!) * (1/sqrt(2*Pi)) */ #include #include main(ac, av) int ac; char **av; { int i, j, cnt; double a, x, y, z, ifact = 1.; if (ac != 2 || sscanf(av[1], "%d", &cnt) != 1) exit(1); for (i=0; i 1) ifact *= i; y /= ifact; if (i%2 == 1) y = -y; printf("\t\t\t%.17lg", y); if (i < cnt-1) putchar(','); putchar('\n'); } } #endif /* foo */ #define FACTORIAL_PRECOMPUTED 36 Nlm_FloatHi LIBCALL Nlm_Factorial(Int4 n) { static Nlm_FloatHi precomputed[FACTORIAL_PRECOMPUTED] = { 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880., 3628800.}; static Int4 nlim = 10; register Int4 m; register Nlm_FloatHi 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(Nlm_LnGamma((Nlm_FloatHi)(n+1))); } return 0.0; /* Undefined! */ } /* Nlm_LnGammaInt(n) -- return log(Gamma(n)) for integral n */ Nlm_FloatHi LIBCALL Nlm_LnGammaInt(Int4 n) { static Nlm_FloatHi precomputed[FACTORIAL_PRECOMPUTED]; static Int4 nlim = 1; /* first two entries are 0 */ register 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(Nlm_Factorial(m)); } return precomputed[nlim = m]; } } return Nlm_LnGamma((Nlm_FloatHi)n); } /* Combined Newton-Raphson and Bisection root-finder Original Function Name: Inv_Xnrbis() Author: Dr. John Spouge Location: NCBI Received: July 16, 1991 */ #define F(x) ((*f)(x)-y) #define DF(x) ((*df)(x)) #define NRBIS_ITMAX 100 Nlm_FloatHi LIBCALL Nlm_NRBis(Nlm_FloatHi y, Nlm_FloatHi (LIBCALL *f )PROTO ((Nlm_FloatHi )), Nlm_FloatHi (LIBCALL *df )PROTO ((Nlm_FloatHi )), Nlm_FloatHi p, Nlm_FloatHi x, Nlm_FloatHi q, Nlm_FloatHi tol) /* tolerance */ { Nlm_FloatHi temp; /* for swapping end-points if necessary */ Nlm_FloatHi dx; /* present interval length */ Nlm_FloatHi dxold; /* old interval length */ Nlm_FloatHi fx; /* f(x)-y */ Nlm_FloatHi dfx; /* Df(x) */ Int4 j; /* loop index */ Nlm_FloatHi fp, fq; /* Preliminary checks for improper bracketing and end-point root. */ if ((fp = F(p)) == 0.) return p; if ((fq = F(q)) == 0.) return q; if ((fp > 0. && fq > 0.) || (fp < 0. && fq < 0.)) { ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_INVAL,"NRBis: root not bracketed"); /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_INVAL, "root not bracketed"));**/ return HUGE_VAL; } /* Swaps end-points if necessary to make F(p)<0. 0.) { temp = p; p = q; q = temp; } /* Set up the Bisection & Newton-Raphson iteration. */ if ((x-p) * (x-q) > 0.) x = 0.5*(p+q); dxold = dx = p-q; for (j = 1; j <= NRBIS_ITMAX; ++j) { fx = F(x); if (fx == 0.) /* Check for termination. */ return x; if (fx < 0.) p = x; else q = x; dfx = DF(x); /* Check: root out of bounds or bisection faster than Newton-Raphson? */ if ((dfx*(x-p)-fx)*(dfx*(x-q)-fx) >= 0. || 2.*ABS(fx) > ABS(dfx*dx)) { dx = dxold; /* Bisect */ dxold = 0.5*(p-q); x = 0.5*(p+q); if (ABS(dxold) <= tol) return x; } else { dx = dxold; /* Newton-Raphson */ dxold = fx/dfx; x -= dxold; if (ABS(dxold) < tol && F(x-SIGN(dxold)*tol)*fx < 0.) return x; } } ErrPostEx(SEV_WARNING,E_Math,ERR_NCBIMATH_ITER,"NRBis: iterations > NRBIS_ITMAX"); /**ERRPOST((CTX_NCBIMATH, ERR_NCBIMATH_ITER, "iterations > NRBIS_ITMAX"));**/ return HUGE_VAL; } #undef F /* clean-up */ #undef DF /* clean-up */ /* 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 Nlm_FloatHi LIBCALL Nlm_RombergIntegrate(Nlm_FloatHi (LIBCALL *f) (Nlm_FloatHi,Nlm_VoidPtr), Nlm_VoidPtr fargs, Nlm_FloatHi p, Nlm_FloatHi q, Nlm_FloatHi eps, Int4 epsit, Int4 itmin) { Nlm_FloatHi romb[ROMBERG_ITMAX]; /* present list of Romberg values */ Nlm_FloatHi h; /* mesh-size */ Int4 i, j, k, npts; long n; /* 4^(error order in romb[i]) */ Int4 epsit_cnt = 0, epsck; register Nlm_FloatHi y; register Nlm_FloatHi x; register Nlm_FloatHi 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; } /* Nlm_Gcd(a, b) Return the greatest common divisor of a and b. Adapted 8-15-90 by WRG from code by S. Altschul. */ long LIBCALL Nlm_Gcd(register long a, register long b) { register 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 LIBCALL Nlm_Nint(register Nlm_FloatHi 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 */ Nlm_FloatHi LIBCALL Nlm_Powi(Nlm_FloatHi x, Int4 n) /* power */ { Nlm_FloatHi 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; } /* Additive random number generator Modelled after "Algorithm A" in Knuth, D. E. (1981). The art of computer programming, volume 2, page 27. 7/26/90 WRG */ static long state[33] = { (long)0xd53f1852, (long)0xdfc78b83, (long)0x4f256096, (long)0xe643df7, (long)0x82c359bf, (long)0xc7794dfa, (long)0xd5e9ffaa, (long)0x2c8cb64a, (long)0x2f07b334, (long)0xad5a7eb5, (long)0x96dc0cde, (long)0x6fc24589, (long)0xa5853646, (long)0xe71576e2, (long)0xdae30df, (long)0xb09ce711, (long)0x5e56ef87, (long)0x4b4b0082, (long)0x6f4f340e, (long)0xc5bb17e8, (long)0xd788d765, (long)0x67498087, (long)0x9d7aba26, (long)0x261351d4, (long)0x411ee7ea, (long)0x393a263, (long)0x2c5a5835, (long)0xc115fcd8, (long)0x25e9132c, (long)0xd0c6e906, (long)0xc2bc5b2d, (long)0x6c065c98, (long)0x6e37bd55 }; #define r_off 12 static long *rJ = &state[r_off], *rK = &state[DIM(state)-1]; void LIBCALL Nlm_RandomSeed(long x) { register Int4 i; state[0] = x; /* linear congruential initializer */ for (i=1; i>1)&0x7fffffff; /* discard the least-random bit */ }