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


© 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)