single-math.c

     
   1  //! @file single-math.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 math stuff supplementing libc.
  25  
  26  // References:
  27  //
  28  //   Milton Abramowitz and Irene Stegun, Handbook of Mathematical Functions,
  29  //   Dover Publications, New York [1970]
  30  //   https://en.wikipedia.org/wiki/Abramowitz_and_Stegun
  31  
  32  #include "a68g.h"
  33  #include "a68g-genie.h"
  34  #include "a68g-prelude.h"
  35  #include "a68g-double.h"
  36  #include "a68g-numbers.h"
  37  #include "a68g-math.h"
  38  
  39  inline REAL_T a68_max (REAL_T x, REAL_T y)
  40  {
  41    return (x > y ? x : y);
  42  }
  43  
  44  inline REAL_T a68_min (REAL_T x, REAL_T y)
  45  {
  46    return (x < y ? x : y);
  47  }
  48  
  49  inline REAL_T a68_sign (REAL_T x)
  50  {
  51    return (x == 0 ? 0 : (x > 0 ? 1 : -1));
  52  }
  53  
  54  inline REAL_T a68_int (REAL_T x)
  55  {
  56    return (x >= 0 ? (INT_T) x : -(INT_T) - x);
  57  }
  58  
  59  inline INT_T a68_round (REAL_T x)
  60  {
  61    return (INT_T) (x >= 0 ? x + 0.5 : x - 0.5);
  62  }
  63  
  64  #define IS_INTEGER(n) (n == a68_int (n))
  65  
  66  inline REAL_T a68_abs (REAL_T x)
  67  {
  68    return (x >= 0 ? x : -x);
  69  }
  70  
  71  REAL_T a68_fdiv (REAL_T x, REAL_T y)
  72  {
  73  // This is for generating +-INF.
  74    return x / y;
  75  }
  76  
  77  REAL_T a68_nan (void)
  78  {
  79    return a68_fdiv (0.0, 0.0);
  80  }
  81  
  82  REAL_T a68_posinf (void)
  83  {
  84    return a68_fdiv (+1.0, 0.0);
  85  }
  86  
  87  REAL_T a68_neginf (void)
  88  {
  89    return a68_fdiv (-1.0, 0.0);
  90  }
  91  
  92  // REAL infinity
  93  
  94  void genie_infinity_real (NODE_T * p)
  95  {
  96    PUSH_VALUE (p, a68_posinf (), A68_REAL);
  97  }
  98  
  99  // REAL minus infinity
 100  
 101  void genie_minus_infinity_real (NODE_T * p)
 102  {
 103    PUSH_VALUE (p, a68_neginf (), A68_REAL);
 104  }
 105  
 106  int a68_finite (REAL_T x)
 107  {
 108  #if defined (HAVE_ISFINITE)
 109    return isfinite (x);
 110  #elif defined (HAVE_FINITE)
 111    return finite (x);
 112  #else
 113    (void) x;
 114    return A68_TRUE;
 115  #endif
 116  }
 117  
 118  int a68_isnan (REAL_T x)
 119  {
 120  #if defined (HAVE_ISNAN)
 121    return isnan (x);
 122  #elif defined (HAVE_IEEE_COMPARISONS)
 123    int status = (x != x);
 124    return status;
 125  #else
 126    return A68_FALSE;
 127  #endif
 128  }
 129  
 130  int a68_isinf (REAL_T x)
 131  {
 132  #if defined (HAVE_ISINF)
 133    if (isinf (x)) {
 134      return (x > 0) ? 1 : -1;
 135    } else {
 136      return 0;
 137    }
 138  #else
 139    if (!a68_finite (x) && !a68_isnan (x)) {
 140      return (x > 0 ? 1 : -1);
 141    } else {
 142      return 0;
 143    }
 144  #endif
 145  }
 146  
 147  // INT operators
 148  
 149  INT_T a68_add_int (INT_T j, INT_T k)
 150  {
 151    if (j >= 0) {
 152      A68_OVERFLOW (A68_MAX_INT - j < k);
 153    } else {
 154      A68_OVERFLOW (k < (-A68_MAX_INT) - j);
 155    }
 156    return j + k;
 157  }
 158  
 159  INT_T a68_sub_int (INT_T j, INT_T k)
 160  {
 161    return a68_add_int (j, -k);
 162  }
 163  
 164  INT_T a68_mul_int (INT_T j, INT_T k)
 165  {
 166    if (j == 0 || k == 0) {
 167      return 0;
 168    } else {
 169      INT_T u = (j > 0 ? j : -j), v = (k > 0 ? k : -k);
 170      A68_OVERFLOW (u > A68_MAX_INT / v);
 171      return j * k;
 172    }
 173  }
 174  
 175  INT_T a68_over_int (INT_T j, INT_T k)
 176  {
 177    A68_INVALID (k == 0);
 178    return j / k;
 179  }
 180  
 181  INT_T a68_mod_int (INT_T j, INT_T k)
 182  {
 183    A68_INVALID (k == 0);
 184    INT_T r = j % k;
 185    return (r < 0 ? (k > 0 ? r + k : r - k) : r);
 186  }
 187  
 188  // OP ** = (INT, INT) INT 
 189  
 190  INT_T a68_m_up_n (INT_T m, INT_T n)
 191  {
 192  // Only positive n.
 193    A68_INVALID (n < 0);
 194  // Special cases.
 195    if (m == 0 || m == 1) {
 196      return m;
 197    } else if (m == -1) {
 198      return (EVEN (n) ? 1 : -1);
 199    }
 200  // General case with overflow check.
 201    UNSIGNED_T bit = 1;
 202    INT_T M = m, P = 1;
 203    do {
 204      if (n & bit) {
 205        P = a68_mul_int (P, M);
 206      }
 207      bit <<= 1;
 208      if (bit <= n) {
 209        M = a68_mul_int (M, M);
 210      }
 211    } while (bit <= n);
 212    return P;
 213  }
 214  
 215  // OP ** = (REAL, INT) REAL 
 216  
 217  REAL_T a68_x_up_n (REAL_T x, INT_T n)
 218  {
 219  // Only positive n.
 220    if (n < 0) {
 221      return 1 / a68_x_up_n (x, -n);
 222    }
 223  // Special cases.
 224    if (x == 0 || x == 1) {
 225      return x;
 226    } else if (x == -1) {
 227      return (EVEN (n) ? 1 : -1);
 228    }
 229  // General case.
 230    UNSIGNED_T bit = 1;
 231    REAL_T M = x, P = 1;
 232    do {
 233      if (n & bit) {
 234        P *= M;
 235      }
 236      bit <<= 1;
 237      if (bit <= n) {
 238        M *= M;
 239      }
 240    } while (bit <= n);
 241    A68_OVERFLOW (!finite (P));
 242    return P;
 243  }
 244  
 245  REAL_T a68_div_int (INT_T j, INT_T k)
 246  {
 247    A68_INVALID (k == 0);
 248    return (REAL_T) j / (REAL_T) k;
 249  }
 250  
 251  // Sqrt (x^2 + y^2) that does not needlessly overflow.
 252  
 253  REAL_T a68_hypot (REAL_T x, REAL_T y)
 254  {
 255    REAL_T xabs = ABS (x), yabs = ABS (y), min, max;
 256    if (xabs < yabs) {
 257      min = xabs;
 258      max = yabs;
 259    } else {
 260      min = yabs;
 261      max = xabs;
 262    }
 263    if (min == 0) {
 264      return max;
 265    } else {
 266      REAL_T u = min / max;
 267      return max * sqrt (1 + u * u);
 268    }
 269  }
 270  
 271  //! @brief Compute Chebyshev series to requested accuracy.
 272  
 273  REAL_T a68_chebyshev (REAL_T x, const REAL_T c[], REAL_T acc)
 274  {
 275  // Iteratively compute the recursive Chebyshev series.
 276  // c[1..N] are coefficients, c[0] is N, and acc is relative accuracy.
 277    acc *= MATH_EPSILON;
 278    if (acc < c[1]) {
 279      diagnostic (A68_MATH_WARNING, A68 (f_entry), WARNING_MATH_ACCURACY, NULL);
 280    }
 281    INT_T i, N = a68_round (c[0]);
 282    REAL_T err = 0, z = 2 * x, u = 0, v = 0, w = 0;
 283    for (i = 1; i <= N; i++) {
 284      if (err > acc) {
 285        w = v;
 286        v = u;
 287        u = z * v - w + c[i];
 288      }
 289      err += a68_abs (c[i]);
 290    }
 291    return 0.5 * (u - w);
 292  }
 293  
 294  // Compute ln (1 + x) accurately. 
 295  // Some C99 platforms implement this incorrectly.
 296  
 297  REAL_T a68_ln1p (REAL_T x)
 298  {
 299  // Based on GNU GSL's gsl_sf_log_1plusx_e.
 300    A68_INVALID (x <= -1);
 301    if (a68_abs (x) < pow (DBL_EPSILON, 1 / 6.0)) {
 302      const REAL_T c1 = -0.5, c2 = 1 / 3.0, c3 = -1 / 4.0, c4 = 1 / 5.0, c5 = -1 / 6.0, c6 = 1 / 7.0, c7 = -1 / 8.0, c8 = 1 / 9.0, c9 = -1 / 10.0;
 303      const REAL_T t = c5 + x * (c6 + x * (c7 + x * (c8 + x * c9)));
 304      return x * (1 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * t)))));
 305    } else if (a68_abs (x) < 0.5) {
 306      REAL_T t = (8 * x + 1) / (x + 2) / 2;
 307      return x * a68_chebyshev (t, c_ln1p, 0.1);
 308    } else {
 309      return ln (1 + x);
 310    }
 311  }
 312  
 313  // Compute ln (x), if possible accurately when x ~ 1.
 314  
 315  REAL_T a68_ln (REAL_T x)
 316  {
 317    A68_INVALID (x <= 0);
 318  #if (A68_LEVEL >= 3)
 319    if (a68_abs (x - 1) < 0.375) {
 320  // Double precision x-1 mitigates cancellation error.
 321      return a68_ln1p (DBLEQ (x) - 1.0q);
 322    } else {
 323      return ln (x);
 324    }
 325  #else
 326    return ln (x);
 327  #endif
 328  }
 329  
 330  // PROC (REAL) REAL exp
 331  
 332  REAL_T a68_exp (REAL_T x)
 333  {
 334    A68_INVALID (x < LOG_DBL_MIN || x > LOG_DBL_MAX);
 335    return exp (x);
 336  }
 337  
 338  // OP ** = (REAL, REAL) REAL
 339  
 340  REAL_T a68_x_up_y (REAL_T x, REAL_T y)
 341  {
 342    return a68_exp (y * a68_ln (x));
 343  }
 344  
 345  // PROC (REAL) REAL csc
 346  
 347  REAL_T a68_csc (REAL_T x)
 348  {
 349    REAL_T z = sin (x);
 350    A68_OVERFLOW (z == 0);
 351    return 1 / z;
 352  }
 353  
 354  // PROC (REAL) REAL acsc
 355  
 356  REAL_T a68_acsc (REAL_T x)
 357  {
 358    A68_OVERFLOW (x == 0);
 359    return asin (1 / x);
 360  }
 361  
 362  // PROC (REAL) REAL sec
 363  
 364  REAL_T a68_sec (REAL_T x)
 365  {
 366    REAL_T z = cos (x);
 367    A68_OVERFLOW (z == 0);
 368    return 1 / z;
 369  }
 370  
 371  // PROC (REAL) REAL asec
 372  
 373  REAL_T a68_asec (REAL_T x)
 374  {
 375    A68_OVERFLOW (x == 0);
 376    return acos (1 / x);
 377  }
 378  
 379  // PROC (REAL) REAL cot
 380  
 381  REAL_T a68_cot (REAL_T x)
 382  {
 383    REAL_T z = sin (x);
 384    A68_OVERFLOW (z == 0);
 385    return cos (x) / z;
 386  }
 387  
 388  // PROC (REAL) REAL acot
 389  
 390  REAL_T a68_acot (REAL_T x)
 391  {
 392    A68_OVERFLOW (x == 0);
 393    return atan (1 / x);
 394  }
 395  
 396  // PROC atan2 (REAL, REAL) REAL
 397  
 398  REAL_T a68_atan2 (REAL_T x, REAL_T y)
 399  {
 400    if (x == 0) {
 401      A68_INVALID (y == 0);
 402      return (y > 0 ? CONST_PI_2 : -CONST_PI_2);
 403    } else {
 404      REAL_T z = atan (ABS (y / x));
 405      if (x < 0) {
 406        z = CONST_PI - z;
 407      }
 408      return (y >= 0 ? z : -z);
 409    }
 410  }
 411  
 412  //! brief PROC (REAL) REAL sindg
 413  
 414  REAL_T a68_sindg (REAL_T x)
 415  {
 416    return sin (x * CONST_PI_OVER_180);
 417  }
 418  
 419  //! brief PROC (REAL) REAL cosdg
 420  
 421  REAL_T a68_cosdg (REAL_T x)
 422  {
 423    return cos (x * CONST_PI_OVER_180);
 424  }
 425  
 426  //! brief PROC (REAL) REAL tandg
 427  
 428  REAL_T a68_tandg (REAL_T x)
 429  {
 430    return tan (x * CONST_PI_OVER_180);
 431  }
 432  
 433  //! brief PROC (REAL) REAL asindg
 434  
 435  REAL_T a68_asindg (REAL_T x)
 436  {
 437    return asin (x) * CONST_180_OVER_PI;
 438  }
 439  
 440  //! brief PROC (REAL) REAL acosdg
 441  
 442  REAL_T a68_acosdg (REAL_T x)
 443  {
 444    return acos (x) * CONST_180_OVER_PI;
 445  }
 446  
 447  //! brief PROC (REAL) REAL atandg
 448  
 449  REAL_T a68_atandg (REAL_T x)
 450  {
 451    return atan (x) * CONST_180_OVER_PI;
 452  }
 453  
 454  // PROC (REAL) REAL cotdg
 455  
 456  REAL_T a68_cotdg (REAL_T x)
 457  {
 458    REAL_T z = a68_sindg (x);
 459    A68_OVERFLOW (z == 0);
 460    return a68_cosdg (x) / z;
 461  }
 462  
 463  // PROC (REAL) REAL acotdg
 464  
 465  REAL_T a68_acotdg (REAL_T z)
 466  {
 467    A68_OVERFLOW (z == 0);
 468    return a68_atandg (1 / z);
 469  }
 470  
 471  // @brief PROC (REAL) REAL sinpi
 472  
 473  REAL_T a68_sinpi (REAL_T x)
 474  {
 475    x = fmod (x, 2);
 476    if (x <= -1) {
 477      x += 2;
 478    } else if (x > 1) {
 479      x -= 2;
 480    }
 481  // x in <-1, 1].
 482    if (x == 0 || x == 1) {
 483      return 0;
 484    } else if (x == 0.5) {
 485      return 1;
 486    }
 487    if (x == -0.5) {
 488      return -1;
 489    } else {
 490      return sin (CONST_PI * x);
 491    }
 492  }
 493  
 494  // @brief PROC (REAL) REAL cospi
 495  
 496  REAL_T a68_cospi (REAL_T x)
 497  {
 498    x = fmod (fabs (x), 2);
 499  // x in [0, 2>.
 500    if (x == 0.5 || x == 1.5) {
 501      return 0;
 502    } else if (x == 0) {
 503      return 1;
 504    } else if (x == 1) {
 505      return -1;
 506    } else {
 507      return cos (CONST_PI * x);
 508    }
 509  }
 510  
 511  // @brief PROC (REAL) REAL tanpi
 512  
 513  REAL_T a68_tanpi (REAL_T x)
 514  {
 515    x = fmod (x, 1);
 516    if (x <= -0.5) {
 517      x += 1;
 518    } else if (x > 0.5) {
 519      x -= 1;
 520    }
 521  // x in <-1/2, 1/2].
 522    A68_OVERFLOW (x == 0.5);
 523    if (x == -0.25) {
 524      return -1;
 525    } else if (x == 0) {
 526      return 0;
 527    } else if (x == 0.25) {
 528      return 1;
 529    } else {
 530      return a68_sinpi (x) / a68_cospi (x);
 531    }
 532  }
 533  
 534  // @brief PROC (REAL) REAL cotpi
 535  
 536  REAL_T a68_cotpi (REAL_T x)
 537  {
 538    x = fmod (x, 1);
 539    if (x <= -0.5) {
 540      x += 1;
 541    } else if (x > 0.5) {
 542      x -= 1;
 543    }
 544  // x in <-1/2, 1/2].
 545    A68_OVERFLOW (x == 0);
 546    if (x == -0.25) {
 547      return -1;
 548    } else if (x == 0.25) {
 549      return 1;
 550    } else if (x == 0.5) {
 551      return 0;
 552    } else {
 553      return a68_cospi (x) / a68_sinpi (x);
 554    }
 555  }
 556  
 557  // @brief PROC (REAL) REAL asinh
 558  
 559  REAL_T a68_asinh (REAL_T x)
 560  {
 561    REAL_T a = ABS (x), s = (x < 0 ? -1.0 : 1);
 562    if (a > 1 / sqrt (DBL_EPSILON)) {
 563      return (s * (a68_ln (a) + a68_ln (2)));
 564    } else if (a > 2) {
 565      return (s * a68_ln (2 * a + 1 / (a + sqrt (a * a + 1))));
 566    } else if (a > sqrt (DBL_EPSILON)) {
 567      REAL_T a2 = a * a;
 568      return (s * a68_ln1p (a + a2 / (1 + sqrt (1 + a2))));
 569    } else {
 570      return (x);
 571    }
 572  }
 573  
 574  // @brief PROC (REAL) REAL acosh
 575  
 576  REAL_T a68_acosh (REAL_T x)
 577  {
 578    if (x > 1 / sqrt (DBL_EPSILON)) {
 579      return (a68_ln (x) + a68_ln (2));
 580    } else if (x > 2) {
 581      return (a68_ln (2 * x - 1 / (sqrt (x * x - 1) + x)));
 582    } else if (x > 1) {
 583      REAL_T t = x - 1;
 584      return (a68_ln1p (t + sqrt (2 * t + t * t)));
 585    } else if (x == 1) {
 586      return (0);
 587    } else {
 588      A68_INVALID (A68_TRUE);
 589    }
 590  }
 591  
 592  // @brief PROC (REAL) REAL atanh
 593  
 594  REAL_T a68_atanh (REAL_T x)
 595  {
 596    REAL_T a = ABS (x);
 597    A68_INVALID (a >= 1);
 598    REAL_T s = (x < 0 ? -1 : 1);
 599    if (a >= 0.5) {
 600      return (s * 0.5 * a68_ln1p (2 * a / (1 - a)));
 601    } else if (a > DBL_EPSILON) {
 602      return (s * 0.5 * a68_ln1p (2 * a + 2 * a * a / (1 - a)));
 603    } else {
 604      return (x);
 605    }
 606  }
 607  
 608  //! @brief Inverse complementary error function.
 609  
 610  REAL_T a68_inverfc (REAL_T y)
 611  {
 612    A68_INVALID (y < 0 || y > 2);
 613    if (y == 0) {
 614      return DBL_MAX;
 615    } else if (y == 1) {
 616      return 0;
 617    } else if (y == 2) {
 618      return -DBL_MAX;
 619    } else {
 620  // Next is based on code that originally contained following statement:
 621  //   Copyright (c) 1996 Takuya Ooura. You may use, copy, modify this 
 622  //   code for any purpose and without fee.
 623      REAL_T s, t, u, v, x, z;
 624      z = (y <= 1 ? y : 2 - y);
 625      v = c_inverfc[0] - a68_ln (z);
 626      u = sqrt (v);
 627      s = (a68_ln (u) + c_inverfc[1]) / v;
 628      t = 1 / (u + c_inverfc[2]);
 629      x = u * (1 - s * (s * c_inverfc[3] + 0.5)) - ((((c_inverfc[4] * t + c_inverfc[5]) * t + c_inverfc[6]) * t + c_inverfc[7]) * t + c_inverfc[8]) * t;
 630      t = c_inverfc[9] / (x + c_inverfc[9]);
 631      u = t - 0.5;
 632      s = (((((((((c_inverfc[10] * u + c_inverfc[11]) * u - c_inverfc[12]) * u - c_inverfc[13]) * u + c_inverfc[14]) * u + c_inverfc[15]) * u - c_inverfc[16]) * u - c_inverfc[17]) * u + c_inverfc[18]) * u + c_inverfc[19]) * u + c_inverfc[20];
 633      s = ((((((((((((s * u - c_inverfc[21]) * u - c_inverfc[22]) * u + c_inverfc[23]) * u + c_inverfc[24]) * u + c_inverfc[25]) * u + c_inverfc[26]) * u + c_inverfc[27]) * u + c_inverfc[28]) * u + c_inverfc[29]) * u + c_inverfc[30]) * u + c_inverfc[31]) * u + c_inverfc[32]) * t - z * a68_exp (x * x - c_inverfc[33]);
 634      x += s * (x * s + 1);
 635      return (y <= 1 ? x : -x);
 636    }
 637  }
 638  
 639  //! @brief Inverse error function.
 640  
 641  REAL_T a68_inverf (REAL_T y)
 642  {
 643    return a68_inverfc (1 - y);
 644  }
 645  
 646  //! @brief PROC (REAL, REAL) REAL ln beta
 647  
 648  REAL_T a68_ln_beta (REAL_T a, REAL_T b)
 649  {
 650    return lgamma (a) + lgamma (b) - lgamma (a + b);
 651  }
 652  
 653  //! @brief PROC (REAL, REAL) REAL beta
 654  
 655  REAL_T a68_beta (REAL_T a, REAL_T b)
 656  {
 657    return a68_exp (a68_ln_beta (a, b));
 658  }
 659  
 660  //! brief PROC (INT) REAL fact
 661  
 662  REAL_T a68_fact (INT_T n)
 663  {
 664    A68_INVALID (n < 0 || n > A68_MAX_FAC);
 665    return factable[n];
 666  }
 667  
 668  //! brief PROC (INT) REAL ln fact
 669  
 670  REAL_T a68_ln_fact (INT_T n)
 671  {
 672    A68_INVALID (n < 0);
 673    if (n <= A68_MAX_FAC) {
 674      return ln_factable[n];
 675    } else {
 676      return lgamma (n + 1);
 677    }
 678  }
 679  
 680  //! @brief PROC choose = (INT n, m) REAL
 681  
 682  REAL_T a68_choose (INT_T n, INT_T m)
 683  {
 684    A68_INVALID (n < m);
 685    return factable[n] / (factable[m] * factable[n - m]);
 686  }
 687  
 688  //! @brief PROC ln choose = (INT n, m) REAL
 689  
 690  REAL_T a68_ln_choose (INT_T n, INT_T m)
 691  {
 692    A68_INVALID (n < m);
 693    return a68_ln_fact (n) - (a68_ln_fact (m) + a68_ln_fact (n - m));
 694  }
 695  
 696  REAL_T a68_beta_inc (REAL_T s, REAL_T t, REAL_T x)
 697  {
 698  // Incomplete beta function I{x}(s, t).
 699  // Continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
 700    if (x < 0 || x > 1) {
 701      errno = ERANGE;
 702      return -1;
 703    } else {
 704      const INT_T lim = 16 * sizeof (REAL_T);
 705      BOOL_T cont = A68_TRUE;
 706  // Rapid convergence when x <= (s+1)/(s+t+2) or else recursion.
 707      if (x > (s + 1) / (s + t + 2)) {
 708  // B{x}(s, t) = 1 - B{1-x}(t, s)
 709        return 1 - a68_beta_inc (s, t, 1 - x);
 710      }
 711  // Lentz's algorithm for continued fraction.
 712      REAL_T W = 1, F = 1, c = 1, d = 0;
 713      INT_T N, m;
 714      for (N = 0, m = 0; cont && N < lim; N++) {
 715        REAL_T T;
 716        if (N == 0) {
 717          T = 1;
 718        } else if (N % 2 == 0) {
 719  // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
 720          T = x * m * (t - m) / (s + 2 * m - 1) / (s + 2 * m);
 721        } else {
 722  // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
 723          T = -x * (s + m) * (s + t + m) / (s + 2 * m + 1) / (s + 2 * m);
 724          m++;
 725        }
 726        d = 1 / (T * d + 1);
 727        c = T / c + 1;
 728        F *= c * d;
 729        if (F == W) {
 730          cont = A68_FALSE;
 731        } else {
 732          W = F;
 733        }
 734      }
 735  // I{x}(s,t)=x^s(1-x)^t / s / B(s,t) F
 736      REAL_T beta = a68_exp (lgamma (s) + lgamma (t) - lgamma (s + t));
 737      return a68_x_up_y (x, s) * a68_x_up_y (1 - x, t) / s / beta * (F - 1);
 738    }
 739  }