1 // $Id: mathx.cpp,v 1.53 2011/04/23 02:02:49 bobgian Exp $
2
3 /*
4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein
5
6 This software is distributed free of charge for non-commercial use
7 and is copyrighted. Of course, we do not guarantee that the software
8 works, and are not responsible for any damage you may cause or have.
9 */
10
11 #include <cassert>
12 #include <cmath>
13 #include <iostream>
14 #include <numeric>
15 #include <stdio.h>
16
17 #include "constants.h" // for use of EXPMIN in UnderFlowExp()
18 #include "definitions.h"
19 #include "errhandling.h"
20 #include "mathx.h"
21 #include "registry.h"
22 #include "runreport.h"
23 #include "stringx.h"
24 #include "vectorx.h"
25
26 #ifdef DMALLOC_FUNC_CHECK
27 #include "/usr/local/include/dmalloc.h"
28 #endif
29
30 #ifdef HAVE_LGAMMA
31 #define LGAMMA lgamma
32 #else
33 #define LGAMMA mylgamma
34 #endif
35
36 using namespace std;
37
38 //------------------------------------------------------------------------------------
39
40 const double EigenCalculator::accuracy = 1e-15;
41 const double SMALL_v0 = 0.01; // erynes heuristic used by BesselK and DvBesselK
42 const double NUDGE_AMOUNT = 1.0e-09; // erynes heuristic used by BesselK and DvBesselK
43
44 //------------------------------------------------------------------------------------
45
46 double mylgamma (double z);
47
48 //------------------------------------------------------------------------------------
49 // Calculation of rate values following a gamma distribution for
50 // given probability values.
51
gamma_rates(double alpha,long int ratenum)52 DoubleVec1d gamma_rates(double alpha, long int ratenum)
53 {
54 vector<double> values;
55 double x, low, mid, high, xlow, xhigh, freq, inc, mid0;
56 long int i, j;
57
58 values.reserve(ratenum);
59
60 x = 10.0;
61 inc = 1.0 / (double)ratenum;
62 freq = -inc / 2.0;
63 mid0 = incompletegamma(alpha, 10.0);
64
65 for (i = 0; i < ratenum; i++)
66 {
67 low = 0;
68 mid = mid0;
69 high = 1.0;
70 freq += inc;
71
72 if (freq < mid)
73 {
74 high = mid;
75 xlow = 0;
76 xhigh = 10.0;
77 x = 5.0;
78 }
79
80 else
81 {
82 low = mid;
83 xlow = 10.0;
84 xhigh = 1e10;
85 x = 1e5;
86 }
87
88 for (j = 0; j < 1000 && fabs(low - high) > 0.0001 && x > 0.000000001; j++)
89 {
90 mid = incompletegamma(alpha, x);
91 if (freq < mid)
92 {
93 high = mid;
94 xhigh = x;
95 x = (x + xlow) / 2.0;
96 }
97
98 else
99 {
100 low = mid;
101 xlow = x;
102 x = (x + xhigh) / 2.0;
103 }
104 }
105
106 if (x >= 10e10)
107 {
108 values.clear();
109 break;
110 }
111
112 values.push_back(x / alpha);
113 }
114
115 return values;
116 }
117
118 //------------------------------------------------------------------------------------
119
120 #define MIN(a,b) (((a)<(b)) ? (b) : (a))
121
my_gamma(double x)122 double my_gamma(double x)
123 {
124 double result = log_gamma(x);
125 if (result < EXPMAX)
126 return exp(log_gamma(x));
127 return EXP_OF_EXPMAX;
128 }
129
130 double
alnorm(double x,int up)131 alnorm (double x, int up)
132 {
133 /* Initialized data */
134 /* *** machine dependent constants ????????????? */
135 /*static */ double zero = 0.0;
136 /*static */ double a1 = 5.75885480458;
137 /*static */ double a2 = 2.62433121679;
138 /*static */ double a3 = 5.92885724438;
139 /*static */ double b1 = -29.8213557807;
140 /*static */ double b2 = 48.6959930692;
141 /*static */ double c1 = -3.8052e-8;
142 /*static */ double c2 = 3.98064794e-4;
143 /*static */ double c3 = -.151679116635;
144 /*static */ double c4 = 4.8385912808;
145 /*static */ double c5 = .742380924027;
146 /*static */ double one = 1.0;
147 /*static */ double c6 = 3.99019417011;
148 /*static */ double d1 = 1.00000615302;
149 /*static */ double d2 = 1.98615381364;
150 /*static */ double d3 = 5.29330324926;
151 /*static */ double d4 = -15.1508972451;
152 /*static */ double d5 = 30.789933034;
153 /*static */ double half = 0.5;
154 /*static */ double ltone = 7.0;
155 /*static */ double utzero = 18.66;
156 /*static */ double con = 1.28;
157 /*static */ double p = .398942280444;
158 /*static */ double q = .39990348504;
159 /*static */ double r = .398942280385;
160
161 /*static */ double y, result;
162
163 if (x < zero)
164 {
165 up = !up;
166 x = -x;
167 }
168 if (x <= ltone || (up && x <= utzero))
169 {
170 y = half * x * x;
171 if (x > con)
172 {
173 result =
174 r * exp (-y) / (x + c1 +
175 d1 / (x + c2 +
176 d2 / (x + c3 +
177 d3 / (x + c4 +
178 d4 / (x + c5 +
179 d5 / (x + c6))))));
180 return ((!up) ? one - result : result);
181 }
182 result =
183 half - x * (p - q * y / (y + a1 + b1 / (y + a2 + b2 / (y + a3))));
184 return ((!up) ? one - result : result);
185 }
186 else
187 {
188 return ((!up) ? 1.0 : 0.0);
189 }
190 /*fake */ return -99;
191 } // alnorm
192
193 // The "complementary error function."
erfc(double x)194 double erfc(double x)
195 {
196 if (x >= 0.0)
197 return incompletegamma(0.5, x*x);
198 return 2.0 - incompletegamma(0.5, x*x); // a mathematical fact
199 }
200
201 // ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3
202 // Computation of the Incomplete Gamma Integral
203 // Auxiliary functions required: LogG() = logarithm of the gamma function,
204 // and Tail() = algorithm AS66
205 // In Mathematica, this is GammaRegularized[a,x] == Gamma[a,x]/Gamma[a].
206 // erynes note: The code below implements
207 // incompletegamma(a,x) = (1/G(a)) times the integral of exp(-t)*t^(a-1)
208 // from t = x to t = infinity for a > 0, where G(a) is the gamma function
209 // (if "a" is a positive integer, then G(a) = (a-1)!).
210 // erynes note number 2:
211 // incompletegamma(1/2,x^2) == erfc(x) for x >= 0,
212 // where erfc(x) is the complementary error function.
incompletegamma(double alpha,double x)213 double incompletegamma (double alpha, double x)
214 {
215 double gama, d_1, d_2, d_3;
216 /*static */ double a, b, c, an, rn;
217 /*static */ double pn1, pn2, pn3, pn4, pn5, pn6, arg;
218
219 gama = 0.0;
220 /* Check that we have valid values for X and P */
221 if (alpha <= 0.0 || x < 0.0)
222 {
223 logic_error e("Failure in incomplete-gamma calculation");
224 throw e;
225 }
226 if (fabs (x) < DBL_EPSILON)
227 return gama;
228
229 /* Use a normal approximation if P > PLIMIT */
230 if (alpha > 1e3)
231 {
232 pn1 =
233 sqrt (alpha) * 3.0 * (pow (x / alpha, (1.0 / 3.0)) + 1.0 / (alpha * 9.0) -
234 1.0);
235 gama = alnorm(pn1, false);
236 return gama;
237 }
238
239 /* If X is extremely large compared to P then set GAMMAD = 1 */
240 if (x > 1e8)
241 {
242 gama = 1.0;
243 return gama;
244 }
245
246 if (x <= 1.0 || x < alpha)
247 {
248 /* Use Pearson's series expansion. */
249 /* (Note that P is not large enough to force overflow in lgamma()). */
250 arg = alpha * log (x) - x - LGAMMA (alpha + 1.0);
251 c = 1.0;
252 gama = 1.0;
253 a = alpha;
254 while (c > 1e-14)
255 {
256 a += 1.0;
257 c = c * x / a;
258 gama += c;
259 }
260 arg += log (gama);
261 gama = 0.0;
262 if (arg >= -88.0)
263 {
264 gama = exp(arg);
265 }
266
267 }
268 else
269 {
270 /* Use a continued fraction expansion */
271 arg = alpha * log (x) - x - LGAMMA (alpha);
272 a = 1.0 - alpha;
273 b = a + x + 1.0;
274 c = 0.0;
275 pn1 = 1.0;
276 pn2 = x;
277 pn3 = x + 1.0;
278 pn4 = x * b;
279 gama = pn3 / pn4;
280 for (;;)
281 {
282 a += 1.0;
283 b += 2.0;
284 c += 1.0;
285 an = a * c;
286 pn5 = b * pn3 - an * pn1;
287 pn6 = b * pn4 - an * pn2;
288 if (fabs (pn6) > 0.0)
289 {
290 rn = pn5 / pn6;
291 /* Computing MIN */
292 d_2 = 1e-14, d_3 = rn * 1e-14;
293 if ((d_1 = gama - rn, fabs (d_1)) <= MIN (d_2, d_3))
294 {
295 arg += log (gama);
296 gama = 1.0;
297 if (arg >= -88.0)
298 {
299 gama = 1.0 - exp (arg);
300 }
301 return gama;
302 }
303 gama = rn;
304 }
305 pn1 = pn3;
306 pn2 = pn4;
307 pn3 = pn5;
308 pn4 = pn6;
309 if (fabs (pn5) >= 1e37)
310 {
311 /* Re-scale terms in continued fraction if terms are large */
312 pn1 /= 1e37;
313 pn2 /= 1e37;
314 pn3 /= 1e37;
315 pn4 /= 1e37;
316 }
317 }
318 }
319 return gama;
320 } // incompletegamma()
321
322 //------------------------------------------------------------------------------------
323 // Uses Lanczos-type approximation to ln(gamma) for z > 0.
324 // Reference: Lanczos, C. 'A precision approximation of the gamma function',
325 // J. SIAM Numer. Anal., B, 1, 86-96, 1964.
326 // Accuracy: About 14 significant digits except for small regions
327 // in the vicinity of 1 and 2.
328 // Programmer: Alan Miller CSIRO Division of Mathematics & Statistics
329 // Latest revision - 17 April 1988
330
log_gamma(double z)331 double log_gamma(double z)
332 {
333 if (z <= 0.0)
334 return DBL_MAX; // This will kill the receiving calculation.
335
336 long int i;
337 double result, denom;
338
339 double a[] = { 1.659470187408462e-7, 9.934937113930748e-6,
340 -0.1385710331296526, 12.50734324009056,
341 -176.6150291498386, 771.3234287757674,
342 -1259.139216722289, 676.5203681218835 };
343
344 result = 0.0;
345 denom = z + 7.0;
346 for (i = 0; i < 8; i++)
347 {
348 result += a[i] / denom;
349 denom -= 1.0;
350 }
351
352 result += 0.9999999999995183;
353 result = log(result) - (z + 6.5) + (z - 0.5) * log(z + 6.5) + 0.9189385332046727;
354 return result;
355 }
356
357 // The function for chi is:
358 // f(x) = {1/[2^(df/2) * gamma(df/2)]} * x^{[(df/2)-1] * e^(-x/2)}
359
360 double
find_chi(long int df,double prob)361 find_chi (long int df, double prob)
362 {
363 double a, b, m;
364 double xb = 200.0;
365 double xa = 0.0;
366 double xm = 5.0;
367 a = probchi (df, xa);
368 m = probchi (df, xm);
369 b = probchi (df, xb);
370 while (fabs (m - prob) > EPSILON)
371 {
372 if (m < prob)
373 {
374 b = m;
375 xb = xm;
376 }
377 else
378 {
379 a = m;
380 xa = xm;
381 }
382 xm = (-(b * xa) + prob * xa + a * xb - prob * xb) / (a - b); //(xa + xb)/2.0;
383
384 m = probchi (df, xm);
385 }
386 return xm;
387 }
388
389 double
probchi(long int df,double chi)390 probchi (long int df, double chi)
391 {
392 double prob;
393 double v = ((double) df) / 2.0;
394 if (chi > DBL_EPSILON && v > DBL_EPSILON)
395 {
396 //lg = EXP (LGAMMA (v));
397 prob = 1.0 - incompletegamma (v, chi / 2.0);
398 }
399 else
400 {
401 prob = 1.0;
402 // printf("prob=%f v=%f chi=%f lg(v/2)=%f ig(chi/2,v/2)=%f\n",
403 // prob,v,chi,lg, incompletegamma(chi/2.0,v/2.0));
404 }
405
406 return prob;
407 }
408
409 // Uses Lanczos-type approximation to ln(gamma) for z > 0.
410 // Reference:
411 // Lanczos, C. 'A precision approximation of the gamma
412 // function', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
413 // Accuracy: About 14 significant digits except for small regions
414 // in the vicinity of 1 and 2.
415 // Programmer: Alan Miller
416 // CSIRO Division of Mathematics & Statistics
417 // Latest revision - 17 April 1988
418 // translated and modified into C by Peter Beerli 1997
419
420 double
mylgamma(double z)421 mylgamma (double z)
422 {
423 double a[9] = { 0.9999999999995183, 676.5203681218835,
424 -1259.139216722289, 771.3234287757674, -176.6150291498386,
425 12.50734324009056, -0.1385710331296526, 9.934937113930748e-6,
426 1.659470187408462e-7
427 };
428 double lnsqrt2pi = 0.9189385332046727;
429 double result;
430 long int j;
431 double tmp;
432 if (z <= 0.0)
433 {
434 return DBL_MAX; //this will kill the receiving calculation
435 }
436 result = 0.0;
437 tmp = z + 7.0;
438 for (j = 9; j >= 2; --j)
439 {
440 result += a[j - 1] / tmp;
441 tmp -= 1.0;
442 }
443 result += a[0];
444 result = log (result) + lnsqrt2pi - (z + 6.5) + (z - 0.5) * log (z + 6.5);
445 return result;
446 } // lgamma
447
448 //------------------------------------------------------------------------------------
449
SafeDivide(double num,double denom)450 double SafeDivide(double num, double denom)
451 {
452 if (denom)
453 {
454 return num/denom;
455 }
456 else
457 {
458 return num > 0.0 ? DBL_BIG : -DBL_BIG;
459 }
460
461 } // SafeDivide
462
463 //------------------------------------------------------------------------------------
464
logfac(long int n)465 double logfac (long int n)
466 {
467 /* log(n!) values were calculated with Mathematica
468 with a precision of 30 digits */
469 switch (n)
470 {
471 case 0:
472 return 0.0;
473 case 1:
474 return 0.0;
475 case 2:
476 return 0.693147180559945309417232121458;
477 case 3:
478 return 1.791759469228055000812477358381;
479 case 4:
480 return 3.1780538303479456196469416013;
481 case 5:
482 return 4.78749174278204599424770093452;
483 case 6:
484 return 6.5792512120101009950601782929;
485 case 7:
486 return 8.52516136106541430016553103635;
487 case 8:
488 return 10.60460290274525022841722740072;
489 case 9:
490 return 12.80182748008146961120771787457;
491 case 10:
492 return 15.10441257307551529522570932925;
493 case 11:
494 return 17.50230784587388583928765290722;
495 case 12:
496 return 19.98721449566188614951736238706;
497 case 13:
498 return 22.5521638531234228855708498286;
499 case 14:
500 return 25.1912211827386815000934346935;
501 case 15:
502 return 27.8992713838408915660894392637;
503 case 16:
504 return 30.6718601060806728037583677495;
505 case 17:
506 return 33.5050734501368888840079023674;
507 case 18:
508 return 36.3954452080330535762156249627;
509 case 19:
510 return 39.3398841871994940362246523946;
511 case 20:
512 return 42.3356164607534850296598759707;
513 case 21:
514 return 45.3801388984769080261604739511;
515 case 22:
516 return 48.4711813518352238796396496505;
517 case 23:
518 return 51.6066755677643735704464024823;
519 case 24:
520 return 54.7847293981123191900933440836;
521 case 25:
522 return 58.0036052229805199392948627501;
523 case 26:
524 return 61.2617017610020019847655823131;
525 case 27:
526 return 64.5575386270063310589513180238;
527 case 28:
528 return 67.8897431371815349828911350102;
529 case 29:
530 return 71.2570389671680090100744070426;
531 case 30:
532 return 74.6582363488301643854876437342;
533 default:
534 return log(factorial(static_cast<double>(n)));
535 //return LGAMMA (n + 1.0);
536 }
537 } // logfac
538
factorial(double number)539 double factorial(double number)
540 {
541 double temp;
542
543 if(number <= 1.0) return 1.0;
544
545 temp = number * factorial(number - 1.0);
546 return temp;
547 }
548
549 //------------------------------------------------------------------------------------
550
UnderFlowExp(double pow)551 double UnderFlowExp(double pow)
552 {
553 return ((pow < EXPMIN) ? 0.0 : exp(pow));
554 } // UnderFlowExp
555
556 //------------------------------------------------------------------------------------
557
IsEven(long int n)558 bool IsEven(long int n)
559 {
560 // an even number divided by 2 is not truncated
561 return ((n / 2L) * 2L == n);
562 } // IsEven
563
564 //------------------------------------------------------------------------------------
565
ScaleLargestToZero(DoubleVec1d & unscaled)566 void ScaleLargestToZero(DoubleVec1d& unscaled)
567 {
568 double big(*std::max_element(unscaled.begin(),unscaled.end()));
569 for (unsigned long int wnum = 0; wnum < unscaled.size(); wnum++)
570 {
571 if (big <= EXPMIN)
572 {
573 assert(false); //Why is this? We probably have an error somewhere.
574 unscaled[wnum] = 0.0;
575 }
576 else if (unscaled[wnum] <= EXPMIN)
577 {
578 unscaled[wnum] = EXPMIN;
579 }
580 else
581 {
582 unscaled[wnum] -= big;
583 }
584 }
585 }
586
ScaleToSumToOne(DoubleVec1d & vec)587 void ScaleToSumToOne(DoubleVec1d& vec)
588 {
589 double sum = std::accumulate(vec.begin(), vec.end(), 0.0);
590 std::transform(vec.begin(),
591 vec.end(),
592 vec.begin(),
593 std::bind2nd(std::divides<double>(),sum));
594 }
595
596 // AddValsOfLogs takes vectors that are logs, scales and converts them to
597 // normal values, adds them, takes their logs again, and re-scales them back.
AddValsOfLogs(DoubleVec1d vec1,DoubleVec1d vec2)598 DoubleVec1d AddValsOfLogs(DoubleVec1d vec1, DoubleVec1d vec2)
599 {
600 assert(vec1.size() == vec2.size());
601 double big(*std::max_element(vec1.begin(), vec1.end()));
602 big = std::max(big, *std::max_element(vec2.begin(), vec2.end()));
603 for (unsigned long int wnum=0; wnum<vec1.size(); wnum++)
604 {
605 if (big <= EXPMIN)
606 {
607 assert(false); //Why is this? We probably have an error somewhere.
608 vec1[wnum] = EXPMIN;
609 vec2[wnum] = EXPMIN;
610 }
611 else
612 {
613 if (vec1[wnum] <= EXPMIN)
614 {
615 vec1[wnum] = EXPMIN;
616 }
617 else
618 {
619 vec1[wnum] -= big;
620 }
621 if (vec2[wnum] <= EXPMIN)
622 {
623 vec2[wnum] = EXPMIN;
624 }
625 else
626 {
627 vec2[wnum] -= big;
628 }
629 }
630 }
631 vec1 = SafeExp(vec1);
632 vec2 = SafeExp(vec2);
633 DoubleVec1d retvec = vec1;
634 std::transform(vec2.begin(),
635 vec2.end(),
636 retvec.begin(),
637 retvec.begin(),
638 std::plus<double>());
639 retvec = SafeLog(retvec);
640 for (unsigned long int wnum = 0; wnum<retvec.size(); wnum++)
641 {
642 if ((retvec[wnum] <= EXPMIN) || (retvec[wnum] + big <=EXPMIN))
643 {
644 retvec[wnum] = EXPMIN;
645 }
646 else if ((retvec[wnum] >= EXPMAX) || (retvec[wnum] + big >= EXPMAX))
647 {
648 retvec[wnum] = EXPMAX;
649 }
650 else
651 {
652 retvec[wnum] += big;
653 }
654 }
655 return retvec;
656 }
657
Coeffs(double x,double y)658 std::pair<double,double> EigenCalculator::Coeffs(double x, double y)
659 // cosine and sine of angle between the origin and (x,y)
660 // pair.first is cosine, pair.second is sine
661 {
662 double root = sqrt(pow(x,2.0)+pow(y,2.0));
663 std::pair<double,double> cs;
664 if (root < accuracy)
665 {
666 cs.first = 1.0;
667 cs.second = 0.0;
668 }
669 else
670 {
671 cs.first = x/root;
672 cs.second = y/root;
673 }
674 return cs;
675 } // Coeffs
676
677 //------------------------------------------------------------------------------------
678
Givens(DoubleVec2d & a,long int x,long int y,long int size,const std::pair<double,double> cs,bool userow)679 void EigenCalculator::Givens(DoubleVec2d& a, long int x, long int y, long int size,
680 const std::pair<double,double> cs, bool userow)
681 // Givens method. Modifies input matrix.
682 {
683 long int k;
684 for (k = 0; k < size; ++k)
685 {
686 if (userow)
687 {
688 double d = cs.first * a[x][k] + cs.second * a[y][k];
689 a[y][k] = cs.first * a[y][k] - cs.second * a[x][k];
690 a[x][k] = d;
691 }
692 else
693 {
694 double d = cs.first * a[k][x] + cs.second * a[k][y];
695 a[k][y] = cs.first * a[k][y] - cs.second * a[k][x];
696 a[k][x] = d;
697 }
698 }
699 // cases: we pass array.size()=4
700 // Pascal: k from 1 to 4 (4 times)
701 // C: k from 0 to 3 (4 times) (fine)
702 // case: we pass i (at its top value)
703 // Pascal: k from 1 to 4 (4 times)
704 // C: k from 0 to 2 (3 times) (not fine)
705 } // Givens
706
707 //------------------------------------------------------------------------------------
708
Tridiag(DoubleVec2d & a,DoubleVec2d & eigvecs)709 void EigenCalculator::Tridiag(DoubleVec2d& a, DoubleVec2d& eigvecs)
710 {
711 unsigned long int i, j;
712 std::pair<double,double> cs;
713 for (i = 1; i < a.size()-1; ++i)
714 {
715 for (j = i+1; j < a.size(); ++j)
716 {
717 cs = Coeffs(a[i-1][i], a[i-1][j]);
718 Givens(a, i, j, a.size(), cs, true);
719 Givens(a, i, j, a.size(), cs, false);
720 Givens(eigvecs, i, j, eigvecs.size(), cs, true);
721 }
722 }
723 } // Tridiag
724
725 //------------------------------------------------------------------------------------
726
Shiftqr(DoubleVec2d & a,DoubleVec2d & eigvecs)727 void EigenCalculator::Shiftqr(DoubleVec2d& a, DoubleVec2d& eigvecs)
728 {
729 long int i;
730 for (i = a.size()-1; i > 0; --i)
731 {
732 do {
733 double d = sqrt(pow(a[i-1][i-1] - a[i][i], 2.0) +
734 pow(a[i][i-1], 2.0));
735 double approx = a[i-1][i-1] + a[i][i];
736 if (a[i][i] < a[i-1][i-1])
737 approx = (approx-d)/2.0;
738 else
739 approx = (approx+d)/2.0;
740 long int j;
741 for (j = 0; j <= i; ++j)
742 a[j][j] = a[j][j] - approx;
743 for (j = 0; j <= i-1; ++j)
744 {
745 std::pair<double,double> cs = Coeffs(a[j][j], a[j+1][j]);
746 // in the following two calls, Pascal's i has been
747 // converted to i+1 because it is standing in for a
748 // size.
749 Givens(a, j, j+1, i+1, cs, true);
750 Givens(a, j, j+1, i+1, cs, false);
751 Givens(eigvecs, j, j+1, eigvecs.size(), cs, true);
752 }
753 for (j = 0; j <= i; ++j)
754 a[j][j] += approx;
755 } while (fabs(a[i][i-1]) > accuracy);
756 }
757 } // Shiftqr
758
759 //------------------------------------------------------------------------------------
760
Eigen(DoubleVec2d a)761 std::pair<DoubleVec1d, DoubleVec2d> EigenCalculator::Eigen(DoubleVec2d a)
762 // return a pair containing eigenvalues and eigenvectors
763 {
764 unsigned long int i;
765 DoubleVec2d eigvecs;
766 for (i = 0; i < a.size(); ++i)
767 {
768 // ones along the diagonal, zeroes elsewhere
769 DoubleVec1d row(a.size(), 0.0);
770 row[i] = 1.0;
771 eigvecs.push_back(row);
772 }
773
774 Tridiag(a, eigvecs);
775 Shiftqr(a, eigvecs);
776 DoubleVec1d eigvals(a.size(), 0.0);
777 for (i = 0; i < a.size(); ++i)
778 {
779 eigvals[i] = a[i][i];
780 }
781 return std::make_pair<DoubleVec1d, DoubleVec2d>(eigvals, eigvecs);
782 } // Eigen
783
784 //------------------------------------------------------------------------------------
785
DotProduct(const DoubleVec2d & first,const DoubleVec2d & second,DoubleVec2d & answer)786 void EigenCalculator::DotProduct(const DoubleVec2d& first, const DoubleVec2d& second,
787 DoubleVec2d& answer)
788 // dot product of first and second put into PRE-EXISTING answer!
789 // second has been PRE-TRANSPOSED!!
790 {
791 // should be square
792 assert(first.size() == first[0].size());
793 // and all the same size!
794 assert(first.size() == second.size());
795 assert(first.size() == answer.size());
796 double initial = 0.0;
797
798 long int i, j, n = first.size();
799 for (i = 0; i < n; ++i)
800 {
801 for (j = 0; j < n; ++j)
802 {
803 answer[i][j] = std::inner_product(first[i].begin(), first[i].end(), second[j].begin(), initial);
804 }
805 }
806
807 } // DotProduct
808
809 //------------------------------------------------------------------------------------
810
Invert(const DoubleVec2d & src)811 DoubleVec2d Invert(const DoubleVec2d& src)
812 // Invert square matrix via Gauss-Jordan reduction
813 {
814 // copy the target vector
815 DoubleVec2d a(src);
816
817 // invert it in place (that's what I have code for!)
818 unsigned long int i, j, k;
819 unsigned long int n = a.size();
820 double temp;
821
822 for (i = 0; i < n; ++i)
823 {
824 // DEBUG: need to do something if matrix is singular!
825 temp = 1.0 / a[i][i];
826 a[i][i] = 1.0;
827 for (j = 0; j < n; ++j)
828 {
829 a[i][j] *= temp;
830 }
831 for (j = 0; j < n; ++j)
832 {
833 if (j != i)
834 {
835 temp = a[j][i];
836 a[j][i] = 0.0;
837 for (k = 0; k < n; ++k)
838 {
839 a[j][k] -= temp * a[i][k];
840 }
841 }
842 }
843 }
844 return a;
845 } // Invert
846
847 //------------------------------------------------------------------------------------
848
Transpose(const DoubleVec2d & src)849 DoubleVec2d Transpose(const DoubleVec2d& src)
850 // Transpose a square matrix
851 {
852 DoubleVec2d a(src);
853 unsigned long int i, j, n=a.size();
854 for (i = 0; i < n; ++i)
855 {
856 for (j = 0; j < n; ++j)
857 {
858 a[i][j] = src[j][i];
859 }
860 }
861 return a;
862 } // Transpose
863
864 //------------------------------------------------------------------------------------
865
866 // This is a free function for approximate comparison of doubles.
867 // It is used to help find the correct profile modifier value in
868 // the face of rounding error.
869
CloseEnough(double a,double b)870 bool CloseEnough(double a, double b)
871 {
872 if (a == 0)
873 {
874 if (b == 0) return true;
875 else return false;
876 }
877 if (a<0)
878 {
879 if (b<0)
880 {
881 a = -a;
882 b = -b;
883 }
884 else return false;
885 }
886 double la = log(a);
887 double lb = log(b);
888 double test = fabs(la - lb);
889 if (test < EPSILON) return true;
890 else return false;
891
892 } // CloseEnough
893
SafeExpAndSum(const DoubleVec1d & src)894 double SafeExpAndSum(const DoubleVec1d& src)
895 {
896 assert(!src.empty()); // don't call me on an empty vector!
897 double sum(0.0), biggest(*(std::max_element(src.begin(),src.end())));
898 DoubleVec1d::const_iterator lns;
899 for(lns = src.begin(); lns != src.end(); ++lns)
900 {
901 if ((*lns) - biggest > EXPMIN)
902 sum += exp((*lns) - biggest);
903 }
904
905 return SafeLog(sum) + biggest;
906
907 } // SafeExpAndSum
908
SafeExp(const DoubleVec1d & src)909 DoubleVec1d SafeExp(const DoubleVec1d& src)
910 {
911 assert(!src.empty()); // don't call me on an empty vector!
912 DoubleVec1d retvec;
913 DoubleVec1d::const_iterator lns;
914 for(lns = src.begin(); lns != src.end(); ++lns)
915 {
916 retvec.push_back(SafeExp(*lns));
917 }
918
919 return retvec;
920
921 } // SafeSumOfLogsToBeExp
922
SafeExp(double src)923 double SafeExp(double src)
924 {
925 if (src <= EXPMIN)
926 {
927 return 0.0;
928 }
929 else if (src >= EXPMAX)
930 {
931 return DBL_BIG;
932 }
933 return exp(src);
934
935 } // SafeExp
936
937 // Technique to avoid overflow and underflow of a*exp(x).
938 // In the description that follows, we set a = -a for a < 0,
939 // and we use "M" to represent the maximum feasible number.
940 // This means that numbers > M overflow, and numbers < 1/M underflow.
941 //
942 // First we consider the case of overflow.
943 // a*exp(x) will not overflow if a*exp(x) <= M, meaning exp(x) <= M/a.
944 // Defining EXPMAX = log(M), we see a*exp(x) will not overflow if
945 // x <= EXPMAX - log(a).
946 // For the case in which exp(x) overflows but a*exp(x) does not, we compute
947 // a*exp(x) == (a*exp(EXPMAX)) * exp(x - EXPMAX), so that neither of the two
948 // terms on the right-hand side of this equation overflows.
949 // (We use the pre-computed EXP_OF_EXPMAX instead of computing exp(EXPMAX).)
950 // In practical terms, if EXPMAX = 200 and x = 212 so that exp(212) overflows,
951 // we can compute a*exp(x) if fabs(a) < 6.14e-6.
952 //
953 // Underflow is analogous to overflow. In this case, we have
954 // fabs(a*exp(x)) < 1.0/M becoming identically 0 for sufficiently large M.
955 // If "a" is sufficiently greater than 1.0, then the product will not underflow.
956 // (Again, we set a = -a for a > 0 to avoid writing fabs(a).)
957 // For a*exp(x) to avoid underflow, we have the condition
958 // a*exp(x) >= 1.0/M, meaning exp(x) >= 1.0/(M*a), or x > log(1.0/(M*a)), so
959 // a*exp(x) will not underflow if
960 // x >= EXPMIN - log(a).
961 // For the case in which exp(x) underflows but a*exp(x) does not, we compute
962 // a*exp(x) == (a*exp(EXPMIN)) * exp(x - EXPMIN), so that neither of the two
963 // terms on the right-hand side of this equation underflows.
964 // (We use the pre-computed EXP_OF_EXPMIN instead of computing exp(EXPMIN).)
965 //
966 // erynes 2003/11/21 -- devised and implemented this function
SafeProductWithExp(double a,double x)967 double SafeProductWithExp(double a, double x)
968 {
969 if (0.0 == a) // must reject this at the outset to avoid log(0) and 1/0
970 return 0.0; // keep going if 0 < a < DBL_EPSILON, in case x is large
971 bool a_is_negative = a < 0 ? true : false;
972 if (a_is_negative)
973 a *= -1.0;
974
975 if (x > 0.0)
976 {
977 if (x <= EXPMAX - log(a))
978 {
979 // computable
980 if (x <= EXPMAX)
981 return a_is_negative ? -a*exp(x) : a*exp(x);
982 // else 0.0 < a < 1.0
983 double a_exp_x1 = EXP_OF_EXPMAX / (1.0/a);
984 return a_is_negative ? -a_exp_x1*exp(x - EXPMAX)
985 : a_exp_x1*exp(x - EXPMAX);
986 }
987 return OverflowOfProductWithExp(a_is_negative ? -a : a, x);
988 }
989
990 // else x < 0.
991 if (x >= EXPMIN - log(a))
992 {
993 // computable
994 if (x >= EXPMIN)
995 return a_is_negative ? -a*exp(x) : a*exp(x);
996 // else a > 1.0
997 double a_exp_x1 = a * EXP_OF_EXPMIN;
998 return a_is_negative ? -a_exp_x1*exp(x - EXPMIN)
999 : a_exp_x1*exp(x - EXPMIN);
1000 }
1001 return UnderflowOfProductWithExp(a_is_negative ? -a : a, x);
1002
1003 } // double SafeProductWithExp(double a, double x)
1004
1005 // This function, when possible, computes exp(x1) - exp(x2)
1006 // without overflow or underflow, although exp(x1) and/or exp(x2)
1007 // may overflow or underflow individually.
1008 // The technique is as follows: Define "result" = exp(x1) - exp(x2).
1009 // Then
1010 // result == exp(x2) * ( exp(x1)/exp(x2) - 1.0 )
1011 // == exp(x2) * ( exp(x1 - x2) - 1.0 ),
1012 // and
1013 // log(result) == x2 + log(exp(x1 - x2) - 1.0),
1014 // so that the right-hand side of this equation is computable.
1015 // If EXPMIN < log(result) < EXPMAX, then we can return
1016 // result = exp(x2 + log(exp(x1 - x2) - 1.0).
1017 // Note that, in practice, for x1 > EXPMAX, x2 must be extremely
1018 // close to x1 in order for this function to be useful.
1019 // For example, for EXPMAX = 200 and x2 = 212 with exp(212) overflowing,
1020 // SafeExpDiff(x1, x2) will overflow and return EXP_OF_EXPMAX
1021 // unless 212.000006 >= x1 >= 211.999994.
1022 //
1023 // erynes 2003/11/21 -- implemented and half-devised this function
SafeExpDiff(double x1,double x2)1024 double SafeExpDiff(double x1, double x2)
1025 {
1026 if (x1 >= EXPMIN && x2 >= EXPMIN &&
1027 x1 <= EXPMAX && x2 <= EXPMAX)
1028 return exp(x1) - exp(x2);
1029
1030 bool result_is_negative = x2 > x1 ? true : false;
1031 double x_larger = result_is_negative ? x2 : x1,
1032 x_smaller = result_is_negative ? x1 : x2,
1033 arg_diff = x_larger - x_smaller;
1034
1035 if (arg_diff < DBL_EPSILON)
1036 return 0.0;
1037
1038 if (arg_diff > EXPMAX)
1039 {
1040 if (x_smaller < 0.0 && x_larger < EXPMAX)
1041 return result_is_negative ? -exp(x_larger) : exp(x_larger);
1042 return OverflowOfSafeExpDiff(x1, x2);
1043 }
1044
1045 if (arg_diff < EXPMIN)
1046 return UnderflowOfSafeExpDiff(x1, x2); // recall x_smaller < x_larger
1047
1048 double log_of_result = x_smaller + log(exp(arg_diff) - 1);
1049
1050 if (log_of_result > EXPMAX)
1051 return OverflowOfSafeExpDiff(x1, x2);
1052
1053 if (log_of_result < EXPMIN)
1054 return UnderflowOfSafeExpDiff(x1, x2);
1055
1056 return result_is_negative ? -exp(log_of_result) : exp(log_of_result);
1057 }
1058
1059 // This function is called when SafeExpDiff(x1, x2) overflows,
1060 // which means exp(x1) - exp(x2) overflows.
1061 // It can easily be modified to also track how often this overflows.
1062 // 2003/11/21 added by erynes
OverflowOfSafeExpDiff(double x1,double x2)1063 double OverflowOfSafeExpDiff(double x1, double x2)
1064 {
1065 RunReport& runreport = registry.GetRunReport();
1066 string msg="Overflow error: Attempted to compute exp(";
1067 msg += ToString(x1) + ") - exp(" + ToString(x2) + "). Returning "
1068 + ToString(x1 < x2 ? -EXP_OF_EXPMAX : EXP_OF_EXPMAX)
1069 + " = " + (x1 < x2 ? "-" : "") + "EXP_OF_EXPMAX. "
1070 + "(Further overflow errors of this type will not be reported.)\n";
1071 runreport.ReportOnce(msg, oncekey_OverflowOfSafeExpDiff, true);
1072
1073 return x1 < x2 ? -EXP_OF_EXPMAX : EXP_OF_EXPMAX;
1074 }
1075
1076 // This function is called when SafeExpDiff(x1, x2) underflows,
1077 // which means exp(x1) - exp(x2) underflows.
1078 // It can easily be modified to track how often this underflows.
1079 // 2003/11/21 added by erynes
UnderflowOfSafeExpDiff(double x1,double x2)1080 double UnderflowOfSafeExpDiff(double x1, double x2)
1081 {
1082 RunReport& runreport = registry.GetRunReport();
1083 string msg = "Underflow error: Attempted to compute exp(";
1084 msg += ToString(x1) + ") - exp(" + ToString(x2) + "). Returning 0. "
1085 + "(Further underflow errors of this type will not be reported.)\n";
1086 runreport.ReportOnce(msg, oncekey_UnderflowOfSafeExpDiff, false);
1087
1088 return 0.0;
1089 }
1090
1091 // This function is called when SafeProductWithExp(a, x) overflows,
1092 // which means a*exp(x) overflows.
1093 // It can easily be modified to also track how often a*exp(x) overflows.
1094 // 2003/11/21 added by erynes
OverflowOfProductWithExp(double a,double x)1095 double OverflowOfProductWithExp(double a, double x)
1096 {
1097 RunReport& runreport = registry.GetRunReport();
1098 string msg = "Overflow error: Attempted to compute ";
1099 msg += ToString(a) + " * exp(" + ToString(x) + "). Returning "
1100 + ToString(a < 0 ? -EXP_OF_EXPMAX : EXP_OF_EXPMAX)
1101 + " = " + (a < 0 ? "-" : "") + "EXP_OF_EXPMAX. "
1102 + "(Further overflow errors of this type will not be reported.)\n";
1103 runreport.ReportOnce(msg, oncekey_OverflowOfProductWithExp, true);
1104
1105 return a < 0 ? -EXP_OF_EXPMAX : EXP_OF_EXPMAX;
1106 }
1107
1108 // This function is called when SafeProductWithExp(a, x) underflows,
1109 // which means a*exp(x) underflows.
1110 // It can easily be modified to track how often a*exp(x) underflows.
1111 // 2003/11/21 added by erynes
UnderflowOfProductWithExp(double a,double x)1112 double UnderflowOfProductWithExp(double a, double x)
1113 {
1114 RunReport& runreport = registry.GetRunReport();
1115 string msg = "Underflow error: Attempted to compute ";
1116 msg += ToString(a) + " * exp(" + ToString(x) + "). Returning 0. "
1117 + "(Further underflow errors of this type will not be reported.)\n";
1118 runreport.ReportOnce(msg, oncekey_UnderflowOfProductWithExp, false);
1119
1120 return 0;
1121 }
1122
LogZero(const double x)1123 double LogZero(const double x)
1124 {
1125 return ((x > 0.0) ? log(x) : 0.0);
1126 }
1127
SafeLog(const double x)1128 double SafeLog(const double x)
1129 {
1130 return ((x > 0.0) ? log(x) : -DBL_BIG);
1131 }
1132
SafeLog(const DoubleVec1d src)1133 DoubleVec1d SafeLog(const DoubleVec1d src)
1134 {
1135 DoubleVec1d retvec;
1136 for (DoubleVec1d::const_iterator i=src.begin(); i != src.end(); ++i)
1137 {
1138 retvec.push_back(SafeLog(*i));
1139 }
1140 return retvec;
1141 }
1142
ChooseRandomFromWeights(DoubleVec1d & weights)1143 long int ChooseRandomFromWeights(DoubleVec1d& weights)
1144 {
1145 double sumweights = accumulate(weights.begin(),weights.end(),0.0);
1146 transform(weights.begin(),weights.end(),weights.begin(),
1147 bind2nd(divides<double>(),sumweights));
1148
1149 long int chosen = -1;
1150 double chance = registry.GetRandom().Float();
1151 while (chance > 0 && chosen < static_cast<long int>(weights.size()-1))
1152 {
1153 chance -= weights[++chosen];
1154 }
1155 return chosen;
1156 }
1157
1158 // ExpE1(x) -- added 2004/10/06 by erynes to calculate the expectation value
1159 // of "endtime" with positive growth
1160 //
1161 // Receives a positive real number x, and returns exp(x)*E1(x), where E1(x) is
1162 // the exponential integral function En(x) for n = 1. The function returns the
1163 // product of exp(x) and E1(x), rather than E1(x) or En(x) alone, because at
1164 // present we only use the product of exp(x) and E1(x), and since the most
1165 // common form of E1(x) includes an internal factor of exp(-x), we can avoid two
1166 // calls to exp() by allowing exp(x) and exp(-x) to cancel one another first.
1167 //
ExpE1(const double & x)1168 double ExpE1(const double& x)
1169 {
1170 if (x <= 0.0)
1171 {
1172 string msg = " Nonpositive number (" + ToString(x)
1173 + ") received by math function ExpE1(); this is illegal.";
1174 implementation_error e(msg);
1175 throw e;
1176 }
1177
1178 const double epsilon = 1.0e-07; // convergence criterion
1179 const unsigned long int max_iterations = 15UL; // surrender criterion
1180 double result;
1181
1182 if (x >= 1.0)
1183 {
1184 // Continued fraction representation: Lentz's algorithm.
1185 double a, b, c, d, h, delta;
1186 unsigned long int nIter(0UL);
1187 b = x + 1.0;
1188 c = 1.0/DBL_EPSILON;
1189 d = 1.0/b;
1190 h = d;
1191 for (unsigned long int i = 1; i <= max_iterations; i++)
1192 {
1193 a = -1.0*i*i;
1194 b += 2.0;
1195 d = 1.0/(a*d + b);
1196 c = b + a/c;
1197 delta = c*d;
1198 h *= delta;
1199 nIter++;
1200 if (fabs(delta - 1.0) < epsilon)
1201 return h;
1202 }
1203 // Failed to converge after nIter iterations.
1204 // Apparently this can only happen when x is very close to 1.
1205 // In that case, our result will be close enough.
1206 return h;
1207 }
1208
1209 else
1210 {
1211 //series representation
1212 double factor, term;
1213 unsigned long int nIter(0UL);
1214 result = -log(x) - EULERS_CONSTANT;
1215 factor = 1.0;
1216 for (unsigned long int i = 1; i <= max_iterations; i++)
1217 {
1218 factor *= -x/i;
1219 term = -factor/i;
1220 result += term;
1221 nIter++;
1222 if (fabs(term) < fabs(result) * epsilon)
1223 return SafeProductWithExp(result, x);
1224 }
1225 // Failed to converge after nIter iterations.
1226 // Apparently this can only happen when x is very close to 1.
1227 // In that case, our result will be close enough.
1228 return SafeProductWithExp(result, x);
1229 }
1230 }
1231
1232 // Function to compute the modified Bessel function of the second kind,
1233 // of order v, at argument x. It is also sometimes called the MacDonald
1234 // function. It is denoted in the literature as K(v,x), where "v" is
1235 // written as a subscript, rather than a parameter. Mathematica denotes
1236 // this function BesselK, so we'll do the same.
1237 //
1238 // Because K(v,x) is calculated iteratively, and because
1239 // dK(v,x)/dx is a simple function of both K(v,x) and K(v+1,x),
1240 // this function returns both of these results, since the latter
1241 // can be obtained essentially for free once the former is computed.
1242 // The return value of this function is K(v,x). K(v+1,x) is returned
1243 // by reference in argument "resultFor_vPlusOne."
1244 //
1245 // Implemented by erynes in July and October 2005, adapted from:
1246 // Shanjie Zhang and Jianming Jin, _Computation of Special Functions._
1247 // New York: Wiley-Interscience (1996).
1248 //
1249 // Contains an original approximation for x < 9 when v is very close
1250 // to an integer. No special approximation appears to be needed
1251 // for any other case, including when v is equal to an integer.
1252 //
BesselK(double v,double x,double & resultFor_vPlusOne)1253 double BesselK(double v, double x, double& resultFor_vPlusOne)
1254 {
1255 if (x <= 0.0)
1256 throw implementation_error("BesselK() received illegal argument: "
1257 + ToString(x) + ".");
1258
1259 bool orderIsNegative(false), mustUseTwo_v0s(false),
1260 first_v0_calculation(true);
1261 double result(0.0); // used ONLY if -1 < v < 0
1262
1263 if (v < 0.0)
1264 {
1265 if (v > -1.0)
1266 {
1267 // Generally, we calculate K(v0,x) and K(v0+1,x),
1268 // where v0 is the "fractional part" of v, and then
1269 // we calculate K(v,x) and K(v+1,x) iteratively.
1270 // If v <= -1, we use the property K(-v,x) = K(v,x),
1271 // and swap the results at the end. But if -1 < v < 0,
1272 // then v < 0 and v+1 > 0, so we have to compute K(v0,x)
1273 // twice--once for v (say, -0.7), and a separate time
1274 // for v+1 (say, 0.3).
1275 v = -v;
1276 mustUseTwo_v0s = true;
1277 }
1278 else
1279 {
1280 // K(-v, x) == K(v, x) and K(-v + 1, x) == K(v - 1, x).
1281 v = -v - 1.0;
1282 }
1283 orderIsNegative = true; // later, we'll swap K(v+1,x) and K(v,x)
1284 }
1285
1286 // For tiny x, K(v,x) = 0.5*gamma(v)*pow(0.5*x, -v).
1287 double gamma_vPlus1__DividedBy__DBL_BIG = my_gamma(v+1.0)/DBL_BIG;
1288 if (x <= 2.0*pow(gamma_vPlus1__DividedBy__DBL_BIG*0.5, 1.0/(v+1.0)))
1289 {
1290 // K(v+1,x) overflows, and hence K(v,x) overflows too.
1291 resultFor_vPlusOne = DBL_BIG;
1292 return DBL_BIG;
1293 }
1294
1295 // For huge x, K(v,x) = sqrt(PI/(2x))*exp(-x).
1296 if (x > -EXPMIN || 0.0 == exp(-x)*sqrt(0.5*PI/x))
1297 {
1298 // K(v+1,x) underflows, for all v.
1299 resultFor_vPlusOne = 0.0;
1300 return 0.0;
1301 }
1302
1303 const double EPSILON_ZJ = 1.0e-15; // Zhang & Jin's convergence criterion
1304 unsigned long int n = static_cast<unsigned long int>(floor(v)); // closest integral order
1305 double v0 = v - n; // the "fractional part" of the order
1306 const unsigned long int MAX_ITERATIONS = 50; // Zhang & Jin's heuristic cutoff
1307
1308 double half_x_squared = 0.25*x*x; // avoid repeated recalculations
1309 double Kv0 = 0.0; // == K(v0, x)
1310 double Kv0plus1 = 0.0; // == K(v0+1, x)
1311
1312 // First, calculate K(v0, x) and K(v0+1, x).
1313
1314 Calculate_Kv0:
1315 if (x <= 9.0)
1316 {
1317 // Small argument: Use the series expansion.
1318 if (0.0 == v0)
1319 {
1320 // K(0,x) = -(ln(x/2) + EULERS_CONSTANT)*I(0,x)
1321 // + sum(coeff_i*term_i, from i=1 to i=large),
1322 //
1323 // where I(0,x) = modified Bessel fn. of the 1st kind of order 0
1324 // = sum(term_i, from i = 0 to i = large),
1325 // with
1326 //
1327 // term_i = ((x/2)^(2*i)) / (i!)^2
1328 //
1329 // and
1330 //
1331 // coeff_i = sum(1/j, from j=1 to j=i), or 1 if i=0.
1332 //
1333 // K(1,x) = (1/x) + (ln(x/2) + EULERS_CONSTANT)*I(1,x)
1334 // - (x/2)*sum(coeff_i_2*term_i_2, from i=1 to i=large),
1335 // where I(1,x) = modified Bessel fn. of the 1st kind of order 1
1336 // = (x/2)*sum(term_i_2, from i=1 to i=large),
1337 //
1338 // with
1339 //
1340 // term_i_2 = ((x/2)^(2*i)) / (i!*(i+1)!)
1341 //
1342 // and
1343 //
1344 // coeff_i_2 = sum(1/j, from j=1 to j=i) + (1/2)*(1/(i+1)), or 1/2 if i=0.
1345
1346 double prevKv0 = DBL_BIG, prevKv0plus1 = DBL_BIG;
1347 double term_A = log(0.5*x) + EULERS_CONSTANT;
1348 double term_i = 1.0, term_i_2 = 1.0;
1349 double coeff_i = 0.0, coeff_i_2 = 0.0;
1350 Kv0 = -term_A; // i=0 term
1351 Kv0plus1 = 1.0/x + term_A*0.5*x - 0.25*x; // i=0 term
1352
1353 for (unsigned long int i = 1; i <= MAX_ITERATIONS; i++)
1354 {
1355 term_i *= half_x_squared/(i*i);
1356 term_i_2 *= half_x_squared/(i*(i+1));
1357 coeff_i += 1.0/i;
1358 coeff_i_2 = coeff_i + 1.0/(2.0*(i+1.0));
1359 Kv0 += (-term_A + coeff_i)*term_i;
1360 Kv0plus1 += (term_A - coeff_i_2)*0.5*x*term_i_2;
1361 if (fabs((prevKv0 - Kv0)/Kv0) < EPSILON_ZJ &&
1362 fabs((prevKv0plus1 - Kv0plus1)/Kv0plus1) < EPSILON_ZJ)
1363 break; // convergence achieved
1364 prevKv0 = Kv0;
1365 prevKv0plus1 = Kv0plus1;
1366 }
1367 }
1368
1369 else // v0 > 0
1370 {
1371 // K(v0, x) = (PI/2)*(I(-v0, x) - I(v0, x))/sin(v0*PI),
1372 // where
1373 // I(v0,x) = modified Bessel fn. of the 1st kind of order v0
1374 // = sum(((x/2)^(2i+v0))/(i!*gamma(i+v0+1)),
1375 // from i=0 to i=large)
1376 //
1377 // and gamma(x+1) = x*gamma(x).
1378 //
1379 // For v0 close to 0 or 1, the above formula for K(v0,x) approaches 0/0,
1380 // and behaves quasi-randomly. An excellent example is K(1.999999,7).
1381 // Increasing the number of iterations and increasing the strictness
1382 // of the convergence criterion (i.e., decreasing EPSILON_ZJ)
1383 // does not seem to help. We thus employ the following approximation,
1384 // derived by erynes in autumn 2005.
1385
1386 if (v0 <= SMALL_v0 || v0 >= 1.0 - SMALL_v0)
1387 {
1388 double K0(0.0), K1(0.0); // K(0,x) and K(1,x)
1389 K0 = BesselK(0.0, x, K1);
1390 if (v0 >= 1.0 - SMALL_v0)
1391 {
1392 v0 = -(1.0 - v0); // calculate K(-eps, x), K(1 - eps, x), etc.
1393 n += 1; // An example: Redefining 2 + .95 to be 3 - .05, n=2 --> n=3.
1394 }
1395
1396 Kv0 = K0 * (v0*v0/(2.0*x + 1.0) + 1.0);
1397 Kv0plus1 = ((2.0*x-1.0)*(4.0*x+1.0)*K1/(2.0*(x+1.0)*(x+1.0)))*v0*v0 + (K0/x)*v0 + K1;
1398
1399 if (0 == n)
1400 {
1401 resultFor_vPlusOne = Kv0plus1;
1402 if (Kv0 >= DBL_BIG)
1403 {
1404 // Try to avoid returning "inf" in either variable.
1405 resultFor_vPlusOne = DBL_BIG;
1406 return DBL_BIG;
1407 }
1408 if (resultFor_vPlusOne >= DBL_BIG)
1409 resultFor_vPlusOne = DBL_BIG;
1410 return Kv0;
1411 }
1412 if (1 == n)
1413 {
1414 resultFor_vPlusOne = Kv0 + (2.0*(1.0 + v0)/x)*Kv0plus1;
1415 if (Kv0plus1 >= DBL_BIG)
1416 {
1417 // Try to avoid returning "inf" in either variable.
1418 resultFor_vPlusOne = DBL_BIG;
1419 return DBL_BIG;
1420 }
1421 if (resultFor_vPlusOne >= DBL_BIG)
1422 resultFor_vPlusOne = DBL_BIG;
1423 return Kv0plus1;
1424 }
1425 }
1426
1427 else // the size of v0 makes K(v0, x) = (PI/2)*(I(-v0, x) - I(v0, x))/sin(v0*PI) stable
1428 {
1429 // Variables for computing K(v0,x).
1430 double gamma_1_plus_v0 = my_gamma(1.0+v0),
1431 gamma_1_minus_v0 = my_gamma(1.0-v0);
1432 double factor_neg = 1.0/(gamma_1_minus_v0*pow(0.5*x, v0));
1433 double factor_pos = 1.0/(factor_neg*gamma_1_minus_v0*gamma_1_plus_v0);
1434 double neg_term_i = 1.0, pos_term_i = 1.0;
1435 double sum = factor_neg - factor_pos; // i=0 terms
1436 double prev_sum = DBL_BIG;
1437 // Variables for computing K(v0+1,x).
1438 double factor_neg_2 = -v0*factor_neg/(0.5*x);
1439 double factor_pos_2 = factor_pos*0.5*x/(1.0+v0);
1440 double neg_term_i_2 = 1.0, pos_term_i_2 = 1.0;
1441 double sum_2 = factor_neg_2 - factor_pos_2; // i=0 terms
1442 double prev_sum_2 = DBL_BIG;
1443
1444 for (unsigned long int i = 1; i < 2.4*MAX_ITERATIONS; i++) // Zhang & Jin's heuristic cutoff
1445 {
1446 // Compute the sum for K(v0,x).
1447 neg_term_i *= half_x_squared/(i*(i - v0));
1448 pos_term_i *= half_x_squared/(i*(i + v0));
1449 sum += factor_neg*neg_term_i - factor_pos*pos_term_i;
1450
1451 // Compute the sum for K(v0+1,x).
1452 neg_term_i_2 *= half_x_squared/(i*(i - (v0+1.0)));
1453 pos_term_i_2 *= half_x_squared/(i*(i + (v0+1.0)));
1454 sum_2 += factor_neg_2*neg_term_i_2 - factor_pos_2*pos_term_i_2;
1455
1456 if (fabs((prev_sum - sum)/sum) < EPSILON_ZJ &&
1457 fabs((prev_sum_2 - sum_2)/sum_2) < EPSILON_ZJ)
1458 break; // convergence achieved
1459 prev_sum = sum;
1460 prev_sum_2 = sum_2;
1461 }
1462
1463 Kv0 = (0.5*PI)*sum/sin(v0*PI);
1464 Kv0plus1 = (0.5*PI)*sum_2/sin((v0+1.0)*PI);
1465 }
1466 }
1467 }
1468 else // x > 9.
1469 {
1470 // Large argument: Use the asymptotic expansion.
1471 //
1472 // K(v0,x) = sqrt(PI/(2x))*exp(-x)*(1 + sum(product[(u - (2k-1)^2)/(8k*x),
1473 // from u=1 to u=m],
1474 // from m=1 to m=kmax)),
1475 // where u = 4*v0*v0
1476 // and kmax is an integer determined heuristically by Zhang & Jin (1996).
1477 //
1478
1479 double sum(0.0), product(1.0), sum_2(0.0), product_2(1.0);
1480 double four_v0_squared = 4.0*v0*v0; // avoid repeated recalculations
1481 double four_v0plus1_squared = 4.0*(v0+1.0)*(v0+1.0); // ditto
1482 unsigned long int kmax = 14; // kmax is Zhang & Jin's heuristic cutoff
1483
1484 if (x >= 50.0)
1485 kmax = 8; // Zhang & Jin heuristics
1486 else if (x >= 35.0)
1487 kmax = 10; // Zhang & Jin heuristics
1488
1489 for (unsigned long int k = 1; k <= kmax; k++)
1490 {
1491 product *= (four_v0_squared - (2.0*k - 1.0)*(2.0*k - 1.0))/(8.0*k*x);
1492 sum += product;
1493 product_2 *= (four_v0plus1_squared - (2.0*k - 1.0)*(2.0*k - 1.0))/(8.0*k*x);
1494 sum_2 += product_2;
1495 }
1496 Kv0 = sqrt(PI/(2.0*x))*exp(-x)*(1.0 + sum);
1497 Kv0plus1 = sqrt(PI/(2.0*x))*exp(-x)*(1.0 + sum_2);
1498 }
1499
1500 if (mustUseTwo_v0s)
1501 {
1502 if (first_v0_calculation)
1503 {
1504 result = Kv0;
1505 first_v0_calculation = false;
1506 v0 = -v + 1.0;
1507 goto Calculate_Kv0;
1508 }
1509 else
1510 resultFor_vPlusOne = Kv0;
1511 if (result >= DBL_BIG)
1512 {
1513 // Try to avoid returning "inf" in either variable.
1514 result = resultFor_vPlusOne = DBL_BIG;
1515 }
1516 else if (resultFor_vPlusOne >= DBL_BIG)
1517 resultFor_vPlusOne = DBL_BIG;
1518 return result;
1519 }
1520
1521 // Now, use K(v0,x) and K(v0+1,x) to calculate the desired K(v,x) == K(v0+n,x),
1522 // using the recurrence relationship:
1523 //
1524 // K(v+1,x) = K(v-1,x) + (2v/x)K(v,x).
1525
1526 double A = 0.0, B = Kv0, C = Kv0plus1;
1527 for (unsigned long int i = 1; i <= n; i++)
1528 {
1529 A = B;
1530 B = C;
1531 C = A + (2.0*(v0 + i)/x)*B;
1532 }
1533
1534 if (orderIsNegative)
1535 swap(B,C);
1536 if (B >= DBL_BIG)
1537 B = C = DBL_BIG;
1538 else if (C >= DBL_BIG)
1539 C = DBL_BIG;
1540 resultFor_vPlusOne = C;
1541 return B;
1542 }
1543
1544 // Function to compute a numerical estimate of the analytically
1545 // intractable partial derivative of K(v,x) with repect to v,
1546 // where K(v,x) is the modified Bessel function of the second kind
1547 // of order v evaluated at x, sometimes called the MacDonald function.
1548 // [Note: The partial derivative of K(v,x) with repect to x
1549 // is simply (v/x)*K(v,x) - K(v+1,x).]
1550 //
1551 // Arguments: v is the order, x is the argument,
1552 // Kvx is the already-computed value of K(v,x).
1553 //
1554 // WARNING: This simple algorithm may be insufficiently precise!
1555 //
1556 // Contains an original approximation derived by erynes for when
1557 // v is close to an integer.
1558 //
DvBesselK(double v,double x,double Kvx)1559 double DvBesselK(double v, double x, double Kvx)
1560 {
1561 if (x <= 0.0)
1562 throw implementation_error("DvBesselK() received illegal argument: "
1563 + ToString(x) + ".");
1564 if (0.0 == v)
1565 return 0.0; // this is independent of x, for nonzero x
1566
1567 // K(-v,x) == K(v,x), so the sign of the value we return
1568 // must equal the sign of v. Work with v > 0 for convenience.
1569 bool vIsNegative(false);
1570 if (v < 0.0)
1571 {
1572 vIsNegative = true;
1573 v = -v;
1574 }
1575
1576 // For tiny x, K(v,x) = 0.5*gamma(v)*pow(0.5*x, -v).
1577 double gamma_v__DividedBy__DBL_BIG = my_gamma(v)/DBL_BIG;
1578 if (x <= 2.0*pow(0.5*gamma_v__DividedBy__DBL_BIG, 1.0/v))
1579 return (vIsNegative ? -DBL_BIG : DBL_BIG);
1580 else if (DBL_BIG == Kvx)
1581 return (vIsNegative ? -DBL_BIG : DBL_BIG);
1582
1583 unsigned long int n = static_cast<unsigned long int>(floor(v));
1584 double v0 = v - 1.0*n;
1585 const double h = 1.0e-06; // used to compute the derivative from its definition
1586 double dummy; // dummy variable needed for storage; value not used
1587
1588 if (v0 <= SMALL_v0 && v0+h >= SMALL_v0)
1589 {
1590 // need to nudge backwards
1591 v0 = (SMALL_v0 - NUDGE_AMOUNT) - h;
1592 v = n + v0;
1593 Kvx = BesselK(v,x,dummy);
1594 }
1595 else if (v0 < 1.0 - SMALL_v0 && v0+h >= 1.0 - SMALL_v0)
1596 {
1597 // need to nudge backwards
1598 v0 = (1.0 - SMALL_v0 - NUDGE_AMOUNT) - h;
1599 v = n + v0;
1600 Kvx = BesselK(v,x,dummy);
1601 }
1602 else if (v0 > 1.0 - SMALL_v0 && v0+h == 1.0)
1603 {
1604 // need to nudge backwards
1605 v0 = (1.0 - NUDGE_AMOUNT) - h;
1606 v = n + v0;
1607 Kvx = BesselK(v,x,dummy);
1608 }
1609 // Note that v0 > 1.0 - SMALL_v0 && v0+h > 1.0 is safe,
1610 // because h < SMALL_v0, hence v0+h < 1.0 + SMALL_v0.
1611 // Also is v0 == 1.0 - SMALL_v0 is safe, because in that case
1612 // K(v0,x) and K(v0+h,x) will both be calculated using the
1613 // derived approximation.
1614
1615 double K_v_plus_h__x(0.0), result;
1616
1617 K_v_plus_h__x = BesselK(v+h, x, dummy);
1618 result = (vIsNegative ? -(K_v_plus_h__x - Kvx)/h :
1619 (K_v_plus_h__x - Kvx)/h);
1620 if (result <= -DBL_BIG)
1621 result = -DBL_BIG;
1622 else if (result >= DBL_BIG)
1623 result = DBL_BIG;
1624 return result;
1625
1626 #if 0 // Potentially DEAD CODE (bobgian, Feb 2010)
1627 // This is old code by erynes that computes the derivative of the approximation
1628 // that's conditionally employed in BesselK(). For at least the moment (early
1629 // 2007), we're taking it out, and always computing the derivative using an
1630 // approximation to the definition of the derivative (i.e., slope computed
1631 // at points x and x+h where h is small).
1632
1633 if (v0 > 1.0 - SMALL_v0)
1634 {
1635 n += 1;
1636 v0 = v - 1.0*n; // Thus, v0 < 0.0 Example: 2.0 + 0.95 --> 3.0 - 0.05.
1637 }
1638
1639 double K1x(0.0);
1640 double K0x = BesselK(0.0, x, K1x);
1641 double A = 0.0;
1642 // erynes determined that, for small "e", e > 0 or e < 0,
1643 // K(e,x) = K(0,x)*(1 + e^2/(2x + 1)),
1644 // K(1+e,x) = K(1,x) + (K(0,x)/x)*e + ((2x-1)(4x+1)K(1,x)/(2*(x+1)^2))*e^2,
1645 // where the "smallness" of e can vary with x, but is often near 0.01.
1646 // "e" in these equations corresponds to "v0" in the code.
1647 double B = (2.0 * K0x/(2.0 * x + 1.0)) * v0;
1648 if (0 == n)
1649 return (vIsNegative ? -B : B);
1650 double C = K0x/x + ((2.0*x-1.0)*(4.0*x+1.0)*K1x/((x+1.0)*(x+1.0)))*v0; // recall x > 0
1651 Kvx = K1x + (K0x/x)*v0 + ((2.0*x-1.0)*(4.0*x+1.0)*K1x/(2.0*(x+1.0)*(x+1.0)))*v0*v0;
1652
1653 double Kv1x = K0x*(v0*v0/(2.0*x + 1.0) + 1.0) + (2.0*(1.0+v0)/x)*Kvx; // == K(2+v0,x)
1654
1655 for (unsigned long int i = 1; i < n; i++)
1656 {
1657 A = B;
1658 B = C;
1659 if (0 == i % 2)
1660 {
1661 C = A + (2.0/x)*(Kv1x + (i + v0)*B);
1662 if (i != n - 1) // avoid unnecessary computation
1663 Kvx = BesselK(1.0+i+v0, x, Kv1x);
1664 }
1665 else
1666 C = A + (2.0/x)*(Kvx + (i + v0)*B);
1667 }
1668
1669 return (vIsNegative ? -C : C);
1670 #endif // 0
1671 }
1672
1673 // Euler's psi function == d(log_gamma(x))/dx.
1674 // Implemented by erynes in August 2005, adapted from:
1675 // Shanjie Zhang and Jianming Jin, _Computation of Special Functions._
1676 // New York: Wiley-Interscience (1996).
1677 // (Their implementation is taken directly from Abramowitz and Stegun.)
1678 //
psi(double x)1679 double psi(double x)
1680 {
1681 double result(0.0);
1682 double abs_x = x >= 0.0 ? x : -x;
1683 if (floor(x) == x)
1684 {
1685 // x is an integer. Use the expansion for integers.
1686 if (x <= 0.0)
1687 {
1688 throw implementation_error("psi() received illegal argument: "
1689 + ToString(x) + ".");
1690 return DBL_BIG; // psi(x) is undefined for 0 and negative integers
1691 }
1692 for (unsigned long int k = 1; k < static_cast<unsigned long int>(floor(x)); k++)
1693 result += 1.0/k;
1694 return -EULERS_CONSTANT + result;
1695 }
1696
1697 if (floor(abs_x + 0.5) == abs_x + 0.5)
1698 {
1699 // x is a half-integer. Use the expansion for half-integers.
1700 for (unsigned long int k = 1;
1701 k <= static_cast<unsigned long int>(floor(abs_x - 0.5)); k++)
1702 result += 2.0/(2.0*k - 1.0);
1703 result += -EULERS_CONSTANT - 2.0*LOG2;
1704 }
1705 else
1706 {
1707 // Not an integer or half-integer.
1708 if (abs_x < 10.0)
1709 {
1710 // Add an integer to abs_x to make abs_x + n > 10,
1711 // then use the asymptotic expansion for large x,
1712 // then use the recurrence relationship between psi(x+n) and psi(x).
1713 // First, compute sum[1/(x+k)], then later, subtract this
1714 // from psi(x+n) to obtain psi(x).
1715 unsigned long int n = static_cast<unsigned long int>(10.0 - floor(abs_x));
1716 for (unsigned long int k = 0; k < n; k++)
1717 result -= 1.0/(abs_x + k);
1718 abs_x += n;
1719 }
1720 // The asymptotic expansion. These numbers are Bernoulli numbers;
1721 // the expansion actually comes from Abramowitz and Stegun's book.
1722 double x2 = 1.0/(abs_x*abs_x),
1723 a1 = -1.0/12.0,
1724 a2 = +1.0/120.0,
1725 a3 = -1.0/252.0,
1726 a4 = +1.0/240.0,
1727 a5 = -1.0/132.0,
1728 a6 = +691.0/32760.0,
1729 a7 = a1,
1730 a8 = +3617.0/8160.0;
1731 result += log(abs_x) - 0.5/abs_x +
1732 x2*(((((((a8*x2+a7)*x2+a6)*x2+a5)*x2+a4)*x2+a3)*x2+a2)*x2+a1);
1733 }
1734
1735 if (x < 0.0)
1736 result += 1.0/abs_x + PI/tan(PI*abs_x); // psi(-x) = psi(x) + 1/x + PI*cot(PI*x)
1737 return result;
1738 }
1739
1740 //____________________________________________________________________________________
1741