forester/archive/RIO/others/hmmer/squid/sre_math.c

   1 /*****************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  *****************************************************************/
  10
  11 /* sre_math.c
  12  *
  13  * Portability for and extensions to C math library.
  14  * RCS $Id: sre_math.c,v 1.1.1.1 2005/03/22 08:34:32 cmzmasek Exp $
  15  */
  16
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include "squid.h"
  21
  22 #ifdef MEMDEBUG
  23 #include "dbmalloc.h"
  24 #endif
  25
  26 static int sre_reseed = 0;      /* TRUE to reinit sre_random() */
  27 static int sre_randseed = 666;  /* default seed for sre_random()   */
  28
  29 /* Function: ExponentialRandom()
  30  * Date:     SRE, Mon Sep  6 21:24:29 1999 [St. Louis]
  31  *
  32  * Purpose:  Pick an exponentially distributed random variable
  33  *           0 > x >= infinity
  34  *
  35  * Args:     (void)
  36  *
  37  * Returns:  x
  38  */
  39 float
  40 ExponentialRandom(void)
  41 {
  42   float x;
  43
  44   do x = sre_random(); while (x == 0.0);
  45   return -log(x);
  46 }
  47
  48 /* Function: Gaussrandom()
  49  *
  50  * Pick a Gaussian-distributed random variable
  51  * with some mean and standard deviation, and
  52  * return it.
  53  *
  54  * Based on RANLIB.c public domain implementation.
  55  * Thanks to the authors, Barry W. Brown and James Lovato,
  56  * University of Texas, M.D. Anderson Cancer Center, Houston TX.
  57  * Their implementation is from Ahrens and Dieter, "Extensions
  58  * of Forsythe's method for random sampling from the normal
  59  * distribution", Math. Comput. 27:927-937 (1973).
  60  *
  61  * Impenetrability of the code is to be blamed on its FORTRAN/f2c lineage.
  62  *
  63  */
  64 float
  65 Gaussrandom(float mean, float stddev)
  66 {
  67   static float a[32] = {
  68     0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
  69     0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
  70     1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
  71     1.862732,2.153875
  72   };
  73   static float d[31] = {
  74     0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
  75     0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
  76     0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
  77     0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
  78   };
  79   static float t[31] = {
  80     7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
  81     1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
  82     2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
  83     4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
  84     9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
  85   };
  86   static float h[31] = {
  87     3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
  88     4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
  89     4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
  90     5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
  91     8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
  92   };
  93   static long i;
  94   static float snorm,u,s,ustar,aa,w,y,tt;
  95
  96   u = sre_random();
  97   s = 0.0;
  98   if(u > 0.5) s = 1.0;
  99   u += (u-s);
 100   u = 32.0*u;
 101   i = (long) (u);
 102   if(i == 32) i = 31;
 103   if(i == 0) goto S100;
 104   /*
 105    * START CENTER
 106    */
 107   ustar = u-(float)i;
 108   aa = *(a+i-1);
 109 S40:
 110   if(ustar <= *(t+i-1)) goto S60;
 111   w = (ustar-*(t+i-1))**(h+i-1);
 112 S50:
 113   /*
 114    * EXIT   (BOTH CASES)
 115    */
 116   y = aa+w;
 117   snorm = y;
 118   if(s == 1.0) snorm = -y;
 119   return (stddev*snorm + mean);
 120 S60:
 121   /*
 122    * CENTER CONTINUED
 123    */
 124   u = sre_random();
 125   w = u*(*(a+i)-aa);
 126   tt = (0.5*w+aa)*w;
 127   goto S80;
 128 S70:
 129   tt = u;
 130   ustar = sre_random();
 131 S80:
 132   if(ustar > tt) goto S50;
 133   u = sre_random();
 134   if(ustar >= u) goto S70;
 135   ustar = sre_random();
 136   goto S40;
 137 S100:
 138   /*
 139    * START TAIL
 140    */
 141   i = 6;
 142   aa = *(a+31);
 143   goto S120;
 144 S110:
 145   aa += *(d+i-1);
 146   i += 1;
 147 S120:
 148   u += u;
 149   if(u < 1.0) goto S110;
 150   u -= 1.0;
 151 S140:
 152   w = u**(d+i-1);
 153   tt = (0.5*w+aa)*w;
 154   goto S160;
 155 S150:
 156   tt = u;
 157 S160:
 158   ustar = sre_random();
 159   if(ustar > tt) goto S50;
 160   u = sre_random();
 161   if(ustar >= u) goto S150;
 162   u = sre_random();
 163   goto S140;
 164 }
 165
 166
 167 /* Function: Linefit()
 168  *
 169  * Purpose:  Given points x[0..N-1] and y[0..N-1], fit to
 170  *           a straight line y = a + bx.
 171  *           a, b, and the linear correlation coefficient r
 172  *           are filled in for return.
 173  *
 174  * Args:     x     - x values of data
 175  *           y     - y values of data
 176  *           N     - number of data points
 177  *           ret_a - RETURN: intercept
 178  *           ret_b - RETURN: slope
 179  *           ret_r - RETURN: correlation coefficient
 180  *
 181  * Return:   1 on success, 0 on failure.
 182  */
 183 int
 184 Linefit(float *x, float *y, int N, float *ret_a, float *ret_b, float *ret_r)
 185 {
 186   float xavg, yavg;
 187   float sxx, syy, sxy;
 188   int   i;
 189
 190   /* Calculate averages, xavg and yavg
 191    */
 192   xavg = yavg = 0.0;
 193   for (i = 0; i < N; i++)
 194     {
 195       xavg += x[i];
 196       yavg += y[i];
 197     }
 198   xavg /= (float) N;
 199   yavg /= (float) N;
 200
 201   sxx = syy = sxy = 0.0;
 202   for (i = 0; i < N; i++)
 203     {
 204       sxx    += (x[i] - xavg) * (x[i] - xavg);
 205       syy    += (y[i] - yavg) * (y[i] - xavg);
 206       sxy    += (x[i] - xavg) * (y[i] - yavg);
 207     }
 208   *ret_b = sxy / sxx;
 209   *ret_a = yavg - xavg*(*ret_b);
 210   *ret_r = sxy / (sqrt(sxx) * sqrt(syy));
 211   return 1;
 212 }
 213
 214
 215 /* Function: WeightedLinefit()
 216  *
 217  * Purpose:  Given points x[0..N-1] and y[0..N-1] with
 218  *           variances (measurement errors) var[0..N-1],
 219  *           fit to a straight line y = mx + b.
 220  *
 221  * Method:   Algorithm from Numerical Recipes in C, [Press88].
 222  *
 223  * Return:   (void)
 224  *           ret_m contains slope; ret_b contains intercept
 225  */
 226 void
 227 WeightedLinefit(float *x, float *y, float *var, int N, float *ret_m, float *ret_b)
 228 {
 229   int    i;
 230   double s;
 231   double sx, sy;
 232   double sxx, sxy;
 233   double delta;
 234   double m, b;
 235
 236   s = sx = sy = sxx = sxy = 0.;
 237   for (i = 0; i < N; i++)
 238     {
 239       s   += 1./var[i];
 240       sx  += x[i] / var[i];
 241       sy  += y[i] / var[i];
 242       sxx += x[i] * x[i] / var[i];
 243       sxy += x[i] * y[i] / var[i];
 244     }
 245
 246   delta = s * sxx - (sx * sx);
 247   b = (sxx * sy - sx * sxy) / delta;
 248   m = (s * sxy - sx * sy) / delta;
 249
 250   *ret_m = m;
 251   *ret_b = b;
 252 }
 253
 254
 255 /* Function: Gammln()
 256  *
 257  * Returns the natural log of the gamma function of x.
 258  * x is > 0.0.
 259  *
 260  * Adapted from a public domain implementation in the
 261  * NCBI core math library. Thanks to John Spouge and
 262  * the NCBI. (According to the NCBI, that's Dr. John
 263  * "Gammas Galore" Spouge to you, pal.)
 264  */
 265 double
 266 Gammln(double x)
 267 {
 268   int i;
 269   double xx, tx;
 270   double tmp, value;
 271   static double cof[11] = {
 272     4.694580336184385e+04,
 273     -1.560605207784446e+05,
 274     2.065049568014106e+05,
 275     -1.388934775095388e+05,
 276     5.031796415085709e+04,
 277     -9.601592329182778e+03,
 278     8.785855930895250e+02,
 279     -3.155153906098611e+01,
 280     2.908143421162229e-01,
 281     -2.319827630494973e-04,
 282     1.251639670050933e-10
 283   };
 284
 285   /* Protect against x=0. We see this in Dirichlet code,
 286    * for terms alpha = 0. This is a severe hack but it is effective
 287    * and (we think?) safe. (due to GJM)
 288    */
 289   if (x <= 0.0) return 999999.;
 290
 291   xx       = x - 1.0;
 292   tx = tmp = xx + 11.0;
 293   value    = 1.0;
 294   for (i = 10; i >= 0; i--)     /* sum least significant terms first */
 295     {
 296       value += cof[i] / tmp;
 297       tmp   -= 1.0;
 298     }
 299   value  = log(value);
 300   tx    += 0.5;
 301   value += 0.918938533 + (xx+0.5)*log(tx) - tx;
 302   return value;
 303 }
 304
 305
 306 /* Vector operations for doubles and floats.
 307  * DNorm(), FNorm()   -- normalize a probability vector of length n.
 308  *                       return 0 if all values were zero.
 309  * DScale(), FScale() -- multiply all items in vector by scale
 310  * DSet(), FSet()     -- set all items in vector to value.
 311  * DAdd(), FAdd()     -- add vec2 to vec1.
 312  * DDot(), FDot()     -- calculate dot product of two vectors.
 313  * DCopy(), FCopy()   -- set vec1 to be same as vec2.
 314  * DMax(), FMax()     -- return index of maximum element in vec
 315  */
 316 int
 317 DNorm(double *vec, int n)
 318 {
 319   int    x;
 320   double sum;
 321
 322   sum = 0.0;
 323   for (x = 0; x < n; x++) sum += vec[x];
 324   if (sum != 0.0)
 325     for (x = 0; x < n; x++) vec[x] /= sum;
 326   else
 327     { squid_errno = SQERR_DIVZERO; return 0; }
 328   return 1;
 329 }
 330 int
 331 FNorm(float *vec, int n)
 332 {
 333   int    x;
 334   float  sum;
 335
 336   sum = 0.0;
 337   for (x = 0; x < n; x++) sum += vec[x];
 338   if (sum != 0.0)
 339     for (x = 0; x < n; x++) vec[x] /= sum;
 340   else
 341     { squid_errno = SQERR_DIVZERO; return 0; }
 342   return 1;
 343 }
 344
 345 void
 346 DScale(double *vec, int n, double scale)
 347 {
 348   int x;
 349   for (x = 0; x < n; x++)
 350     vec[x] *= scale;
 351 }
 352 void
 353 FScale(float *vec, int n, float scale)
 354 {
 355   int x;
 356   for (x = 0; x < n; x++)
 357     vec[x] *= scale;
 358 }
 359
 360 void
 361 DSet(double *vec, int n, double value)
 362 {
 363   int x;
 364   for (x = 0; x < n; x++)
 365     vec[x] = value;
 366 }
 367 void
 368 FSet(float *vec, int n, float value)
 369 {
 370   int x;
 371   for (x = 0; x < n; x++)
 372     vec[x] = value;
 373 }
 374
 375 double
 376 DSum(double *vec, int n)
 377 {
 378   double sum = 0.;
 379   int    x;
 380   for (x = 0; x < n; x++)
 381     sum += vec[x];
 382   return sum;
 383 }
 384 float
 385 FSum(float *vec, int n)
 386 {
 387   float sum = 0.;
 388   int   x;
 389   for (x = 0; x < n; x++)
 390     sum += vec[x];
 391   return sum;
 392 }
 393
 394 void
 395 DAdd(double *vec1, double *vec2, int n)
 396 {
 397   int x;
 398   for (x = 0; x < n; x++)
 399     vec1[x] += vec2[x];
 400 }
 401 void
 402 FAdd(float *vec1, float *vec2, int n)
 403 {
 404   int x;
 405   for (x = 0; x < n; x++)
 406     vec1[x] += vec2[x];
 407 }
 408
 409 void
 410 DCopy(double *vec1, double *vec2, int n)
 411 {
 412   int x;
 413   for (x = 0; x < n; x++)
 414     vec1[x] = vec2[x];
 415 }
 416 void
 417 FCopy(float *vec1, float *vec2, int n)
 418 {
 419   int x;
 420   for (x = 0; x < n; x++)
 421     vec1[x] = vec2[x];
 422 }
 423
 424 double
 425 DDot(double *vec1, double *vec2, int n)
 426 {
 427   double result = 0.;
 428   int x;
 429
 430   for (x = 0; x < n; x++)
 431     result += vec1[x] * vec2[x];
 432   return result;
 433 }
 434 float
 435 FDot(float *vec1, float *vec2, int n)
 436 {
 437   float result = 0.;
 438   int x;
 439
 440   for (x = 0; x < n; x++)
 441     result += vec1[x] * vec2[x];
 442   return result;
 443 }
 444
 445 /* Functions: DMax(), FMax()
 446  * Date:      SRE, Fri Aug 29 11:14:08 1997 (Denver CO)
 447  *
 448  * Purpose:   return index of maximum element in vec.
 449  */
 450 int
 451 DMax(double *vec, int n)
 452 {
 453   int i;
 454   int best = 0;
 455
 456   for (i = 1; i < n; i++)
 457     if (vec[i] > vec[best]) best = i;
 458   return best;
 459 }
 460 int
 461 FMax(float *vec, int n)
 462 {
 463   int i;
 464   int best = 0;
 465
 466   for (i = 1; i < n; i++)
 467     if (vec[i] > vec[best]) best = i;
 468   return best;
 469 }
 470
 471
 472 /* 2D matrix operations
 473  */
 474 float **
 475 FMX2Alloc(int rows, int cols)
 476 {
 477   float **mx;
 478   int     r;
 479
 480   mx    = (float **) MallocOrDie(sizeof(float *) * rows);
 481   mx[0] = (float *)  MallocOrDie(sizeof(float) * rows * cols);
 482   for (r = 1; r < rows; r++)
 483     mx[r] = mx[0] + r*cols;
 484   return mx;
 485 }
 486 void
 487 FMX2Free(float **mx)
 488 {
 489   free(mx[0]);
 490   free(mx);
 491 }
 492 double **
 493 DMX2Alloc(int rows, int cols)
 494 {
 495   double **mx;
 496   int      r;
 497
 498   mx    = (double **) MallocOrDie(sizeof(double *) * rows);
 499   mx[0] = (double *)  MallocOrDie(sizeof(double) * rows * cols);
 500   for (r = 1; r < rows; r++)
 501     mx[r] = mx[0] + r*cols;
 502   return mx;
 503 }
 504 void
 505 DMX2Free(double **mx)
 506 {
 507   free(mx[0]);
 508   free(mx);
 509 }
 510 /* Function: FMX2Multiply()
 511  *
 512  * Purpose:  Matrix multiplication.
 513  *           Multiply an m x p matrix A by a p x n matrix B,
 514  *           giving an m x n matrix C.
 515  *           Matrix C must be a preallocated matrix of the right
 516  *           size.
 517  */
 518 void
 519 FMX2Multiply(float **A, float **B, float **C, int m, int p, int n)
 520 {
 521   int i, j, k;
 522
 523   for (i = 0; i < m; i++)
 524     for (j = 0; j < n; j++)
 525       {
 526         C[i][j] = 0.;
 527         for (k = 0; k < p; k++)
 528           C[i][j] += A[i][p] * B[p][j];
 529       }
 530 }
 531
 532 /* Function: sre_random()
 533  *
 534  * Purpose:  Return a uniform deviate from 0.0 to 1.0.
 535  *           sre_randseed is a static variable, set
 536  *           by sre_srandom(). sre_reseed is a static flag
 537  *           raised by sre_srandom(), saying that we need
 538  *           to re-initialize.
 539  *           [0.0 <= x < 1.0]
 540  *
 541  *           Uses a simple linear congruential generator with
 542  *           period 2^28. Based on discussion in Robert Sedgewick's
 543  *           _Algorithms in C_, Addison-Wesley, 1990.
 544  *
 545  *           Requires that long int's have at least 32 bits.
 546  *
 547  * Reliable and portable, but slow. Benchmarks on wol,
 548  * using IRIX cc and IRIX C library rand() and random():
 549  *     sre_random():    0.8 usec/call
 550  *     random():        0.3 usec/call
 551  *     rand():          0.3 usec/call
 552  */
 553 #define RANGE 268435456         /* 2^28        */
 554 #define DIV   16384             /* sqrt(RANGE) */
 555 #define MULT  72530821          /* my/Cathy's birthdays, x21, x even (Knuth)*/
 556 float
 557 sre_random(void)
 558 {
 559   static long  rnd;
 560   static int   firsttime = 1;
 561   long         high1, low1;
 562   long         high2, low2;
 563
 564   if (sre_reseed || firsttime)
 565     {
 566       sre_reseed = firsttime = 0;
 567       if (sre_randseed <= 0) sre_randseed = 666; /* seeds of zero break me */
 568       high1 = sre_randseed / DIV;  low1  = sre_randseed % DIV;
 569       high2 = MULT / DIV;          low2  = MULT % DIV;
 570       rnd = (((high2*low1 + high1*low2) % DIV)*DIV + low1*low2) % RANGE;
 571     }
 572    high1 = rnd / DIV;  low1  = rnd % DIV;
 573    high2 = MULT / DIV; low2  = MULT % DIV;
 574    rnd = (((high2*low1 + high1*low2) % DIV)*DIV + low1*low2) % RANGE;
 575
 576    return ((float) rnd / (float) RANGE);
 577 }
 578 #undef RANGE
 579 #undef DIV
 580 #undef MULT
 581
 582
 583 /* Function: sre_srandom()
 584  *
 585  * Purpose:  Initialize with a random seed. Seed can be
 586  *           any integer.
 587  */
 588 void
 589 sre_srandom(int seed)
 590 {
 591   if (seed < 0) seed = -1 * seed;
 592   sre_reseed   = 1;
 593   sre_randseed = seed;
 594 }
 595
 596
 597 /* Functions: DChoose(), FChoose()
 598  *
 599  * Purpose:   Make a random choice from a normalized distribution.
 600  *            DChoose() is for double-precision vectors;
 601  *            FChoose() is for single-precision float vectors.
 602  *            Returns the number of the choice.
 603  */
 604 int
 605 DChoose(double *p, int N)
 606 {
 607   double roll;                  /* random fraction */
 608   double sum;                   /* integrated prob */
 609   int    i;                     /* counter over the probs */
 610
 611   roll    = sre_random();
 612   sum     = 0.0;
 613   for (i = 0; i < N; i++)
 614     {
 615       sum += p[i];
 616       if (roll < sum) return i;
 617     }
 618   SQD_DASSERT2((fabs(1.0 - sum) < 1e-14)); /* a verification at level 2 */
 619   return (int) (sre_random() * N);         /* bulletproof */
 620 }
 621 int
 622 FChoose(float *p, int N)
 623 {
 624   float roll;                   /* random fraction */
 625   float sum;                    /* integrated prob */
 626   int   i;                      /* counter over the probs */
 627
 628   roll    = sre_random();
 629   sum     = 0.0;
 630   for (i = 0; i < N; i++)
 631     {
 632       sum += p[i];
 633       if (roll < sum) return i;
 634     }
 635   SQD_DASSERT2((fabs(1.0f - sum) < 1e-6f));  /* a verification at level 2 */
 636   return (int) (sre_random() * N);           /* bulletproof */
 637 }
 638
 639 /* Functions: DLogSum(), FLogSum()
 640  *
 641  * Calculate the sum of a log vector
 642  * *in normal space*, and return the log of the sum.
 643  */
 644 double
 645 DLogSum(double *logp, int n)
 646 {
 647   int x;
 648   double max, sum;
 649
 650   max = logp[0];
 651   for (x = 1; x < n; x++)
 652     if (logp[x] > max) max = logp[x];
 653   sum = 0.0;
 654   for (x = 0; x < n; x++)
 655     if (logp[x] > max - 50.)
 656       sum += exp(logp[x] - max);
 657   sum = log(sum) + max;
 658   return sum;
 659 }
 660 float
 661 FLogSum(float *logp, int n)
 662 {
 663   int x;
 664   float max, sum;
 665
 666   max = logp[0];
 667   for (x = 1; x < n; x++)
 668     if (logp[x] > max) max = logp[x];
 669   sum = 0.0;
 670   for (x = 0; x < n; x++)
 671     if (logp[x] > max - 50.)
 672       sum += exp(logp[x] - max);
 673   sum = log(sum) + max;
 674   return sum;
 675 }
 676
 677
 678 /* Function: IncompleteGamma()
 679  *
 680  * Purpose:  Returns 1 - P(a,x) where:
 681  *           P(a,x) = \frac{1}{\Gamma(a)} \int_{0}^{x} t^{a-1} e^{-t} dt
 682  *                  = \frac{\gamma(a,x)}{\Gamma(a)}
 683  *                  = 1 - \frac{\Gamma(a,x)}{\Gamma(a)}
 684  *
 685  *           Used in a chi-squared test: for a X^2 statistic x
 686  *           with v degrees of freedom, call:
 687  *                  p = IncompleteGamma(v/2., x/2.)
 688  *           to get the probability p that a chi-squared value
 689  *           greater than x could be obtained by chance even for
 690  *           a correct model. (i.e. p should be large, say
 691  *           0.95 or more).
 692  *
 693  * Method:   Based on ideas from Numerical Recipes in C, Press et al.,
 694  *           Cambridge University Press, 1988.
 695  *
 696  * Args:     a  - for instance, degrees of freedom / 2     [a > 0]
 697  *           x  - for instance, chi-squared statistic / 2  [x >= 0]
 698  *
 699  * Return:   1 - P(a,x).
 700  */
 701 double
 702 IncompleteGamma(double a, double x)
 703 {
 704   int iter;                     /* iteration counter */
 705
 706   if (a <= 0.) Die("IncompleteGamma(): a must be > 0");
 707   if (x <  0.) Die("IncompleteGamma(): x must be >= 0");
 708
 709   /* For x > a + 1 the following gives rapid convergence;
 710    * calculate 1 - P(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)}:
 711    *     use a continued fraction development for \Gamma(a,x).
 712    */
 713   if (x > a+1)
 714     {
 715       double oldp;              /* previous value of p    */
 716       double nu0, nu1;          /* numerators for continued fraction calc   */
 717       double de0, de1;          /* denominators for continued fraction calc */
 718
 719       nu0 = 0.;                 /* A_0 = 0       */
 720       de0 = 1.;                 /* B_0 = 1       */
 721       nu1 = 1.;                 /* A_1 = 1       */
 722       de1 = x;                  /* B_1 = x       */
 723
 724       oldp = nu1;
 725       for (iter = 1; iter < 100; iter++)
 726         {
 727           /* Continued fraction development:
 728            * set A_j = b_j A_j-1 + a_j A_j-2
 729            *     B_j = b_j B_j-1 + a_j B_j-2
 730            * We start with A_2, B_2.
 731            */
 732                                 /* j = even: a_j = iter-a, b_j = 1 */
 733                                 /* A,B_j-2 are in nu0, de0; A,B_j-1 are in nu1,de1 */
 734           nu0 = nu1 + ((double)iter - a) * nu0;
 735           de0 = de1 + ((double)iter - a) * de0;
 736
 737                                 /* j = odd: a_j = iter, b_j = x */
 738                                 /* A,B_j-2 are in nu1, de1; A,B_j-1 in nu0,de0 */
 739           nu1 = x * nu0 + (double) iter * nu1;
 740           de1 = x * de0 + (double) iter * de1;
 741
 742                                 /* rescale */
 743           if (de1)
 744             {
 745               nu0 /= de1;
 746               de0 /= de1;
 747               nu1 /= de1;
 748               de1 =  1.;
 749             }
 750                                 /* check for convergence */
 751           if (fabs((nu1-oldp)/nu1) < 1.e-7)
 752             return nu1 * exp(a * log(x) - x - Gammln(a));
 753
 754           oldp = nu1;
 755         }
 756       Die("IncompleteGamma(): failed to converge using continued fraction approx");
 757     }
 758   else /* x <= a+1 */
 759     {
 760       double p;                 /* current sum               */
 761       double val;               /* current value used in sum */
 762
 763       /* For x <= a+1 we use a convergent series instead:
 764        *   P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)},
 765        * where
 766        *   \gamma(a,x) = e^{-x}x^a \sum_{n=0}{\infty} \frac{\Gamma{a}}{\Gamma{a+1+n}} x^n
 767        * which looks appalling but the sum is in fact rearrangeable to
 768        * a simple series without the \Gamma functions:
 769        *   = \frac{1}{a} + \frac{x}{a(a+1)} + \frac{x^2}{a(a+1)(a+2)} ...
 770        * and it's obvious that this should converge nicely for x <= a+1.
 771        */
 772
 773       p = val = 1. / a;
 774       for (iter = 1; iter < 10000; iter++)
 775         {
 776           val *= x / (a+(double)iter);
 777           p   += val;
 778
 779           if (fabs(val/p) < 1.e-7)
 780             return 1. - p * exp(a * log(x) - x - Gammln(a));
 781         }
 782       Die("IncompleteGamma(): failed to converge using series approx");
 783     }
 784   /*NOTREACHED*/
 785   return 0.;
 786 }
 787