double-gamic.c

     
   1  //! @file double-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  //! LONG 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  
  81  #if (A68_LEVEL >= 3)
  82  
  83  #include "a68g-double.h"
  84  #include "a68g-genie.h"
  85  #include "a68g-lib.h"
  86  #include "a68g-mp.h"
  87  #include "a68g-prelude.h"
  88  #include "a68g-quad.h"
  89  
  90  #define DPMIN FLT128_MIN        // Number near the smallest representable double-point number
  91  #define EPS FLT128_EPSILON      // Machine epsilon
  92  #define ITMAX 1000000000        // Maximum allowed number of iterations
  93  #define NITERMAX_ROMBERG 15     // Maximum allowed number of Romberg iterations
  94  #define TOL_DIFF 0.2q           // Tolerance factor used for the approximation of I_{x,y}^{mu,p} using differences
  95  #define TOL_ROMBERG 0.1q        // Tolerance factor used to stop the Romberg iterations
  96  
  97  // double_plim: compute plim (x), the limit of the partition of the domain (p,x)
  98  // detailed in the paper.
  99  //
 100  //            |      x              if   0 < x
 101  //            |
 102  // plim (x) = <      0              if -9 <= x <= 0
 103  //            |
 104  //            | 5.*sqrt (|x|)-5.    otherwise
 105  
 106  DOUBLE_T double_plim (DOUBLE_T x)
 107  {
 108    return (x >= 0.0q) ? x : ((x >= -9.0q) ? 0.0q : 5.0q * sqrt (-x) - 5.0q);
 109  }
 110  
 111  //! @brief compute G(p,x) in the domain x <= p using a continued fraction
 112  //
 113  // p >= 0
 114  // x <= p
 115  
 116  void double_G_cfrac_lower (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
 117  {
 118  // deal with special case
 119    if (x == 0.0q) {
 120      *Gcfrac = 0.0q;
 121      return;
 122    }
 123  // Evaluate the continued fraction using Modified Lentz's method. However,
 124  // as detailed in the paper, perform manually the first pass (n=1), of the
 125  // initial Modified Lentz's method.
 126    DOUBLE_T an = 1.0q, bn = p, del;
 127    DOUBLE_T f = an / bn, c = an / DPMIN, d = 1.0q / bn;
 128    INT_T n = 2;
 129    do {
 130      INT_T k = n / 2;
 131      an = (n & 1 ? k : -(p - 1.0q + k)) * x;
 132      bn++;
 133      d = an * d + bn;
 134      if (d == 0.0q) {
 135        d = DPMIN;
 136      }
 137      c = bn + an / c;
 138      if (c == 0.0q) {
 139        c = DPMIN;
 140      }
 141      d = 1.0q / d;
 142      del = d * c;
 143      f *= del;
 144      n++;
 145    }
 146    while ((fabsq (del - 1.0q) >= EPS) && (n < ITMAX));
 147    *Gcfrac = f;
 148  }
 149  
 150  //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
 151  // using a recursive integration by parts relation.
 152  // This function cannot be used when mu > 0.
 153  //
 154  // p > 0, integer
 155  // x < 0, |x| < max (1,p-1)
 156  
 157  void double_G_ibp (DOUBLE_T * Gibp, DOUBLE_T p, DOUBLE_T x)
 158  {
 159    BOOL_T odd = (INT_T) (p) % 2 != 0;
 160    DOUBLE_T t = fabsq (x), del;
 161    DOUBLE_T tt = 1.0q / (t * t), c = 1.0q / t, d = (p - 1.0q);
 162    DOUBLE_T s = c * (t - d);
 163    INT_T l = 0;
 164    BOOL_T stop;
 165    do {
 166      c *= d * (d - 1.0q) * tt;
 167      d -= 2.0q;
 168      del = c * (t - d);
 169      s += del;
 170      l++;
 171      stop = fabsq (del) < fabsq (s) * EPS;
 172    }
 173    while ((l < floorq ((p - 2.0q) / 2.0q)) && !stop);
 174    if (odd && !stop) {
 175      s += d * c / t;
 176    }
 177    *Gibp = ((odd ? -1.0q : 1.0q) * expq (-t + lgammaq (p) - (p - 1.0q) * logq (t)) + s) / t;
 178  }
 179  
 180  //! @brief compute the G-function in the domain x > p using a
 181  // continued fraction.
 182  //
 183  // p > 0
 184  // x > p, or x = +infinity
 185  
 186  void double_G_cfrac_upper (DOUBLE_T * Gcfrac, DOUBLE_T p, DOUBLE_T x)
 187  {
 188    DOUBLE_T c, d, del, f, an, bn;
 189    INT_T i, n;
 190  // Special case
 191    if (isinfq (x)) {
 192      *Gcfrac = 0.0q;
 193      return;
 194    }
 195  // Evaluate the continued fraction using Modified Lentz's method. However,
 196  // as detailed in the paper, perform manually the first pass (n=1), of the
 197  // initial Modified Lentz's method.
 198    an = 1.0q;
 199    bn = x + 1.0q - p;
 200    BOOL_T t = (bn != 0.0q);
 201    if (t) {
 202  // b{1} is non-zero
 203      f = an / bn;
 204      c = an / DPMIN;
 205      d = 1.0q / bn;
 206      n = 2;
 207    } else {
 208  // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
 209      an = -(1.0q - p);
 210      bn = x + 3.0q - p;
 211      f = an / bn;
 212      c = an / DPMIN;
 213      d = 1.0q / bn;
 214      n = 3;
 215    }
 216    i = n - 1;
 217    do {
 218      an = -i * (i - p);
 219      bn += 2.0q;
 220      d = an * d + bn;
 221      if (d == 0.0q) {
 222        d = DPMIN;
 223      }
 224      c = bn + an / c;
 225      if (c == 0.0q) {
 226        c = DPMIN;
 227      }
 228      d = 1.0q / d;
 229      del = d * c;
 230      f *= del;
 231      i++;
 232      n++;
 233    }
 234    while ((fabsq (del - 1.0q) >= EPS) && (n < ITMAX));
 235    *Gcfrac = t ? f : 1.0q / f;
 236  }
 237  
 238  //! @brief compute G : (p,x) --> R defined as follows
 239  //
 240  // if x <= p:
 241  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
 242  // otherwise:
 243  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
 244  //
 245  //   p > 0
 246  //   x is a real number or +infinity.
 247  
 248  void double_G_func (DOUBLE_T * G, DOUBLE_T p, DOUBLE_T x)
 249  {
 250    if (p >= double_plim (x)) {
 251      double_G_cfrac_lower (G, p, x);
 252    } else if (x < 0.0q) {
 253      double_G_ibp (G, p, x);
 254    } else {
 255      double_G_cfrac_upper (G, p, x);
 256    }
 257  }
 258  
 259  //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
 260  
 261  void double_romberg_iterations (DOUBLE_T * R, DOUBLE_T sigma, INT_T n, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p, DOUBLE_T h, DOUBLE_T pow2)
 262  {
 263    INT_T adr0_prev = ((n - 1) * n) / 2, adr0 = (n * (n + 1)) / 2;
 264    QUAD_T sum = QUAD_REAL_ZERO;
 265    for (INT_T j = 1; j <= pow2; j++) {
 266      DOUBLE_T xx = x + ((y - x) * (2.0q * j - 1.0q)) / (2.0q * pow2);
 267      QUAD_T f = double_real_to_quad_real (expq (-mu * xx + (p - 1.0q) * logq (xx) - sigma));
 268      sum = _add_quad_real_ (sum, f);
 269    }
 270    R[adr0] = 0.5q * R[adr0_prev] + h * quad_real_to_double_real (sum);
 271    DOUBLE_T pow4 = 4.0q;
 272    for (INT_T m = 1; m <= n; m++) {
 273      R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1.0q);
 274      pow4 *= 4.0q;
 275    }
 276  }
 277  
 278  //! @ compute I_{x,y}^{mu,p} using a Romberg approximation.
 279  // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
 280  
 281  void double_romberg_estimate (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
 282  {
 283    DOUBLE_T *R = (DOUBLE_T *) get_heap_space (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2 * sizeof (DOUBLE_T));
 284    ASSERT (R != NULL);
 285  // Initialization (n=1)
 286    *sigma = -mu * y + (p - 1.0q) * logq (y);
 287    R[0] = 0.5q * (y - x) * (expq (-mu * x + (p - 1.0q) * logq (x) - (*sigma)) + 1.0q);
 288  // Loop for n > 0
 289    DOUBLE_T relneeded = EPS / TOL_ROMBERG;
 290    INT_T adr0 = 0;
 291    INT_T n = 1;
 292    DOUBLE_T h = (y - x) / 2.0q;       // n=1, h = (y-x)/2^n
 293    DOUBLE_T pow2 = 1.0q;              // n=1; pow2 = 2^(n-1)
 294    if (NITERMAX_ROMBERG >= 1) {
 295      DOUBLE_T relerr;
 296      do {
 297        double_romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
 298        h /= 2.0q;
 299        pow2 *= 2.0q;
 300        adr0 = (n * (n + 1)) / 2;
 301        relerr = fabsq ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
 302        n++;
 303      } while (n <= NITERMAX_ROMBERG && relerr > relneeded);
 304    }
 305  // save Romberg estimate and free memory
 306    *rho = R[adr0 + (n - 1)];
 307    a68_free (R);
 308  }
 309  
 310  //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
 311  //
 312  //   I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
 313  //
 314  // This procedure computes (rho, sigma) described below.
 315  // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
 316  //
 317  //   mu is a real number non equal to zero 
 318  //     (in general we take mu = 1 or -1 but any nonzero real number is allowed)
 319  //
 320  //   x, y are two numbers with 0 <= x <= y <= +infinity,
 321  //     (the setting y=+infinity is allowed only when mu > 0)
 322  //
 323  //   p is a real number > 0, p must be an integer when mu < 0.
 324  
 325  void deltagammainc_double_real (DOUBLE_T * rho, DOUBLE_T * sigma, DOUBLE_T x, DOUBLE_T y, DOUBLE_T mu, DOUBLE_T p)
 326  {
 327    DOUBLE_T mA, mB, mx, my, nA, nB, nx, ny;
 328  // Particular cases
 329    if (isinfq (x) && isinfq (y)) {
 330      *rho = 0.0q;
 331      *sigma = a68_dneginf ();
 332      return;
 333    } else if (x == y) {
 334      *rho = 0.0q;
 335      *sigma = a68_dneginf ();
 336      return;
 337    }
 338    if (x == 0.0q && isinfq (y)) {
 339      *rho = 1.0q;
 340      (*sigma) = lgammaq (p) - p * logq (mu);
 341      return;
 342    }
 343  // Initialization
 344    double_G_func (&mx, p, mu * x);
 345    nx = (isinfq (x) ? a68_dneginf () : -mu * x + p * logq (x));
 346    double_G_func (&my, p, mu * y);
 347    ny = (isinfq (y) ? a68_dneginf () : -mu * y + p * logq (y));
 348  
 349  // Compute (mA,nA) and (mB,nB) such as I_{x,y}^{mu,p} can be
 350  // approximated by the difference A-B, where A >= B >= 0, A = mA*exp (nA) an 
 351  // B = mB*exp (nB). When the difference involves more than one digit loss due to
 352  // cancellation errors, the integral I_{x,y}^{mu,p} is evaluated using the
 353  // Romberg approximation method.
 354  
 355    if (mu < 0.0q) {
 356      mA = my;
 357      nA = ny;
 358      mB = mx;
 359      nB = nx;
 360    } else {
 361      if (p < double_plim (mu * x)) {
 362        mA = mx;
 363        nA = nx;
 364        mB = my;
 365        nB = ny;
 366      } else if (p < double_plim (mu * y)) {
 367        mA = 1.0q;
 368        nA = lgammaq (p) - p * logq (mu);
 369        nB = fmax (nx, ny);
 370        mB = mx * expq (nx - nB) + my * expq (ny - nB);
 371      } else {
 372        mA = my;
 373        nA = ny;
 374        mB = mx;
 375        nB = nx;
 376      }
 377    }
 378  // Compute (rho,sigma) such that rho*exp (sigma) = A-B
 379    *rho = mA - mB * expq (nB - nA);
 380    *sigma = nA;
 381  // If the difference involved a significant loss of precision, compute Romberg estimate.
 382    if (!isinfq (y) && ((*rho) / mA < TOL_DIFF)) {
 383      double_romberg_estimate (rho, sigma, x, y, mu, p);
 384    }
 385  }
 386  
 387  // A68G Driver routines
 388  
 389  //! @brief PROC long gamma inc g = (LONG REAL p, x, y, mu) LONG REAL
 390  
 391  void genie_gamma_inc_g_double_real (NODE_T * n)
 392  {
 393    A68_LONG_REAL x, y, mu, p;
 394    POP_OBJECT (n, &mu, A68_LONG_REAL);
 395    POP_OBJECT (n, &y, A68_LONG_REAL);
 396    POP_OBJECT (n, &x, A68_LONG_REAL);
 397    POP_OBJECT (n, &p, A68_LONG_REAL);
 398    DOUBLE_T rho, sigma;
 399    deltagammainc_double_real (&rho, &sigma, VALUE (&x).f, VALUE (&y).f, VALUE (&mu).f, VALUE (&p).f);
 400    PUSH_VALUE (n, dble (rho * expq (sigma)), A68_LONG_REAL);
 401  }
 402  
 403  //! @brief PROC long gamma inc f = (LONG REAL p, x) LONG REAL
 404  
 405  void genie_gamma_inc_f_double_real (NODE_T * n)
 406  {
 407    A68_LONG_REAL x, p;
 408    POP_OBJECT (n, &x, A68_LONG_REAL);
 409    POP_OBJECT (n, &p, A68_LONG_REAL);
 410    DOUBLE_T rho, sigma;
 411    deltagammainc_double_real (&rho, &sigma, VALUE (&x).f, a68_dposinf (), 1.0q, VALUE (&p).f);
 412    PUSH_VALUE (n, dble (rho * expq (sigma)), A68_LONG_REAL);
 413  }
 414  
 415  //! @brief PROC long gamma inc gf = (LONG REAL p, x) LONG REAL
 416  
 417  void genie_gamma_inc_gf_double_real (NODE_T * q)
 418  {
 419  // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
 420  // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
 421    A68_LONG_REAL x, p;
 422    POP_OBJECT (q, &x, A68_LONG_REAL);
 423    POP_OBJECT (q, &p, A68_LONG_REAL);
 424    DOUBLE_T G;
 425    double_G_func (&G, VALUE (&p).f, VALUE (&x).f);
 426    PUSH_VALUE (q, dble (G), A68_LONG_REAL);
 427  }
 428  
 429  //! @brief PROC long gamma inc = (LONG REAL p, x) LONG REAL
 430  
 431  void genie_gamma_inc_h_double_real (NODE_T * n)
 432  {
 433  #if (A68_LEVEL >= 3) && defined (HAVE_GNU_MPFR)
 434    genie_gamma_inc_double_real_mpfr (n);
 435  #else
 436    genie_gamma_inc_f_double_real (n);
 437  #endif
 438  }
 439  
 440  #endif