Algol 68 Genie

mp-gamic.c

     
   1  //! @file mp-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] 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  //
  38  // DELTAGAMMAINC Fast and Accurate Evaluation of a Generalized Incomplete Gamma
  39  // Function. Copyright (C) 2016 Remy Abergel (remy.abergel AT gmail.com), Lionel
  40  // Moisan (Lionel.Moisan AT parisdescartes.fr).
  41  //
  42  // This file is a part of the DELTAGAMMAINC software, dedicated to the
  43  // computation of a generalized incomplete gammafunction. See the Companion paper
  44  // for a complete description of the algorithm.
  45  // ``Fast and accurate evaluation of a generalized incomplete gamma function''
  46  // (Rémy Abergel, Lionel Moisan), preprint MAP5 nº2016-14, revision 1.
  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  #include "a68g-genie.h"
  75  #include "a68g-mp.h"
  76  #include "a68g-prelude.h"
  77  
  78  // Processing time of Abergel's algorithms rises steeply with precision.
  79  #define MAX_PRECISION (LONG_LONG_MP_DIGITS + LONG_MP_DIGITS)
  80  #define GAM_DIGITS(digs) FUN_DIGITS (digs)
  81  
  82  #define ITMAX(p, digits) lit_mp (p, 1000000, 0, digits);
  83  #define DPMIN(p, digits) lit_mp (p, 1, 10 - MAX_MP_EXPONENT, digits)
  84  #define EPS(p, digits) lit_mp (p, 1, 1 - digits, digits)
  85  #define NITERMAX_ROMBERG 16 // Maximum allowed number of Romberg iterations
  86  #define TOL_ROMBERG(p, digits) lit_mp (p, MP_RADIX / 10, -1, digits)
  87  #define TOL_DIFF(p, digits) lit_mp (p, MP_RADIX / 5, -1, digits)
  88  
  89  //! @brief compute G(p,x) in the domain x <= p >= 0 using a continued fraction
  90  
  91  MP_T *G_cfrac_lower_mp (NODE_T *q, MP_T * Gcfrac, MP_T *p, MP_T *x, int digs)
  92  {
  93    if (IS_ZERO_MP (x)) {
  94      SET_MP_ZERO (Gcfrac, digs);
  95      return Gcfrac;
  96    }
  97    ADDR_T pop_sp = A68_SP;
  98    MP_T *c = nil_mp (q, digs);
  99    MP_T *d = nil_mp (q, digs);
 100    MP_T *del = nil_mp (q, digs);
 101    MP_T *f = nil_mp (q, digs);
 102  // Evaluate the continued fraction using Modified Lentz's method. However,
 103  // as detailed in the paper, perform manually the first pass (n=1), of the
 104  // initial Modified Lentz's method.
 105  // an = 1; bn = p; f = an / bn; c = an / DPMIN; d = 1 / bn; n = 2;
 106    MP_T *an = lit_mp (q, 1, 0, digs);
 107    MP_T *bn = nil_mp (q, digs);
 108    MP_T *trm = nil_mp (q, digs);
 109    MP_T *dpmin = DPMIN (q, digs);
 110    MP_T *eps = EPS (q, digs);
 111    MP_T *itmax = ITMAX (q, digs);
 112    (void) move_mp (bn, p, digs);
 113    (void) div_mp (q, f, an, bn, digs);
 114    (void) div_mp (q, c, an, dpmin, digs);
 115    (void) rec_mp (q, d, bn, digs);
 116    MP_T *n = lit_mp (q, 2, 0, digs);
 117    MP_T *k = nil_mp (q, digs);
 118    MP_T *two = lit_mp (q, 2, 0, digs);
 119    BOOL_T odd = A68_FALSE, cont = A68_TRUE;
 120    while (cont) {
 121      A68_BOOL ge, lt;
 122  //  k = n / 2;
 123      (void) over_mp (q, k, n, two, digs);    
 124  //  an = (n & 1 ? k : -(p - 1 + k)) * x;
 125      if (odd) {
 126        (void) move_mp (an, k, digs);
 127        odd = A68_FALSE;
 128      } else {
 129        (void) minus_one_mp (q, trm, p, digs);
 130        (void) add_mp (q, trm, trm, k, digs);
 131        (void) minus_mp (q, an, trm, digs);
 132        odd = A68_TRUE;
 133      }
 134      (void) mul_mp (q, an, an, x, digs);
 135  //  bn++;
 136      (void) plus_one_mp (q, bn, bn, digs);
 137  //  d = an * d + bn;
 138      (void) mul_mp (q, trm, an, d, digs);
 139      (void) add_mp (q, d, trm, bn, digs);
 140  //  if (d == 0) { d = DPMIN; }
 141      if (IS_ZERO_MP (d)) {
 142        (void) move_mp (d, dpmin, digs);
 143      }
 144  //  c = bn + an / c; mind possible overflow.
 145      (void) div_mp (q, trm, an, c, digs);
 146      (void) add_mp (q, c, bn, trm, digs);
 147  //  if (c == 0) { c = DPMIN; }
 148      if (IS_ZERO_MP (c)) {
 149        (void) move_mp (c, dpmin, digs);
 150      }
 151  //  d = 1 / d;
 152      (void) rec_mp (q, d, d, digs);
 153  //  del = d * c;
 154      (void) mul_mp (q, del, d, c, digs);
 155  //  f *= del;
 156      (void) mul_mp (q, f, f, del, digs);
 157  //  n++;
 158      (void) plus_one_mp (q, n, n, digs);
 159  //  while ((fabs_double (del - 1) >= EPS) && (n < ITMAX));
 160      (void) minus_one_mp (q, trm, del, digs);
 161      (void) abs_mp (q, trm, trm, digs);
 162      (void) ge_mp (q, &ge, trm, eps, digs);
 163      (void) lt_mp (q, &lt, n, itmax, digs);
 164      cont = VALUE (&ge) && VALUE (&lt);
 165    }
 166    (void) move_mp (Gcfrac, f, digs);
 167    A68_SP = pop_sp;
 168    return Gcfrac;
 169  }
 170  
 171  //! @brief compute the G-function in the domain x > p using a
 172  // continued fraction.
 173  // 0 < p < x, or x = +infinity
 174  
 175  MP_T *G_cfrac_upper_mp (NODE_T *q, MP_T * Gcfrac, MP_T *p, MP_T *x, int digs)
 176  {
 177    ADDR_T pop_sp = A68_SP;
 178    if (PLUS_INF_MP (x)) {
 179      SET_MP_ZERO (Gcfrac, digs);
 180      return Gcfrac;
 181    }
 182    MP_T *c = nil_mp (q, digs);
 183    MP_T *d = nil_mp (q, digs);
 184    MP_T *del = nil_mp (q, digs);
 185    MP_T *f = nil_mp (q, digs);
 186    MP_T *trm = nil_mp (q, digs);
 187    MP_T *dpmin = DPMIN (q, digs);
 188    MP_T *eps = EPS (q, digs);
 189    MP_T *itmax = ITMAX (q, digs);
 190    MP_T *n = lit_mp (q, 2, 0, digs);
 191    MP_T *i = nil_mp (q, digs);
 192    MP_T *two = lit_mp (q, 2, 0, digs);
 193  // an = 1;
 194    MP_T *an = lit_mp (q, 1, 0, digs);
 195  // bn = x + 1 - p;
 196    MP_T *bn = lit_mp (q, 1, 0, digs);
 197    (void) add_mp (q, bn, x, bn, digs);
 198    (void) sub_mp (q, bn, bn, p, digs);
 199    BOOL_T t = !IS_ZERO_MP (bn);
 200  // Evaluate the continued fraction using Modified Lentz's method. However,
 201  // as detailed in the paper, perform manually the first pass (n=1), of the
 202  // initial Modified Lentz's method.
 203    if (t) {
 204  // b{1} is non-zero
 205  // f = an / bn;
 206      (void) div_mp (q, f, an, bn, digs);
 207  // c = an / DPMIN;
 208      (void) div_mp (q, c, an, dpmin, digs);
 209  // d = 1 / bn;
 210      (void) rec_mp (q, d, bn, digs);
 211  // n = 2;
 212      set_mp (n, 2, 0, digs);
 213    } else {
 214  // b{1}=0 but b{2} is non-zero, compute Mcfrac = a{1}/f with f = a{2}/(b{2}+) a{3}/(b{3}+) ...
 215  //  an = -(1 - p);
 216      (void) minus_one_mp (q, an, p, digs);
 217  //  bn = x + 3 - p;
 218      (void) set_mp (bn, 3, 0, digs);
 219      (void) add_mp (q, bn, x, bn, digs);
 220      (void) sub_mp (q, bn, bn, p, digs);
 221  //  f = an / bn;
 222      (void) div_mp (q, f, an, bn, digs);
 223  //  c = an / DPMIN;
 224      (void) div_mp (q, c, an, dpmin, digs);
 225  //  d = 1 / bn;
 226      (void) rec_mp (q, d, bn, digs);
 227  //  n = 3;
 228      set_mp (n, 3, 0, digs);
 229    }
 230  //  i = n - 1;
 231    minus_one_mp (q, i, n, digs); 
 232    BOOL_T cont = A68_TRUE;  
 233    while (cont) {
 234      A68_BOOL ge, lt;
 235  //  an = -i * (i - p);
 236      (void) sub_mp (q, trm, p, i, digs);
 237      (void) mul_mp (q, an, i, trm, digs);
 238  //  bn += 2;
 239      (void) add_mp (q, bn, bn, two, digs);
 240  //  d = an * d + bn;
 241      (void) mul_mp (q, trm, an, d, digs);
 242      (void) add_mp (q, d, trm, bn, digs);
 243  //  if (d == 0) { d = DPMIN; }
 244      if (IS_ZERO_MP (d)) {
 245        (void) move_mp (d, dpmin, digs);
 246      }
 247  //  c = bn + an / c; mind possible overflow.
 248      (void) div_mp (q, trm, an, c, digs);
 249      (void) add_mp (q, c, bn, trm, digs);
 250  //  if (c == 0) { c = DPMIN; }
 251      if (IS_ZERO_MP (c)) {
 252        (void) move_mp (c, dpmin, digs);
 253      }
 254  //  d = 1 / d;
 255      (void) rec_mp (q, d, d, digs);
 256  //  del = d * c;
 257      (void) mul_mp (q, del, d, c, digs);
 258  //  f *= del;
 259      (void) mul_mp (q, f, f, del, digs);
 260  //  i++;
 261      (void) plus_one_mp (q, i, i, digs);
 262  //  n++;
 263      (void) plus_one_mp (q, n, n, digs);
 264  //  while ((fabs_double (del - 1) >= EPS) && (n < ITMAX));
 265      (void) minus_one_mp (q, trm, del, digs);
 266      (void) abs_mp (q, trm, trm, digs);
 267      (void) ge_mp (q, &ge, trm, eps, digs);
 268      (void) lt_mp (q, &lt, n, itmax, digs);
 269      cont = VALUE (&ge) && VALUE (&lt);
 270    }
 271    A68_SP = pop_sp;
 272  // *Gcfrac = t ? f : 1 / f;
 273    if (t) {
 274      (void) move_mp (Gcfrac, f, digs);
 275    } else {
 276      (void) rec_mp (q, Gcfrac, f, digs);
 277    }
 278    return Gcfrac;
 279  }
 280  
 281  //! @brief compute the G-function in the domain x < 0 and |x| < max (1,p-1)
 282  // using a recursive integration by parts relation.
 283  // This function cannot be used when mu > 0.
 284  // p > 0, integer; x < 0, |x| < max (1,p-1)
 285  
 286  MP_T *G_ibp_mp (NODE_T *q, MP_T * Gibp, MP_T *p, MP_T *x, int digs)
 287  {
 288    ADDR_T pop_sp = A68_SP;
 289    MP_T *trm = nil_mp (q, digs), *trn = nil_mp (q, digs);
 290    MP_T *eps = EPS (q, digs);
 291  // t = fabs_double (x);
 292    MP_T *t = nil_mp (q, digs);
 293    (void) abs_mp (q, x, x, digs);
 294  // tt = 1 / (t * t);
 295    MP_T *tt = nil_mp (q, digs);
 296    (void) mul_mp (q, tt, t, t, digs);
 297    (void) rec_mp (q, tt, tt, digs);
 298  // odd = (INT_T) (p) % 2 != 0;
 299    MP_T *two = lit_mp (q, 2, 0, digs);
 300    (void) trunc_mp (q, trm, p, digs);
 301    (void) mod_mp (q, trm, trm, two, digs);
 302    BOOL_T odd = !IS_ZERO_MP (trm);
 303  // c = 1 / t;
 304    MP_T *c = nil_mp (q, digs);
 305    (void) rec_mp (q, c, t, digs);
 306  // d = (p - 1);
 307    MP_T *d = nil_mp (q, digs);
 308    (void) minus_one_mp (q, d, p, digs);
 309  // s = c * (t - d);
 310    MP_T *s = nil_mp (q, digs);
 311    (void) sub_mp (q, trm, t, d, digs);
 312    (void) mul_mp (q, s, c, trm, digs);
 313  // l = 0;
 314    MP_T *l = nil_mp (q, digs);
 315    BOOL_T cont = A68_TRUE, stop = A68_TRUE;
 316    MP_T *del = nil_mp (q, digs);
 317    while (cont) {
 318  //  c *= d * (d - 1) * tt;
 319      (void) minus_one_mp (q, trm, d, digs);
 320      (void) mul_mp (q, trm, d, trm, digs);
 321      (void) mul_mp (q, trm, trm, tt, digs);
 322      (void) mul_mp (q, c, c, trm, digs);
 323  //  d -= 2;
 324      (void) sub_mp (q, d, d, two, digs);
 325  //  del = c * (t - d);
 326      (void) sub_mp (q, trm, t, d, digs);
 327      (void) mul_mp (q, del, c, trm, digs);
 328  //  s += del;
 329      (void) add_mp (q, s, s, del, digs);
 330  //  l++;
 331      (void) plus_one_mp (q, l, l, digs);
 332  //  stop = fabs_double (del) < fabs_double (s) * EPS;
 333      (void) abs_mp (q, trm, del, digs);
 334      (void) abs_mp (q, trn, s, digs);
 335      (void) mul_mp (q, trn, trn, eps, digs);
 336      A68_BOOL lt;
 337      (void) lt_mp (q, &lt, trm, trn, digs);
 338      stop = VALUE (&lt);
 339  //while ((l < floor_double ((p - 2) / 2)) && !stop);
 340      (void) sub_mp (q, trm, p, two, digs);
 341      (void) half_mp (q, trm, trm, digs);
 342      (void) floor_mp (q, trm, trm, digs);
 343      (void) lt_mp (q, &lt, l, trm, digs);
 344      cont = VALUE (&lt) && !stop;
 345    }
 346    if (odd && !stop) {
 347  //   s += d * c / t;
 348      (void) div_mp (q, trm, c, t, digs);
 349      (void) mul_mp (q, trm, d, trm, digs);
 350      (void) add_mp (q, s, s, trm, digs);
 351    }
 352  // Gibp = ((odd ? -1 : 1) * exp_double (-t + lgammaq (p) - (p - 1) * log_double (t)) + s) / t;
 353    (void) ln_mp (q, trn, t, digs);
 354    (void) minus_one_mp (q, trm, p, digs);
 355    (void) mul_mp (q, trm, trm, trn, digs);
 356    (void) lngamma_mp (q, trn, p, digs);
 357    (void) sub_mp (q, trm, trn, trm, digs);
 358    (void) sub_mp (q, trm, trm, t, digs);
 359    (void) exp_mp (q, Gibp, trm, digs);
 360    if (odd) {
 361      (void) minus_mp (q, Gibp, Gibp, digs);
 362    }
 363    (void) add_mp (q, Gibp, Gibp, s, digs);
 364    (void) div_mp (q, Gibp, Gibp, t, digs);
 365    A68_SP = pop_sp;
 366    return Gibp;
 367  }
 368  
 369  MP_T *plim_mp (NODE_T *p, MP_T *z, MP_T *x, int digs)
 370  {
 371    ADDR_T pop_sp = A68_SP;
 372    if (MP_DIGIT (x, 1) > 0) {
 373      (void) move_mp (z, x, digs);
 374    } else {
 375      MP_T *five = lit_mp (p, 5, 0, digs);
 376      MP_T *nine = lit_mp (p, -9, 0, digs);
 377      A68_BOOL ge;
 378      (void) ge_mp (p, &ge, x, nine, digs);
 379      if (VALUE (&ge)) {
 380        SET_MP_ZERO (z, digs);
 381      } else {
 382        (void) minus_mp (p, z, x, digs);
 383        (void) sqrt_mp (p, z, z, digs);
 384        (void) mul_mp (p, z, five, z, digs);
 385        (void) sub_mp (p, z, z, five, digs);
 386      }
 387    }
 388    A68_SP = pop_sp;
 389    return z;
 390  }
 391  
 392  //! @brief compute G : (p,x) --> R defined as follows
 393  // if x <= p:
 394  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-sign (x)*s) ds from s = 0 to |x|
 395  // otherwise:
 396  //   G(p,x) = exp (x-p*ln (|x|)) * integral of s^{p-1} * exp (-s) ds from s = x to infinity
 397  //   p > 0; x is a real number or +infinity.
 398  
 399  void G_func_mp (NODE_T *q, MP_T * G, MP_T *p, MP_T *x, int digs)
 400  {
 401    ADDR_T pop_sp = A68_SP;
 402    MP_T *pl = nil_mp (q, digs);
 403    A68_BOOL ge;
 404    (void) plim_mp (q, pl, x, digs);
 405    (void) ge_mp (q, &ge, p, pl, digs);
 406    if (VALUE (&ge)) {
 407      G_cfrac_lower_mp (q, G, p, x, digs);
 408    } else if (MP_DIGIT (x, 1) < 0) {
 409      G_ibp_mp (q, G, p, x, digs);
 410    } else {
 411      G_cfrac_upper_mp (q, G, p, x, digs);
 412    }
 413    A68_SP = pop_sp;
 414  }
 415  
 416  //! @brief compute I_{x,y}^{mu,p} using a Romberg approximation.
 417  // Compute rho and sigma so I_{x,y}^{mu,p} = rho * exp (sigma)
 418  
 419  //! @brief iteration of the Romberg approximation of I_{x,y}^{mu,p}
 420  
 421  #define ROMBERG_N (((NITERMAX_ROMBERG + 1) * (NITERMAX_ROMBERG + 2)) / 2)
 422  
 423  static inline int IX (int n, int digs) {
 424    int offset = n * SIZE_MP (digs);
 425    return offset;
 426  }
 427  
 428  void mp_romberg_iterations 
 429    (NODE_T *q, MP_T *R, MP_T *sigma, INT_T n, MP_T *x, MP_T *y, MP_T *mu, MP_T *p, MP_T *h, MP_T *pow2, int digs)
 430  {
 431    MP_T *trm = nil_mp (q, digs), *trn = nil_mp (q, digs);
 432    MP_T *sum = nil_mp (q, digs), *xx = nil_mp (q, digs);
 433    INT_T adr0_prev = ((n - 1) * n) / 2;
 434    INT_T adr0 = (n * (n + 1)) / 2;
 435    MP_T *j = lit_mp (q, 1, 0, digs);
 436    A68_BOOL le;
 437    VALUE (&le) = A68_TRUE;
 438    while (VALUE (&le)) {
 439  //  xx = x + ((y - x) * (2 * j - 1)) / (2 * pow2);
 440      (void) add_mp (q, trm, j, j, digs);
 441      (void) minus_one_mp (q, trm, trm, digs);
 442      (void) sub_mp (q, trn, y, x, digs);
 443      (void) mul_mp (q, trm, trm, trn, digs);
 444      (void) div_mp (q, trm, trm, pow2, digs);
 445      (void) half_mp (q, trm, trm, digs);
 446      (void) add_mp (q, xx, x, trm, digs);
 447  //  sum += exp (-mu * xx + (p - 1) * a68_ln_real (xx) - sigma);
 448      (void) ln_mp (q, trn, xx, digs);
 449      (void) minus_one_mp (q, trm, p, digs);
 450      (void) mul_mp (q, trm, trm, trn, digs);
 451      (void) mul_mp (q, trn, mu, xx, digs);
 452      (void) sub_mp (q, trm, trm, trn, digs);
 453      (void) sub_mp (q, trm, trm, sigma, digs);
 454      (void) exp_mp (q, trm, trm, digs);
 455      (void) add_mp (q, sum, sum, trm, digs);
 456  // j++;
 457      (void) plus_one_mp (q, j, j, digs);
 458      (void) le_mp (q, &le, j, pow2, digs);
 459    }
 460  // R[adr0] = 0.5 * R[adr0_prev] + h * sum;
 461    (void) half_mp (q, trm, &R[IX (adr0_prev, digs)], digs);
 462    (void) mul_mp (q, trn, h, sum, digs);
 463    (void) add_mp (q, &R[IX (adr0, digs)], trm, trn, digs);
 464  // REAL_T pow4 = 4;
 465    MP_T *pow4 = lit_mp (q, 4, 0, digs);
 466    for (int m = 1; m <= n; m++) {
 467  //  R[adr0 + m] = (pow4 * R[adr0 + (m - 1)] - R[adr0_prev + (m - 1)]) / (pow4 - 1);
 468      (void) mul_mp (q, trm, pow4, &R[IX (adr0 + m - 1, digs)], digs);
 469      (void) sub_mp (q, trm, trm, &R[IX (adr0_prev + m - 1, digs)], digs);
 470      (void) minus_one_mp (q, trn, pow4, digs);
 471      (void) div_mp (q, &R[IX (adr0 + m, digs)], trm, trn, digs);
 472  //  pow4 *= 4;
 473      (void) add_mp (q, trm, pow4, pow4, digs);
 474      (void) add_mp (q, pow4, trm, trm, digs);
 475    }
 476  }
 477  
 478  void mp_romberg_estimate (NODE_T *q, MP_T * rho, MP_T * sigma, MP_T *x, MP_T *y, MP_T *mu, MP_T *p, int digs)
 479  {
 480    ADDR_T pop_sp = A68_SP;
 481    MP_T *R = (MP_T *) get_heap_space (ROMBERG_N * SIZE_MP (digs));
 482  // Initialization (n=1)
 483    MP_T *trm = nil_mp (q, digs), *trn = nil_mp (q, digs);
 484  // sigma = -mu * y + (p - 1) * ln (y);
 485    (void) ln_mp (q, trn, y, digs);
 486    (void) minus_one_mp (q, trm, p, digs);
 487    (void) mul_mp (q, trm, trm, trn, digs);
 488    (void) mul_mp (q, trn, mu, y, digs);
 489    (void) sub_mp (q, sigma, trm, trn, digs);
 490  // R[0] = 0.5 * (y - x) * (exp (-mu * x + (p - 1) * ln (x) - (*sigma)) + 1);
 491    (void) ln_mp (q, trn, x, digs);
 492    (void) minus_one_mp (q, trm, p, digs);
 493    (void) mul_mp (q, trm, trm, trn, digs);
 494    (void) mul_mp (q, trn, mu, x, digs);
 495    (void) sub_mp (q, trm, trm, trn, digs);
 496    (void) sub_mp (q, trm, trm, sigma, digs);
 497    (void) exp_mp (q, trm, trm, digs);
 498    (void) plus_one_mp (q, trm, trm, digs);
 499    (void) sub_mp (q, trn, y, x, digs);
 500    (void) mul_mp (q, trm, trm, trn, digs);
 501    (void) half_mp (q, &R[IX (0, digs)], trm, digs);
 502  // Loop for n > 0
 503    MP_T *relerr = nil_mp (q, digs);
 504    MP_T *relneeded = EPS (q, digs);
 505    (void) div_mp (q, relneeded, relneeded, TOL_ROMBERG (q, digs), digs);
 506    INT_T adr0 = 0;
 507    INT_T n = 1;
 508  // REAL_T h = (y - x) / 2;       // n=1, h = (y-x)/2^n
 509    MP_T *h = nil_mp (q, digs);
 510    (void) sub_mp (q, h, y, x, digs);
 511    (void) half_mp (q, h, h, digs);
 512  //  REAL_T pow2 = 1;              // n=1; pow2 = 2^(n-1)
 513    MP_T *pow2 = lit_mp (q, 1, 0, digs);
 514    BOOL_T cont = A68_TRUE;
 515    while (cont) {
 516  // while (n <= NITERMAX_ROMBERG && relerr > relneeded);
 517  //  mp_romberg_iterations (R, *sigma, n, x, y, mu, p, h, pow2);
 518      ADDR_T pop_sp_2 = A68_SP;
 519      mp_romberg_iterations (q, R, sigma, n, x, y, mu, p, h, pow2, digs);
 520      A68_SP = pop_sp_2;
 521  //  h /= 2;
 522      (void) half_mp (q, h, h, digs);
 523  //  pow2 *= 2;
 524      (void) add_mp (q, pow2, pow2, pow2, digs);
 525  //  adr0 = (n * (n + 1)) / 2;
 526      adr0 = (n * (n + 1)) / 2;
 527  //  relerr = abs ((R[adr0 + n] - R[adr0 + n - 1]) / R[adr0 + n]);
 528      (void) sub_mp (q, trm, &R[IX (adr0 + n, digs)], &R[IX (adr0 + n - 1, digs)], digs);
 529      (void) div_mp (q, trm, trm, &R[IX (adr0 + n, digs)], digs);
 530      (void) abs_mp (q, relerr, trm, digs);
 531  //  n++;
 532      n++; 
 533      A68_BOOL gt;
 534      (void) gt_mp (q, &gt, relerr, relneeded, digs);
 535      cont = (n <= NITERMAX_ROMBERG) && VALUE (&gt);
 536    }
 537  // save Romberg estimate and free memory
 538  // rho = R[adr0 + (n - 1)];
 539    (void) move_mp (rho, &R[IX (adr0 + n - 1, digs)], digs);
 540    a68_free (R);
 541    A68_SP = pop_sp;
 542  }
 543  
 544  //! @brief compute generalized incomplete gamma function I_{x,y}^{mu,p}
 545  //   I_{x,y}^{mu,p} = integral from x to y of s^{p-1} * exp (-mu*s) ds
 546  // This procedure computes (rho, sigma) described below.
 547  // The approximated value of I_{x,y}^{mu,p} is I = rho * exp (sigma)
 548  //   mu is a real number non equal to zero 
 549  //     (in general we take mu = 1 or -1 but any nonzero real number is allowed)
 550  //   x, y are two numbers with 0 <= x <= y <= +infinity,
 551  //     (the setting y=+infinity is allowed only when mu > 0)
 552  //   p is a real number > 0, p must be an integer when mu < 0.
 553  
 554  void Dgamic_mp (NODE_T *q, MP_T * rho, MP_T * sigma, MP_T *x, MP_T *y, MP_T *mu, MP_T *p, int digs)
 555  {
 556    ADDR_T pop_sp = A68_SP;
 557  // Particular cases
 558    if (PLUS_INF_MP (x) && PLUS_INF_MP (y)) {
 559      SET_MP_ZERO (rho, digs);
 560      SET_MP_ZERO (sigma, digs);
 561      MP_STATUS (sigma) = (UNSIGNED_T) MP_STATUS (sigma) | MINUS_INF_MASK;
 562      A68_SP = pop_sp;
 563      return;
 564    } else if (same_mp (q, x, y, digs)) {
 565      SET_MP_ZERO (rho, digs);
 566      SET_MP_ZERO (sigma, digs);
 567      MP_STATUS (sigma) = (UNSIGNED_T) MP_STATUS (sigma) | MINUS_INF_MASK;
 568      A68_SP = pop_sp;
 569      return;
 570    }
 571    if (IS_ZERO_MP (x) && PLUS_INF_MP (y)) {
 572      set_mp (rho, 1, 0, digs);
 573      MP_T *lgam = nil_mp (q, digs);
 574      MP_T *lnmu = nil_mp (q, digs);
 575      (void) lngamma_mp (q, lgam, p, digs);
 576      (void) ln_mp (q, lnmu, mu, digs);
 577      (void) mul_mp (q, lnmu, p, lnmu, digs);
 578      (void) sub_mp (q, sigma, lgam, lnmu, digs);
 579      return;
 580    }
 581  // Initialization
 582    MP_T *mx = nil_mp (q, digs);
 583    MP_T *nx = nil_mp (q, digs);
 584    MP_T *my = nil_mp (q, digs);
 585    MP_T *ny = nil_mp (q, digs);
 586    MP_T *mux = nil_mp (q, digs);
 587    MP_T *muy = nil_mp (q, digs);
 588  // Initialization
 589  // nx = (a68_isinf_real (x) ? a68_neginf_double_real () : -mu * x + p * ln (x));
 590    if (PLUS_INF_MP (x)) {
 591      SET_MP_ZERO (mx, digs);
 592      MP_STATUS (nx) = (UNSIGNED_T) MP_STATUS (nx) | MINUS_INF_MASK;
 593    } else {
 594      (void) mul_mp (q, mux, mu, x, digs);
 595      G_func_mp (q, mx, p, mux, digs);
 596      (void) ln_mp (q, nx, x, digs);
 597      (void) mul_mp (q, nx, p, nx, digs);
 598      (void) sub_mp (q, nx, nx, mux, digs);
 599    }
 600  // ny = (a68_isinf_real (y) ? a68_neginf_double_real () : -mu * y + p * ln (y));
 601    if (PLUS_INF_MP (y)) {
 602      SET_MP_ZERO (my, digs);
 603      MP_STATUS (ny) = (UNSIGNED_T) MP_STATUS (ny) | MINUS_INF_MASK;
 604    } else {
 605      (void) mul_mp (q, muy, mu, y, digs);
 606      G_func_mp (q, my, p, muy, digs);
 607      (void) ln_mp (q, ny, y, digs);
 608      (void) mul_mp (q, ny, p, ny, digs);
 609      (void) sub_mp (q, ny, ny, muy, digs);
 610    }
 611  // Compute (mA,nA) and (mB,nB) such as I_{x,y}^{mu,p} can be
 612  // approximated by the difference A-B, where A >= B >= 0, A = mA*exp (nA) an 
 613  // B = mB*exp (nB). When the difference involves more than one digit loss due to
 614  // cancellation errors, the integral I_{x,y}^{mu,p} is evaluated using the
 615  // Romberg approximation method.
 616    MP_T *mA = nil_mp (q, digs);
 617    MP_T *mB = nil_mp (q, digs);
 618    MP_T *nA = nil_mp (q, digs);
 619    MP_T *nB = nil_mp (q, digs);
 620    MP_T *trm = nil_mp (q, digs);
 621    if (MP_DIGIT (mu, 1) < 0) {
 622      (void) move_mp (mA, my, digs);
 623      (void) move_mp (nA, ny, digs);
 624      (void) move_mp (mB, mx, digs);
 625      (void) move_mp (nB, nx, digs);
 626      goto compute;
 627    }
 628    MP_T *pl = nil_mp (q, digs);
 629    A68_BOOL lt;
 630    if (PLUS_INF_MP (x)) {
 631      VALUE (&lt) = A68_TRUE;
 632    } else {
 633      (void) mul_mp (q, mux, mu, x, digs);
 634      (void) plim_mp (q, pl, mux, digs);
 635      (void) lt_mp (q, &lt, p, pl, digs);
 636    }
 637    if (VALUE (&lt)) {
 638      (void) move_mp (mA, mx, digs);
 639      (void) move_mp (nA, nx, digs);
 640      (void) move_mp (mB, my, digs);
 641      (void) move_mp (nB, ny, digs);
 642      goto compute;
 643    }
 644    if (PLUS_INF_MP (y)) {
 645      VALUE (&lt) = A68_TRUE;
 646    } else {
 647      (void) mul_mp (q, muy, mu, y, digs);
 648      (void) plim_mp (q, pl, muy, digs);
 649      (void) lt_mp (q, &lt, p, pl, digs);
 650    }
 651    if (VALUE (&lt)) {
 652  //  mA = 1;
 653      set_mp (mA, 1, 0, digs);
 654  //  nA = lgammaq (p) - p * log_double (mu);
 655      MP_T *lgam = nil_mp (q, digs);
 656      MP_T *lnmu = nil_mp (q, digs);
 657      (void) lngamma_mp (q, lgam, p, digs);
 658      (void) ln_mp (q, lnmu, mu, digs);
 659      (void) mul_mp (q, lnmu, p, lnmu, digs);
 660      (void) sub_mp (q, nA, lgam, lnmu, digs);
 661  //  nB = fmax (nx, ny);
 662      A68_BOOL ge;
 663      if (MINUS_INF_MP (ny)) {
 664        VALUE (&ge) = A68_TRUE;
 665      } else {
 666        (void) ge_mp (q, &ge, nx, ny, digs);
 667      }
 668      if (VALUE (&ge)) {
 669        (void) move_mp (nB, nx, digs);
 670      } else {
 671        (void) move_mp (nB, ny, digs);
 672      }
 673  //  mB = mx * exp_double (nx - nB) + my * exp_double (ny - nB);
 674      (void) sub_mp (q, trm, nx, nB, digs);
 675      (void) exp_mp (q, trm, trm, digs);
 676      (void) mul_mp (q, mB, mx, trm, digs);
 677      if (! MINUS_INF_MP (ny)) {
 678        (void) sub_mp (q, trm, ny, nB, digs);
 679        (void) exp_mp (q, trm, trm, digs);
 680        (void) mul_mp (q, trm, my, trm, digs);
 681        (void) add_mp (q, mB, nB, trm, digs);
 682      }
 683      goto compute;
 684    }
 685    (void) move_mp (mA, my, digs);
 686    (void) move_mp (nA, ny, digs);
 687    (void) move_mp (mB, mx, digs);
 688    (void) move_mp (nB, nx, digs);
 689    compute:
 690  // Compute (rho,sigma) such that rho*exp (sigma) = A-B
 691  // 1. rho = mA - mB * exp_double (nB - nA);
 692    (void) sub_mp (q, trm, nB, nA, digs);
 693    (void) exp_mp (q, trm, trm, digs);
 694    (void) mul_mp (q, trm, mB, trm, digs);
 695    (void) sub_mp (q, rho, mA, trm, digs);
 696  // 2. sigma = nA;
 697    (void) move_mp (sigma, nA, digs);
 698  // If the difference involved a significant loss of precision, compute Romberg estimate.
 699  // if (!isinf_double (y) && ((*rho) / mA < TOL_DIFF)) {
 700    (void) div_mp (q, trm, rho, mA, digs);
 701    (void) lt_mp (q, &lt, trm, TOL_DIFF (q, digs), digs);
 702    if (!PLUS_INF_MP (y) && VALUE (&lt)) {
 703      mp_romberg_estimate (q, rho, sigma, x, y, mu, p, digs);
 704    }
 705    A68_SP = pop_sp;
 706  }
 707  
 708  void Dgamic_wrap_mp (NODE_T *q, MP_T * s, MP_T * rho, MP_T * sigma, MP_T *x, MP_T *y, MP_T *mu, MP_T *p, int digs)
 709  {
 710    ADDR_T pop_sp = A68_SP;
 711    int gdigs = GAM_DIGITS (MAX_PRECISION);
 712    errno = 0;
 713    if (digs <= gdigs) {
 714      gdigs = GAM_DIGITS (digs);
 715      MP_T *rho_g = len_mp (q, rho, digs, gdigs);
 716      MP_T *sigma_g = len_mp (q, sigma, digs, gdigs);
 717      MP_T *x_g = len_mp (q, x, digs, gdigs);
 718      MP_T *y_g = len_mp (q, y, digs, gdigs);
 719      MP_T *mu_g = len_mp (q, mu, digs, gdigs);
 720      MP_T *p_g = len_mp (q, p, digs, gdigs);
 721      Dgamic_mp (q, rho_g, sigma_g, x_g, y_g, mu_g, p_g, gdigs);
 722      if (IS_ZERO_MP (rho_g) || MINUS_INF_MP (sigma_g)) {
 723        SET_MP_ZERO (s, digs);
 724      } else {
 725        (void) exp_mp (q, sigma_g, sigma_g, gdigs);
 726        (void) mul_mp (q, rho_g, rho_g, sigma_g, gdigs);
 727        (void) shorten_mp (q, s, digs, rho_g, gdigs);
 728      }
 729    } else {
 730      diagnostic (A68_MATH_WARNING, q, WARNING_MATH_PRECISION, MOID (q), CALL, NULL);
 731      MP_T *rho_g = cut_mp (q, rho, digs, gdigs);
 732      MP_T *sigma_g = cut_mp (q, sigma, digs, gdigs);
 733      MP_T *x_g = cut_mp (q, x, digs, gdigs);
 734      MP_T *y_g = cut_mp (q, y, digs, gdigs);
 735      MP_T *mu_g = cut_mp (q, mu, digs, gdigs);
 736      MP_T *p_g = cut_mp (q, p, digs, gdigs);
 737      Dgamic_mp (q, rho_g, sigma_g, x_g, y_g, mu_g, p_g, gdigs);
 738      if (IS_ZERO_MP (rho_g) || MINUS_INF_MP (sigma_g)) {
 739        SET_MP_ZERO (s, digs);
 740      } else {
 741        (void) exp_mp (q, sigma_g, sigma_g, gdigs);
 742        (void) mul_mp (q, rho_g, rho_g, sigma_g, gdigs);
 743        MP_T *tmp = nil_mp (q, MAX_PRECISION);
 744        (void) shorten_mp (q, tmp, MAX_PRECISION, rho_g, gdigs);
 745        (void) lengthen_mp (q, s, digs, tmp, MAX_PRECISION);
 746      }
 747    }
 748    PRELUDE_ERROR (errno != 0, q, ERROR_MATH, MOID (q));
 749    A68_SP = pop_sp;
 750  }
 751  
 752  //! @brief PROC long long gamma inc f = (LONG LONG REAL p, x) LONG LONG REAL
 753  
 754  void genie_gamma_inc_f_mp (NODE_T *p)
 755  {
 756    int digs = DIGITS (MOID (p)), size = SIZE (MOID (p));
 757    ADDR_T pop_sp = A68_SP;
 758    MP_T *x = (MP_T *) STACK_OFFSET (-size);
 759    MP_T *s = (MP_T *) STACK_OFFSET (-2 * size);
 760    MP_T *mu = lit_mp (p, 1, 0, digs);
 761    MP_T *y = nil_mp (p, digs);
 762    MP_T *rho = nil_mp (p, digs);
 763    MP_T *sigma = nil_mp (p, digs);
 764    MP_STATUS (y) = (UNSIGNED_T) MP_STATUS (y) | PLUS_INF_MASK;
 765    Dgamic_wrap_mp (p, s, rho, sigma, x, y, mu, s, digs);
 766    A68_SP = pop_sp;
 767    A68_SP -= size;
 768  }
 769  
 770  //! @brief PROC long long gamma inc g = (LONG LONG REAL p, x, y, mu) LONG LONG REAL
 771  
 772  void genie_gamma_inc_g_mp (NODE_T *p)
 773  {
 774    int digs = DIGITS (MOID (p)), size = SIZE (MOID (p));
 775    ADDR_T pop_sp = A68_SP;
 776    MP_T *mu = (MP_T *) STACK_OFFSET (-size);
 777    MP_T *y = (MP_T *) STACK_OFFSET (-2 * size);
 778    MP_T *x = (MP_T *) STACK_OFFSET (-3 * size);
 779    MP_T *s = (MP_T *) STACK_OFFSET (-4 * size);
 780    MP_T *rho = nil_mp (p, digs);
 781    MP_T *sigma = nil_mp (p, digs);
 782    Dgamic_wrap_mp (p, s, rho, sigma, x, y, mu, s, digs);
 783    A68_SP = pop_sp;
 784    A68_SP -= 3 * size;
 785  }
 786  
 787  //! @brief PROC long long gamma inc gf = (LONG LONG REAL p, x, y, mu) LONG LONG REAL
 788  
 789  void genie_gamma_inc_gf_mp (NODE_T *p)
 790  {
 791  // if x <= p: G(p,x) = exp (x-p*ln (|x|)) * integral over [0,|x|] of s^{p-1} * exp (-sign (x)*s) ds
 792  // otherwise: G(p,x) = exp (x-p*ln (x)) * integral over [x,inf] of s^{p-1} * exp (-s) ds
 793    int digs = DIGITS (MOID (p)), size = SIZE (MOID (p));
 794    ADDR_T pop_sp = A68_SP;
 795    MP_T *x = (MP_T *) STACK_OFFSET (-size);
 796    MP_T *s = (MP_T *) STACK_OFFSET (-2 * size);
 797    int gdigs = GAM_DIGITS (MAX_PRECISION);
 798    errno = 0;
 799    if (digs <= gdigs) {
 800      gdigs = GAM_DIGITS (digs);
 801      MP_T *x_g = len_mp (p, x, digs, gdigs);
 802      MP_T *s_g = len_mp (p, s, digs, gdigs);
 803      MP_T *G = nil_mp (p, gdigs);
 804      G_func_mp (p, G, s_g, x_g, gdigs);
 805      PRELUDE_ERROR (errno != 0, p, ERROR_MATH, MOID (p));
 806      (void) shorten_mp (p, s, digs, G, gdigs);
 807    } else {
 808      diagnostic (A68_MATH_WARNING, p, WARNING_MATH_PRECISION, MOID (p), CALL, NULL);
 809      MP_T *x_g = cut_mp (p, x, digs, gdigs);
 810      MP_T *s_g = cut_mp (p, s, digs, gdigs);
 811      MP_T *G = nil_mp (p, gdigs);
 812      G_func_mp (p, G, s_g, x_g, gdigs);
 813      PRELUDE_ERROR (errno != 0, p, ERROR_MATH, MOID (p));
 814      MP_T *tmp = nil_mp (p, MAX_PRECISION);
 815      (void) shorten_mp (p, tmp, MAX_PRECISION, G, gdigs);
 816      (void) lengthen_mp (p, s, digs, tmp, MAX_PRECISION);
 817    }
 818    A68_SP = pop_sp - size;
 819  }
 820  
 821  //! @brief PROC long long gamma inc = (LONG LONG REAL p, x) LONG LONG REAL
 822  
 823  void genie_gamma_inc_h_mp (NODE_T *p)
 824  {
 825  // #if defined (HAVE_GNU_MPFR) && (A68_LEVEL >= 3)
 826  //   genie_gamma_inc_mpfr (p);
 827  // #else
 828  //   genie_gamma_inc_f_mp (p);
 829  // #endif
 830     genie_gamma_inc_f_mp (p);
 831  }