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