mp-gamma.c

     
   1  //! @file mp-gamma.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 error, gamma and beta functions.
  25  
  26  #include "a68g.h"
  27  #include "a68g-genie.h"
  28  #include "a68g-prelude.h"
  29  #include "a68g-double.h"
  30  #include "a68g-mp.h"
  31  #include "a68g-lib.h"
  32  
  33  //! @brief PROC (LONG REAL) LONG REAL erf
  34  
  35  MP_T *erf_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
  36  {
  37    if (IS_ZERO_MP (x)) {
  38      SET_MP_ZERO (z, digs);
  39      return z;
  40    } else {
  41      ADDR_T pop_sp = A68_SP;
  42  // Note we need double precision!
  43      int gdigs = FUN_DIGITS (2 * digs), sign;
  44      BOOL_T go_on = A68_TRUE;
  45      MP_T *y_g = nil_mp (p, gdigs);
  46      MP_T *z_g = len_mp (p, x, digs, gdigs);
  47      sign = MP_SIGN (x);
  48      (void) abs_mp (p, z_g, z_g, gdigs);
  49      (void) set_mp (y_g, gdigs * LOG_MP_RADIX, 0, gdigs);
  50      (void) sqrt_mp (p, y_g, y_g, gdigs);
  51      (void) sub_mp (p, y_g, z_g, y_g, gdigs);
  52      if (MP_SIGN (y_g) >= 0) {
  53        SET_MP_ONE (z, digs);
  54      } else {
  55  // Taylor expansion.
  56        MP_T *p_g = nil_mp (p, gdigs);
  57        MP_T *r_g = nil_mp (p, gdigs);
  58        MP_T *s_g = nil_mp (p, gdigs);
  59        MP_T *t_g = nil_mp (p, gdigs);
  60        MP_T *u_g = nil_mp (p, gdigs);
  61        (void) mul_mp (p, y_g, z_g, z_g, gdigs);
  62        SET_MP_ONE (s_g, gdigs);
  63        SET_MP_ONE (t_g, gdigs);
  64        for (unt k = 1; go_on; k++) {
  65          (void) mul_mp (p, t_g, y_g, t_g, gdigs);
  66          (void) div_mp_digit (p, t_g, t_g, (MP_T) k, gdigs);
  67          (void) div_mp_digit (p, u_g, t_g, (MP_T) (2 * k + 1), gdigs);
  68          if (EVEN (k)) {
  69            (void) add_mp (p, s_g, s_g, u_g, gdigs);
  70          } else {
  71            (void) sub_mp (p, s_g, s_g, u_g, gdigs);
  72          }
  73          go_on = (MP_EXPONENT (s_g) - MP_EXPONENT (u_g)) < gdigs;
  74        }
  75        (void) mul_mp (p, r_g, z_g, s_g, gdigs);
  76        (void) mul_mp_digit (p, r_g, r_g, 2, gdigs);
  77        (void) mp_pi (p, p_g, MP_SQRT_PI, gdigs);
  78        (void) div_mp (p, r_g, r_g, p_g, gdigs);
  79        (void) shorten_mp (p, z, digs, r_g, gdigs);
  80      }
  81      if (sign < 0) {
  82        (void) minus_mp (p, z, z, digs);
  83      }
  84      A68_SP = pop_sp;
  85      return z;
  86    }
  87  }
  88  
  89  //! @brief PROC (LONG REAL) LONG REAL erfc
  90  
  91  MP_T *erfc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
  92  {
  93    (void) erf_mp (p, z, x, digs);
  94    (void) one_minus_mp (p, z, z, digs);
  95    return z;
  96  }
  97  
  98  //! @brief PROC (LONG REAL) LONG REAL inverf
  99  
 100  MP_T *inverf_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 101  {
 102    ADDR_T pop_sp = A68_SP;
 103    int gdigs, sign;
 104    BOOL_T go_on = A68_TRUE;
 105  // Precision adapts to the argument, but not too much.
 106  // If this is not precise enough, you need more digs
 107  // in your entire calculation, not just in this routine.
 108  // Calculate an initial Newton-Raphson estimate while at it.
 109  #if (A68_LEVEL >= 3)
 110    DOUBLE_T y = ABS (mp_to_double_real (p, x, digs));
 111    if (y < erfq (5.0q)) {
 112      y = inverf_double_real (y);
 113      gdigs = FUN_DIGITS (digs);
 114    } else {
 115      y = 5.0q;
 116      gdigs = FUN_DIGITS (2 * digs);
 117    }
 118    MP_T *z_g = nil_mp (p, gdigs);
 119    (void) double_real_to_mp (p, z_g, y, gdigs);
 120  #else
 121    REAL_T y = ABS (mp_to_real (p, x, digs));
 122    if (y < erf (4.0)) {
 123      y = a68_inverf (y);
 124      gdigs = FUN_DIGITS (digs);
 125    } else {
 126      y = 4.0;
 127      gdigs = FUN_DIGITS (2 * digs);
 128    }
 129    MP_T *z_g = nil_mp (p, gdigs);
 130    (void) real_to_mp (p, z_g, y, gdigs);
 131  #endif
 132    MP_T *x_g = len_mp (p, x, digs, gdigs);
 133    MP_T *y_g = nil_mp (p, gdigs);
 134    sign = MP_SIGN (x);
 135    (void) abs_mp (p, x_g, x_g, gdigs);
 136    SET_MP_ONE (y_g, gdigs);
 137    (void) sub_mp (p, y_g, x_g, y_g, gdigs);
 138    if (MP_SIGN (y_g) >= 0) {
 139      errno = EDOM;
 140      A68_SP = pop_sp;
 141      return NaN_MP;
 142    } else {
 143  // Newton-Raphson.
 144      MP_T *d_g = nil_mp (p, gdigs);
 145      MP_T *f_g = nil_mp (p, gdigs);
 146      MP_T *p_g = nil_mp (p, gdigs);
 147  // sqrt (pi) / 2
 148      (void) mp_pi (p, p_g, MP_SQRT_PI, gdigs);
 149      (void) half_mp (p, p_g, p_g, gdigs);
 150  // nrdigs prevents end-less iteration
 151      int nrdigs;
 152      for (nrdigs = 0; nrdigs < digs && go_on; nrdigs++) {
 153        (void) move_mp (y_g, z_g, gdigs);
 154        (void) erf_mp (p, f_g, z_g, gdigs);
 155        (void) sub_mp (p, f_g, f_g, x_g, gdigs);
 156        (void) mul_mp (p, d_g, z_g, z_g, gdigs);
 157        (void) minus_mp (p, d_g, d_g, gdigs);
 158        (void) exp_mp (p, d_g, d_g, gdigs);
 159        (void) div_mp (p, f_g, f_g, d_g, gdigs);
 160        (void) mul_mp (p, f_g, f_g, p_g, gdigs);
 161        (void) sub_mp (p, z_g, z_g, f_g, gdigs);
 162        (void) sub_mp (p, y_g, z_g, y_g, gdigs);
 163        if (IS_ZERO_MP (y_g)) {
 164          go_on = A68_FALSE;
 165        } else {
 166          go_on = ABS (MP_EXPONENT (y_g)) < digs;
 167        }
 168      }
 169      (void) shorten_mp (p, z, digs, z_g, gdigs);
 170    }
 171    if (sign < 0) {
 172      (void) minus_mp (p, z, z, digs);
 173    }
 174    A68_SP = pop_sp;
 175    return z;
 176  }
 177  
 178  //! @brief PROC (LONG REAL) LONG REAL inverfc
 179  
 180  MP_T *inverfc_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 181  {
 182    (void) one_minus_mp (p, z, x, digs);
 183    (void) inverf_mp (p, z, z, digs);
 184    return z;
 185  }
 186  
 187  // Reference: 
 188  //   John L. Spouge. "Computation of the Gamma, Digamma, and Trigamma Functions". 
 189  //   SIAM Journal on Numerical Analysis. 31 (3) [1994]
 190  //
 191  // Spouge's algorithm sums terms of greatly varying magnitude.
 192  
 193  #define GAMMA_PRECISION(z) (2 * (z))
 194  
 195  //! brief Set up gamma coefficient table
 196  
 197  void mp_gamma_table (NODE_T *p, int digs)
 198  {
 199    if (A68_MP (mp_gamma_size) <= 0) {
 200      int b = 1;
 201      int gdigs = GAMMA_PRECISION (digs);
 202      REAL_T log_lim = -digs * LOG_MP_RADIX, log_error; 
 203      do {
 204        ABEND (b >= MP_RADIX, ERROR_HIGH_PRECISION, __func__);
 205  // error = 1 / (sqrt (b) * a68_x_up_y (2 * M_PI, b + 0.5));
 206        log_error = -(log10 (b) / 2 + (b + 0.5) * log10 (2 *M_PI));
 207        b += 1;
 208      } while (log_error > log_lim);
 209      A68_MP (mp_gamma_size) = b;
 210      A68_MP (mp_gam_ck) = (MP_T **) get_heap_space ((size_t) ((b + 1) * sizeof (MP_T *)));
 211      A68_MP (mp_gam_ck)[0] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
 212      (void) mp_pi (p, (A68_MP (mp_gam_ck)[0]), MP_SQRT_TWO_PI, gdigs);
 213      ADDR_T pop_sp = A68_SP;
 214      MP_T *ak = nil_mp (p, gdigs);
 215      MP_T *db = lit_mp (p, b, 0, gdigs);
 216      MP_T *ck = nil_mp (p, gdigs);
 217      MP_T *dk = nil_mp (p, gdigs);
 218      MP_T *dz = nil_mp (p, gdigs);
 219      MP_T *hlf = nil_mp (p, gdigs);
 220      MP_T *fac = lit_mp (p, 1, 0, gdigs);
 221      SET_MP_HALF (hlf, gdigs);
 222      for (unt k = 1; k < b; k++) {
 223        set_mp (dk, k, 0, gdigs);
 224        (void) sub_mp (p, ak, db, dk, gdigs);
 225        (void) sub_mp (p, dz, dk, hlf, gdigs);
 226        (void) pow_mp (p, ck, ak, dz, gdigs);
 227        (void) exp_mp (p, dz, ak, gdigs);
 228        (void) mul_mp (p, ck, ck, dz, gdigs);
 229        (void) div_mp (p, ck, ck, fac, gdigs);
 230        A68_MP (mp_gam_ck)[k] = (MP_T *) get_heap_space (SIZE_MP (gdigs));
 231        (void) move_mp ((A68_MP(mp_gam_ck)[k]), ck, gdigs);
 232        (void) mul_mp (p, fac, fac, dk, gdigs);
 233        (void) minus_mp (p, fac, fac, gdigs);
 234      }
 235      A68_SP = pop_sp;
 236    }
 237  }
 238  
 239  MP_T *mp_spouge_sum (NODE_T *p, MP_T *sum, MP_T *x_g, int gdigs)
 240  {
 241    ADDR_T pop_sp = A68_SP;
 242    int a = A68_MP (mp_gamma_size);
 243    MP_T *da = nil_mp (p, gdigs);
 244    MP_T *dz = nil_mp (p, gdigs);
 245    (void) move_mp (sum, A68_MP (mp_gam_ck)[0], gdigs);
 246  // Sum small to large to preserve precision.
 247    for (int k = a - 1; k > 0; k--) {
 248      set_mp (da, k, 0, gdigs);
 249      (void) add_mp (p, dz, x_g, da, gdigs);
 250      (void) div_mp (p, dz, A68_MP (mp_gam_ck)[k], dz, gdigs);
 251      (void) add_mp (p, sum, sum, dz, gdigs);
 252    }
 253    A68_SP = pop_sp;
 254    return sum;
 255  }
 256  
 257  //! @brief PROC (LONG REAL) LONG REAL gamma
 258  
 259  MP_T *gamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 260  {
 261    mp_gamma_table (p, digs);
 262    int gdigs = GAMMA_PRECISION (digs);
 263  // Set up coefficient table.
 264    ADDR_T pop_sp = A68_SP;
 265    if (MP_DIGIT (x, 1) < 0) {
 266  // G(1-x)G(x) = pi / sin (pi*x)
 267      MP_T *pi = nil_mp (p, digs);
 268      MP_T *sz = nil_mp (p, digs);
 269      MP_T *xm = nil_mp (p, digs);
 270      (void) mp_pi (p, pi, MP_PI, digs);
 271      (void) one_minus_mp (p, xm, x, digs);
 272      (void) gamma_mp (p, xm, xm, digs);
 273      (void) sinpi_mp (p, sz, x, digs);
 274      (void) mul_mp (p, sz, sz, xm, digs);
 275      (void) div_mp (p, z, pi, sz, digs);
 276      A68_SP = pop_sp;
 277      return z;
 278    }
 279    int a = A68_MP (mp_gamma_size);
 280  // Compute Spouge's Gamma
 281    MP_T *sum = nil_mp (p, gdigs);
 282    MP_T *x_g = len_mp (p, x, digs, gdigs);
 283    (void) minus_one_mp (p, x_g, x_g, gdigs);
 284    (void) mp_spouge_sum (p, sum, x_g, gdigs);
 285  // (z+a)^(z+0.5)*exp(-(z+a)) * Sum
 286    MP_T *fac = nil_mp (p, gdigs);
 287    MP_T *dz = nil_mp (p, gdigs);
 288    MP_T *az = nil_mp (p, gdigs);
 289    MP_T *da = nil_mp (p, gdigs);
 290    MP_T *hlf = nil_mp (p, gdigs);
 291    SET_MP_HALF (hlf, gdigs);
 292    set_mp (da, a, 0, gdigs);
 293    (void) add_mp (p, az, x_g, da, gdigs);
 294    (void) add_mp (p, dz, x_g, hlf, gdigs);
 295    (void) pow_mp (p, fac, az, dz, gdigs);
 296    (void) minus_mp (p, az, az, gdigs);
 297    (void) exp_mp (p, dz, az, gdigs);
 298    (void) mul_mp (p, fac, fac, dz, gdigs);
 299    (void) mul_mp (p, fac, sum, fac, gdigs);
 300    (void) shorten_mp (p, z, digs, fac, gdigs);
 301    A68_SP = pop_sp;
 302    return z;
 303  }
 304  
 305  //! @brief PROC (LONG REAL) LONG REAL ln gamma
 306  
 307  MP_T *lngamma_mp (NODE_T * p, MP_T * z, MP_T * x, int digs)
 308  {
 309    mp_gamma_table (p, digs);
 310    int gdigs = GAMMA_PRECISION (digs);
 311  // Set up coefficient table.
 312    ADDR_T pop_sp = A68_SP;
 313    if (MP_DIGIT (x, 1) < 0) {
 314  // G(1-x)G(x) = pi / sin (pi*x)
 315      MP_T *sz = nil_mp (p, digs);
 316      MP_T *dz = nil_mp (p, digs);
 317      MP_T *xm = nil_mp (p, digs);
 318      (void) mp_pi (p, dz, MP_LN_PI, digs);
 319      (void) sinpi_mp (p, sz, x, digs);
 320      (void) ln_mp (p, sz, sz, digs);
 321      (void) sub_mp (p, dz, dz, sz, digs);
 322      (void) one_minus_mp (p, xm, x, digs);
 323      (void) lngamma_mp (p, xm, xm, digs);
 324      (void) sub_mp (p, z, dz, xm, digs);
 325      A68_SP = pop_sp;
 326      return z;
 327    }
 328    int a = A68_MP (mp_gamma_size);
 329  // Compute Spouge's ln Gamma
 330    MP_T *sum = nil_mp (p, gdigs);
 331    MP_T *x_g = len_mp (p, x, digs, gdigs);
 332    (void) minus_one_mp (p, x_g, x_g, gdigs);
 333    (void) mp_spouge_sum (p, sum, x_g, gdigs);
 334  // (x+0.5) * ln (x+a) - (x+a) + ln Sum
 335    MP_T *da = nil_mp (p, gdigs);
 336    MP_T *hlf = nil_mp (p, gdigs);
 337    SET_MP_HALF (hlf, gdigs);
 338    MP_T *fac = nil_mp (p, gdigs);
 339    MP_T *dz = nil_mp (p, gdigs);
 340    MP_T *az = nil_mp (p, gdigs);
 341    set_mp (da, a, 0, gdigs);
 342    (void) add_mp (p, az, x_g, da, gdigs);
 343    (void) ln_mp (p, fac, az, gdigs);
 344    (void) add_mp (p, dz, x_g, hlf, gdigs);
 345    (void) mul_mp (p, fac, fac, dz, gdigs);
 346    (void) sub_mp (p, fac, fac, az, gdigs);
 347    (void) ln_mp (p, dz, sum, gdigs);
 348    (void) add_mp (p, fac, fac, dz, gdigs);
 349    (void) shorten_mp (p, z, digs, fac, gdigs);
 350    A68_SP = pop_sp;
 351    return z;
 352  }
 353  
 354  //! @brief PROC (LONG REAL, LONG REAL) LONG REAL ln beta
 355  
 356  MP_T *lnbeta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
 357  {
 358    ADDR_T pop_sp = A68_SP;
 359    MP_T *aa = nil_mp (p, digs);
 360    MP_T *bb = nil_mp (p, digs);
 361    MP_T *ab = nil_mp (p, digs);
 362    (void) lngamma_mp (p, aa, a, digs);
 363    (void) lngamma_mp (p, bb, b, digs);
 364    (void) add_mp (p, ab, a, b, digs);
 365    (void) lngamma_mp (p, ab, ab, digs);
 366    (void) add_mp (p, z, aa, bb, digs);
 367    (void) sub_mp (p, z, z, ab, digs);
 368    A68_SP = pop_sp;
 369    return z;
 370  }
 371  
 372  //! @brief PROC (LONG REAL, LONG REAL) LONG REAL beta
 373  
 374  MP_T *beta_mp (NODE_T * p, MP_T * z, MP_T * a, MP_T *b, int digs)
 375  {
 376    ADDR_T pop_sp = A68_SP;
 377    MP_T *u = nil_mp (p, digs);
 378    lnbeta_mp (p, u, a, b, digs);
 379    exp_mp (p, z, u, digs);
 380    A68_SP = pop_sp;
 381    return z;
 382  }
 383  
 384  //! @brief PROC (LONG REAL, LONG REAL, LONG REAL) LONG REAL beta inc
 385  
 386  MP_T *beta_inc_mp (NODE_T * p, MP_T * z, MP_T * s, MP_T *t, MP_T *x, int digs)
 387  {
 388  // Incomplete beta function I{x}(s, t).
 389  // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
 390    ADDR_T pop_sp = A68_SP;
 391    A68_BOOL gt;
 392    MP_T *one = lit_mp (p, 1, 0, digs);
 393    gt_mp (p, &gt, x, one, digs);
 394    if (MP_DIGIT (x, 1) < 0 || VALUE (&gt)) {
 395      errno = EDOM;
 396      A68_SP = pop_sp;
 397      return NaN_MP;
 398    }
 399    if (same_mp (p, x, one, digs)) {
 400      SET_MP_ONE (z, digs);
 401      A68_SP = pop_sp;
 402      return z;
 403    }
 404  // Rapid convergence when x <= (s+1)/((s+1)+(t+1)) or else recursion.
 405    {
 406      MP_T *u = nil_mp (p, digs), *v = nil_mp (p, digs), *w = nil_mp (p, digs);
 407      (void) plus_one_mp (p, u, s, digs);
 408      (void) plus_one_mp (p, v, t, digs);
 409      (void) add_mp (p, w, u, v, digs);
 410      (void) div_mp (p, u, u, w, digs);
 411      gt_mp (p, &gt, x, u, digs);
 412      if (VALUE (&gt)) {
 413  // B{x}(s, t) = 1 - B{1-x}(t, s)
 414        (void) one_minus_mp (p, x, x, digs);
 415        PRELUDE_ERROR (beta_inc_mp (p, z, s, t, x, digs) == NaN_MP, p, ERROR_INVALID_ARGUMENT, MOID (p));
 416        (void) one_minus_mp (p, z, z, digs);
 417        A68_SP = pop_sp;
 418        return z;
 419      }
 420    }
 421  // Lentz's algorithm for continued fraction.
 422    A68_SP = pop_sp;
 423    int gdigs = FUN_DIGITS (digs);
 424    const INT_T lim = gdigs * LOG_MP_RADIX;
 425    BOOL_T cont = A68_TRUE;
 426    MP_T *F = lit_mp (p, 1, 0, gdigs);
 427    MP_T *T = lit_mp (p, 1, 0, gdigs);
 428    MP_T *W = lit_mp (p, 1, 0, gdigs);
 429    MP_T *c = lit_mp (p, 1, 0, gdigs);
 430    MP_T *d = nil_mp (p, gdigs);
 431    MP_T *m = nil_mp (p, gdigs);
 432    MP_T *s_g = len_mp (p, s, digs, gdigs);
 433    MP_T *t_g = len_mp (p, t, digs, gdigs);
 434    MP_T *x_g = len_mp (p, x, digs, gdigs);
 435    MP_T *u = lit_mp (p, 1, 0, gdigs);
 436    MP_T *v = nil_mp (p, gdigs);
 437    MP_T *w = nil_mp (p, gdigs);
 438    for (unt N = 0; cont && N < lim; N++) {
 439      if (N == 0) {
 440        SET_MP_ONE (T, gdigs);
 441      } else if (N % 2 == 0) {
 442  // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
 443  // T = x * m * (t - m) / (s + 2.0q * m - 1.0q) / (s + 2.0q * m);
 444        (void) sub_mp (p, u, t_g, m, gdigs);
 445        (void) mul_mp (p, u, u, m, gdigs);
 446        (void) mul_mp (p, u, u, x_g, gdigs);
 447        (void) add_mp (p, v, m, m, gdigs);
 448        (void) add_mp (p, v, s_g, v, gdigs);
 449        (void) minus_one_mp (p, v, v, gdigs);
 450        (void) add_mp (p, w, m, m, gdigs);
 451        (void) add_mp (p, w, s_g, w, gdigs);
 452        (void) div_mp (p, T, u, v, gdigs);
 453        (void) div_mp (p, T, T, w, gdigs);
 454      } else {
 455  // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
 456  // T = -x * (s + m) * (s + t + m) / (s + 2.0q * m + 1.0q) / (s + 2.0q * m);
 457        (void) add_mp (p, u, s_g, m, gdigs);
 458        (void) add_mp (p, v, u, t_g, gdigs);
 459        (void) mul_mp (p, u, u, v, gdigs);
 460        (void) mul_mp (p, u, u, x_g, gdigs);
 461        (void) minus_mp (p, u, u, gdigs);
 462        (void) add_mp (p, v, m, m, gdigs);
 463        (void) add_mp (p, v, s_g, v, gdigs);
 464        (void) plus_one_mp (p, v, v, gdigs);
 465        (void) add_mp (p, w, m, m, gdigs);
 466        (void) add_mp (p, w, s_g, w, gdigs);
 467        (void) div_mp (p, T, u, v, gdigs);
 468        (void) div_mp (p, T, T, w, gdigs);
 469        (void) plus_one_mp (p, m, m, gdigs);
 470      }
 471  // d = 1.0q / (T * d + 1.0q);
 472      (void) mul_mp (p, d, T, d, gdigs);
 473      (void) plus_one_mp (p, d, d, gdigs);
 474      (void) rec_mp (p, d, d, gdigs);
 475  // c = T / c + 1.0q;
 476      (void) div_mp (p, c, T, c, gdigs);
 477      (void) plus_one_mp (p, c, c, gdigs);
 478  // F *= c * d;
 479      (void) mul_mp (p, F, F, c, gdigs);
 480      (void) mul_mp (p, F, F, d, gdigs);
 481      if (same_mp (p, F, W, gdigs)) {
 482        cont = A68_FALSE;
 483      } else {
 484        (void) move_mp (W, F, gdigs);
 485      }
 486    }
 487    minus_one_mp (p, F, F, gdigs);
 488  // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
 489    (void) pow_mp (p, u, x_g, s_g, gdigs);
 490    (void) one_minus_mp (p, v, x_g, gdigs);
 491    (void) pow_mp (p, v, v, t_g, gdigs);
 492    (void) beta_mp (p, w, s_g, t_g, gdigs);
 493    (void) mul_mp (p, m, u, v, gdigs);
 494    (void) div_mp (p, m, m, s_g, gdigs);
 495    (void) div_mp (p, m, m, w, gdigs);
 496    (void) mul_mp (p, m, m, F, gdigs);
 497    (void) shorten_mp (p, z, digs, m, gdigs);
 498    A68_SP = pop_sp;
 499    return z;
 500  }