single-gamic.c

     
   1  //! @file single-gamic.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! REAL generalised incomplete gamma function.
  25  
  26  // Generalised incomplete gamma code in this file was downloaded from 
  27  //   [http://helios.mi.parisdescartes.fr/~rabergel/]
  28  // and adapted for Algol 68 Genie.
  29  //
  30  // Reference:
  31  //   Rémy Abergel, Lionel Moisan. Fast and accurate evaluation of a
  32  //   generalized incomplete gamma function. 2019. hal-01329669v2
  33  //
  34  // Original source code copyright and license:
  35  //
  36  // DELTAGAMMAINC Fast and Accurate Evaluation of a Generalized Incomplete Gamma
  37  // Function. Copyright (C) 2016 Remy Abergel (remy.abergel AT gmail.com), Lionel
  38  // Moisan (Lionel.Moisan AT parisdescartes.fr).
  39  //
  40  // This file is a part of the DELTAGAMMAINC software, dedicated to the
  41  // computation of a generalized incomplete gammafunction. See the Companion paper
  42  // for a complete description of the algorithm.
  43  //
  44  // ``Fast and accurate evaluation of a generalized incomplete gamma function''
  45  // (Rémy Abergel, Lionel Moisan), preprint MAP5 nº2016-14, revision 1.
  46  //
  47  // This program is free software: you can redistribute it and/or modify it under
  48  // the terms of the GNU General Public License as published by the Free Software
  49  // Foundation, either version 3 of the License, or (at your option) any later
  50  // version.
  51  //
  52  // This program is distributed in the hope that it will be useful, but WITHOUT
  53  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  54  // FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  55  // details.
  56  //
  57  // You should have received a copy of the GNU General Public License along with
  58  // this program.  If not, see [http://www.gnu.org/licenses/].
  59  
  60  // References
  61  //
  62  //   R. Abergel and L. Moisan. 2016. Fast and accurate evaluation of a
  63  //   generalized incomplete gamma function, preprint MAP5 nº2016-14, revision 1
  64  //
  65  //   Rémy Abergel, Lionel Moisan. Fast and accurate evaluation of a
  66  //   generalized incomplete gamma function. 2019. hal-01329669v2
  67  //
  68  //   F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark
  69  //   (Eds.). 2010. NIST Handbook of Mathematical Functions. Cambridge University
  70  //   Press. (see online version at [http://dlmf.nist.gov/])
  71  //
  72  //   W. H. Press, S. A. Teukolsky, W. T. Vetterling, and
  73  //   B. P. Flannery. 1992. Numerical recipes in C: the art of scientific
  74  //   computing (2nd ed.).
  75  //
  76  //   G. R. Pugh, 2004. An analysis of the Lanczos Gamma approximation (phd
  77  //   thesis)
  78  
  79  #include "a68g.h"
  80  #include "a68g-genie.h"
  81  #include "a68g-prelude.h"
  82  #include "a68g-lib.h"
  83  #include "a68g-double.h"
  84  #include "a68g-mp.h"
  85  
  86  #define ITMAX 1000000000        // Maximum allowed number of iterations
  87  #define DPMIN DBL_MIN           // Number near the smallest representable double-point number
  88  #define EPS DBL_EPSILON         // Machine epsilon
  89  #define NITERMAX_ROMBERG 15     // Maximum allowed number of Romberg iterations
  90  #define TOL_ROMBERG 0.1         // Tolerance factor used to stop the Romberg iterations
  91  #define TOL_DIFF 0.2            // Tolerance factor used for the approximation of I_{x,y}^{mu,p} using differences
  92  
  93  // plim: compute plim (x), the limit of the partition of the domain (p,x)
  94  // detailed in the paper.
  95  //
  96  //            |      x              if   0 < x
  97  //            |
  98  // plim (x) = <      0              if -9 <= x <= 0
  99  //            |
 100  //            | 5.*sqrt (|x|)-5.    otherwise
 101  
 102  REAL_T plim (REAL_T x)
 103  {
 104    return (x >= 0) ? x : ((x >= -9) ? 0 : 5 * sqrt (-x) - 5);
 105  }
 106  
 107  //! @brief compute G(p,x) in the domain x <= p using a continued fraction
 108  //
 109  // p >= 0
 110  // x <= p
 111  
 112  void G_cfrac_lower (REAL_T * Gcfrac, REAL_T p, REAL_T x)
 113  {
 114    REAL_T c, d, del, f, an, bn;
 115    INT_T k, n;
 116  // deal with special case
 117    if (x == 0) {
 118      *Gcfrac = 0;
 119      return;
 120    }
 121  // Evaluate the continued fraction using Modified Lentz's method. However,
 122  // as detailed in the paper, perform manually the first pass (n=1), of the
 123  // initial Modified Lentz's method.
 124    an = 1;
 125    bn = p;
 126    f = an / bn;
 127    c = an / DPMIN;
 128    d = 1 / bn;
 129    n = 2;
 130    do {
 131      k = n / 2;
 132      an = (n & 1 ? k : -(p - 1 + k)) * x;
 133      bn++;
 134      d = an * d + bn;
 135      if (d == 0) {
 136        d = DPMIN;
 137      }
 138      c = bn + an / c;
 139      if (c == 0) {
 140        c = DPMIN;
 141      }
 142      d = 1 / d;
 143      del = d * c;
 144      f *= del;
 145      n++;
 146    }
 147    while ((a68_abs (del - 1.0) >= EPS) && (n < ITMAX));
 148    *Gcfrac = f;
 149  }
 150  
 151  //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
 152  // using a recursive integration by parts relation.
 153  // This function cannot be used when mu > 0.
 154  //
 155  // p > 0, integer
 156  // x < 0, |x| < max (1,p-1)
 157  
 158  void G_ibp (REAL_T * Gibp, REAL_T p, REAL_T x)
 159  {
 160    REAL_T t, tt, c, d, s, del;
 161    INT_T l;
 162    BOOL_T odd, stop;
 163    t = a68_abs (x);
 164    tt = 1 / (t * t);
 165    odd = (INT_T) a68_int (p) % 2 != 0;
 166    c = 1 / t;
 167    d = (p - 1);
 168    s = c * (t - d);
 169    l = 0;
 170    do {
 171      c *= d * (d - 1) * tt;
 172      d -= 2;
 173      del = c * (t - d);
 174      s += del;
 175      l++;
 176      stop = a68_abs (del) < a68_abs (s) * EPS;
 177    }
 178    while ((l < floor ((p - 2) / 2)) && !stop);
 179    if (odd && !stop) {
 180      s += d * c / t;
 181    }
 182    *Gibp = ((odd ? -1 : 1) * a68_exp (-t + lgamma (p) - (p - 1) * a68_ln (t)) + s) / t;
 183  }
 184  
 185  //! @brief compute the G-function in the domain x > p using a
 186  // continued fraction.
 187  //
 188  // p > 0
 189  // x > p, or x = +infinity
 190  
 191  void G_cfrac_upper (REAL_T * Gcfrac, REAL_T p, REAL_T x)
 192  {
 193  // Special case
 194    if (a68_isinf (x)) {
 195      *Gcfrac = 0;
 196      return;
 197    }
 198  // Evaluate the continued fraction using Modified Lentz's method. However,
 199  // as detailed in the paper, perform manually the first pass (n=1), of the
 200  // initial Modified Lentz's method.
 201    REAL_T an = 1, bn = x + 1 - p;
 202    REAL_T c, d, del, f;
 203    BOOL_T t = (bn != 0);
 204    INT_T i, n;
 205    if (t) {
 206  // b{1} is non-zero
 207      f = an / bn;
 208      c = an / DPMIN;
 209      d = 1 / bn;
 210      n = 2;
 211    } else {
 212  // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
 213      an = -(1 - p);
 214      bn = x + 3 - p;
 215      f = an / bn;
 216      c = an / DPMIN;
 217      d = 1 / bn;
 218      n = 3;
 219    }
 220    i = n - 1;
 221    do {
 222      an = -i * (i - p);
 223      bn += 2;
 224      d = an * d + bn;
 225      if (d == 0) {
 226        d = DPMIN;
 227      }
 228      c = bn + an / c;
 229      if (c == 0) {
 230        c = DPMIN;
 231      }
 232      d = 1 / d;
 233      del = d * c;
 234      f *= del;
 235      i++;
 236      n++;
 237    }
 238    while ((a68_abs (del - 1.0) >= EPS) && (n < ITMAX));
 239    *Gcfrac = t ? f : 1 / f;
 240  }
 241  
 242  //! @brief compute G : (p,x) --> R defined as follows
 243  //
 244  // if x <= p:
 245  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
 246  // otherwise:
 247  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
 248  //
 249  //   p > 0
 250  //   x is a real number or +infinity.
 251  
 252  void G_func (REAL_T * G, REAL_T p, REAL_T x)
 253  {
 254    if (p >= plim (x)) {
 255      G_cfrac_lower (G, p, x);
 256    } else if (x < 0) {
 257      G_ibp (G, p, x);
 258    } else {
 259      G_cfrac_upper (G, p, x);
 260    }
 261  }
 262  
 263  //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
 264  
 265  void romberg_iterations (REAL_T * R, REAL_T sigma, INT_T n, REAL_T x, REAL_T y, REAL_T mu, REAL_T p, REAL_T h, REAL_T pow2)
 266  {
 267    INT_T j, m;
 268    REAL_T sum, xx;
 269    INT_T adr0_prev = ((n - 1) * n) / 2;
 270    INT_T adr0 = (n * (n + 1)) / 2;
 271    for (sum = 0, j = 1; j <= pow2; j++) {
 272      xx = x + ((y - x) * (2 * j - 1)) / (2 * pow2);
 273      sum += a68_exp (-mu * xx + (p - 1) * a68_ln (xx) - sigma);
 274    }
 275    R[adr0] = 0.5 * R[adr0_prev] + h * sum;
 276    REAL_T pow4 = 4;
 277    for (m = 1; m <= n; m++) {
 278      R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1);
 279      pow4 *= 4;
 280    }
 281  }
 282  
 283  //! @ compute I_{x,y}^{mu,p} using a Romberg approximation.
 284  // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
 285  
 286  void romberg_estimate (REAL_T * rho, REAL_T * sigma, REAL_T x, REAL_T y, REAL_T mu, REAL_T p)
 287  {
 288    REAL_T *R = (REAL_T *) get_heap_space (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2 * sizeof (REAL_T));
 289    ASSERT (R != NULL);
 290  // Initialization (n=1)
 291    *sigma = -mu * y + (p - 1) * a68_ln (y);
 292    R[0] = 0.5 * (y - x) * (a68_exp (-mu * x + (p - 1) * a68_ln (x) - (*sigma)) + 1);
 293  // Loop for n > 0
 294    REAL_T relneeded = EPS / TOL_ROMBERG;
 295    INT_T adr0 = 0;
 296    INT_T n = 1;
 297    REAL_T h = (y - x) / 2;       // n=1, h = (y-x)/2^n
 298    REAL_T pow2 = 1;              // n=1; pow2 = 2^(n-1)
 299    if (NITERMAX_ROMBERG >= 1) {
 300      REAL_T relerr;
 301      do {
 302        romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
 303        h /= 2;
 304        pow2 *= 2;
 305        adr0 = (n * (n + 1)) / 2;
 306        relerr = a68_abs ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
 307        n++;
 308      } while (n <= NITERMAX_ROMBERG && relerr > relneeded);
 309    }
 310  // save Romberg estimate and free memory
 311    *rho = R[adr0 + (n - 1)];
 312    a68_free (R);
 313  }
 314  
 315  //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
 316  //
 317  //   I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
 318  //
 319  // This procedure computes (rho, sigma) described below.
 320  // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
 321  //
 322  //   mu is a real number non equal to zero 
 323  //     (in general we take mu = 1 or -1 but any nonzero real number is allowed)
 324  //
 325  //   x, y are two numbers with 0 <= x <= y <= +infinity,
 326  //     (the setting y=+infinity is allowed only when mu > 0)
 327  //
 328  //   p is a real number > 0, p must be an integer when mu < 0.
 329  
 330  void deltagammainc (REAL_T * rho, REAL_T * sigma, REAL_T x, REAL_T y, REAL_T mu, REAL_T p)
 331  {
 332  // Particular cases
 333    if (a68_isinf (x) && a68_isinf (y)) {
 334      *rho = 0;
 335      *sigma = a68_neginf ();
 336      return;
 337    } else if (x == y) {
 338      *rho = 0;
 339      *sigma = a68_neginf ();
 340      return;
 341    }
 342    if (x == 0 && a68_isinf (y)) {
 343      *rho = 1;
 344      (*sigma) = lgamma (p) - p * a68_ln (mu);
 345      return;
 346    }
 347  // Initialization
 348    REAL_T mA, mB, mx, my, nA, nB, nx, ny;
 349    G_func (&mx, p, mu * x);
 350    nx = (a68_isinf (x) ? a68_neginf () : -mu * x + p * a68_ln (x));
 351    G_func (&my, p, mu * y);
 352    ny = (a68_isinf (y) ? a68_neginf () : -mu * y + p * a68_ln (y));
 353  // Compute (mA,nA) and (mB,nB) such as I_{x,y}^{mu,p} can be
 354  // approximated by the difference A-B, where A >= B >= 0, A = mA*exp (nA) an 
 355  // B = mB*exp (nB). When the difference involves more than one digit loss due to
 356  // cancellation errors, the integral I_{x,y}^{mu,p} is evaluated using the
 357  // Romberg approximation method.
 358  
 359    if (mu < 0) {
 360      mA = my;
 361      nA = ny;
 362      mB = mx;
 363      nB = nx;
 364    } else {
 365      if (p < plim (mu * x)) {
 366        mA = mx;
 367        nA = nx;
 368        mB = my;
 369        nB = ny;
 370      } else if (p < plim (mu * y)) {
 371        mA = 1;
 372        nA = lgamma (p) - p * a68_ln (mu);
 373        nB = fmax (nx, ny);
 374        mB = mx * a68_exp (nx - nB) + my * a68_exp (ny - nB);
 375      } else {
 376        mA = my;
 377        nA = ny;
 378        mB = mx;
 379        nB = nx;
 380      }
 381    }
 382  // Compute (rho,sigma) such that rho*exp (sigma) = A-B
 383    *rho = mA - mB * a68_exp (nB - nA);
 384    *sigma = nA;
 385  // If the difference involved a significant loss of precision, compute Romberg estimate.
 386    if (!a68_isinf (y) && ((*rho) / mA < TOL_DIFF)) {
 387      romberg_estimate (rho, sigma, x, y, mu, p);
 388    }
 389  }
 390  
 391  // A68G Driver routines
 392  
 393  //! @brief PROC gamma inc g = (REAL p, x, y, mu) REAL
 394  
 395  void genie_gamma_inc_g_real (NODE_T * n)
 396  {
 397    A68_REAL x, y, mu, p;
 398    POP_OBJECT (n, &mu, A68_REAL);
 399    POP_OBJECT (n, &y, A68_REAL);
 400    POP_OBJECT (n, &x, A68_REAL);
 401    POP_OBJECT (n, &p, A68_REAL);
 402    REAL_T rho, sigma;
 403    deltagammainc (&rho, &sigma, VALUE (&x), VALUE (&y), VALUE (&mu), VALUE (&p));
 404    PUSH_VALUE (n, rho * a68_exp (sigma), A68_REAL);
 405  }
 406  
 407  //! @brief PROC gamma inc f = (REAL p, x) REAL
 408  
 409  void genie_gamma_inc_f_real (NODE_T * n)
 410  {
 411    A68_REAL x, p;
 412    POP_OBJECT (n, &x, A68_REAL);
 413    POP_OBJECT (n, &p, A68_REAL);
 414    REAL_T rho, sigma;
 415    deltagammainc (&rho, &sigma, VALUE (&x), a68_posinf (), 1, VALUE (&p));
 416    PUSH_VALUE (n, rho * a68_exp (sigma), A68_REAL);
 417  }
 418  
 419  //! @brief PROC gamma inc = (REAL p, x) REAL
 420  
 421  void genie_gamma_inc_h_real (NODE_T * n)
 422  {
 423  #if (A68_LEVEL >= 3) && defined (HAVE_GNU_MPFR)
 424    genie_gamma_inc_real_mpfr (n);
 425  #else
 426    genie_gamma_inc_f_real (n);
 427  #endif
 428  }
 429  
 430  //! @brief PROC gamma inc gf = (REAL p, x) REAL
 431  
 432  void genie_gamma_inc_gf_real (NODE_T * q)
 433  {
 434  // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
 435  // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
 436    A68_REAL x, p;
 437    POP_OBJECT (q, &x, A68_REAL);
 438    POP_OBJECT (q, &p, A68_REAL);
 439    REAL_T G;
 440    G_func (&G, VALUE (&p), VALUE (&x));
 441    PUSH_VALUE (q, G, A68_REAL);
 442  }