mp-mpfr.c

     
   1  //! @file mp-mpfr.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 routines using GNU MPFR.
  25  
  26  #include "a68g.h"
  27  #include "a68g-genie.h"
  28  #include "a68g-prelude.h"
  29  #include "a68g-mp.h"
  30  #include "a68g-double.h"
  31  
  32  #if (A68_LEVEL >= 3)
  33  
  34  #if defined (HAVE_GNU_MPFR)
  35  
  36  #define DEFAULT GMP_RNDN
  37  
  38  #define MPFR_REAL_BITS (REAL_MANT_DIG)
  39  #define MPFR_LONG_REAL_BITS (FLT128_MANT_DIG)
  40  #define MPFR_MP_BITS (MANT_BITS (mpfr_digits ()))
  41  
  42  #define NO_MPFR ((mpfr_ptr) NULL)
  43  
  44  #define CHECK_MPFR(p, z) PRELUDE_ERROR (mpfr_number_p (z) == 0, (p), ERROR_MATH, M_LONG_LONG_REAL)
  45  
  46  void zeroin_mpfr (NODE_T *, mpfr_t *, mpfr_t, mpfr_t, mpfr_t, int (*)(mpfr_t, const mpfr_t, mpfr_rnd_t));
  47  
  48  //! @brief Decimal digits in mpfr significand.
  49  
  50  size_t mpfr_digits (void)
  51  {
  52    return (long_mp_digits () * LOG_MP_RADIX);
  53  }
  54  
  55  //! @brief Convert mp to mpfr.
  56  
  57  void mp_to_mpfr (NODE_T * p, MP_T * z, mpfr_t * x, int digits)
  58  {
  59  // This routine looks a lot like "strtod".
  60    (void) p;
  61    mpfr_set_ui (*x, 0, DEFAULT);
  62    if (MP_EXPONENT (z) * (MP_T) LOG_MP_RADIX > (MP_T) REAL_MIN_10_EXP) {
  63      int j, expo;
  64      BOOL_T neg = MP_DIGIT (z, 1) < 0;
  65      mpfr_t term, W;
  66      mpfr_inits2 (MPFR_MP_BITS, term, W, NO_MPFR);
  67      MP_DIGIT (z, 1) = ABS (MP_DIGIT (z, 1));
  68      expo = (int) (MP_EXPONENT (z) * LOG_MP_RADIX);
  69      mpfr_set_ui (W, 10, DEFAULT);
  70      mpfr_pow_si (W, W, expo, DEFAULT);
  71      for (j = 1; j <= digits; j++) {
  72        mpfr_set_d (term, MP_DIGIT (z, j), DEFAULT);
  73        mpfr_mul (term, term, W, DEFAULT);
  74        mpfr_add (*x, *x, term, DEFAULT);
  75        mpfr_div_ui (W, W, MP_RADIX, DEFAULT);
  76      }
  77      if (neg) {
  78        mpfr_neg (*x, *x, DEFAULT);
  79      }
  80    }
  81  }
  82  
  83  //! @brief Convert mpfr to mp number.
  84  
  85  MP_T *mpfr_to_mp (NODE_T * p, MP_T * z, mpfr_t * x, int digits)
  86  {
  87    mpfr_t u, v, t;
  88    SET_MP_ZERO (z, digits);
  89    if (mpfr_zero_p (*x)) {
  90      return z;
  91    }
  92    int sign_x = mpfr_sgn (*x);
  93    mpfr_inits2 (MPFR_MP_BITS, t, u, v, NO_MPFR);
  94  // Scale to [0, 0.1>.
  95  // a = ABS (x);
  96    mpfr_set (u, *x, DEFAULT);
  97    mpfr_abs (u, u, DEFAULT);
  98  // expo = (int) log10q (a);
  99    mpfr_log10 (v, u, DEFAULT);
 100    INT_T expo = mpfr_get_si (v, DEFAULT);
 101  // v /= ten_up_double (expo);
 102    mpfr_set_ui (v, 10, DEFAULT);
 103    mpfr_pow_si (v, v, expo, DEFAULT);
 104    mpfr_div (u, u, v, DEFAULT);
 105    expo--;
 106    if (mpfr_cmp_ui (u, 1) >= 0) {
 107      mpfr_div_ui (u, u, 10, DEFAULT);
 108      expo++;
 109    }
 110  // Transport digits of x to the mantissa of z.
 111    INT_T sum = 0, W = (MP_RADIX / 10); int j, k;
 112    for (k = 0, j = 1; j <= digits && k < mpfr_digits (); k++) {
 113      mpfr_mul_ui (t, u, 10, DEFAULT);
 114      mpfr_floor (v, t);
 115      mpfr_frac (u, t, DEFAULT);
 116      sum += W * mpfr_get_d (v, DEFAULT);
 117      W /= 10;
 118      if (W < 1) {
 119        MP_DIGIT (z, j++) = (MP_T) sum;
 120        sum = 0;
 121        W = (MP_RADIX / 10);
 122      }
 123    }
 124  // Store the last digits.
 125    if (j <= digits) {
 126      MP_DIGIT (z, j) = (MP_T) sum;
 127    }
 128    (void) align_mp (z, &expo, digits);
 129    MP_EXPONENT (z) = (MP_T) expo;
 130    MP_DIGIT (z, 1) *= sign_x;
 131    check_mp_exp (p, z);
 132    mpfr_clear (t);
 133    mpfr_clear (u);
 134    mpfr_clear (v);
 135    return z;
 136  }
 137  
 138  //! @brief PROC mpfr_mp = (LONG LONG REAL) LONG LONG REAL
 139  
 140  void genie_mpfr_mp (NODE_T * p)
 141  {
 142    MOID_T *mode = MOID (p);
 143    int digits = DIGITS (mode);
 144    int size = SIZE (mode);
 145    mpfr_t u;
 146    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 147    mpfr_init2 (u, MPFR_MP_BITS);
 148    mp_to_mpfr (p, z, &u, digits);
 149    mpfr_out_str (stdout, 10, 0, u, DEFAULT);
 150    CHECK_MPFR (p, u);
 151    mpfr_to_mp (p, z, &u, digits);
 152    mpfr_clear (u);
 153  }
 154  
 155  //! @brief mpfr_beta_inc
 156  
 157  void mpfr_beta_inc (mpfr_t I, mpfr_t s, mpfr_t t, mpfr_t x, mpfr_rnd_t rnd)
 158  {
 159  // Incomplete beta function I{x}(s, t).
 160  // From a continued fraction, see dlmf.nist.gov/8.17; Lentz's algorithm.
 161    errno = EDOM;                 // Until proven otherwise
 162  //mpfr_printf ("%.128Rf", x);
 163    if (mpfr_cmp_d (x, 0) < 0 || mpfr_cmp_d (x, 1) > 0) {
 164      mpfr_set_nan (I);
 165    } else {
 166      mpfr_t a, b, c, d, e, F, T, W;
 167      int N, m, cont = A68_TRUE;
 168      mpfr_prec_t lim = 2 * mpfr_get_prec (x);
 169      mpfr_inits2 (MPFR_MP_BITS, a, b, c, d, e, F, T, W, NO_MPFR);
 170  // Rapid convergence when x < (s+1)/(s+t+2)
 171      mpfr_add_d (a, s, 1, rnd);
 172      mpfr_add (b, s, t, rnd);
 173      mpfr_add_d (b, b, 2, rnd);
 174      mpfr_div (c, a, b, rnd);
 175  // Recursion when x > (s+1)/(s+t+2)
 176      if (mpfr_cmp (x, c) > 0) {
 177  // B{x}(s, t) = 1 - B{1-x}(t, s)
 178        mpfr_d_sub (d, 1, x, rnd);
 179        mpfr_beta_inc (I, t, s, d, rnd);
 180        mpfr_d_sub (I, 1, I, rnd);
 181        mpfr_clears (a, b, c, d, e, F, T, W, NO_MPFR);
 182        return;
 183      }
 184  // Lentz's algorithm for continued fraction.
 185      mpfr_set_d (W, 1, rnd);
 186      mpfr_set_d (F, 1, rnd);
 187      mpfr_set_d (c, 1, rnd);
 188      mpfr_set_d (d, 0, rnd);
 189      for (N = 0, m = 0; cont && N < lim; N++) {
 190        if (N == 0) {
 191  // d := 1
 192          mpfr_set_d (T, 1, rnd);
 193        } else if (N % 2 == 0) {
 194  // d{2m} := x m(t-m)/((s+2m-1)(s+2m))
 195          mpfr_sub_si (a, t, m, rnd);
 196          mpfr_mul_si (a, a, m, rnd);
 197          mpfr_mul (a, a, x, rnd);
 198          mpfr_add_si (b, s, m, rnd);
 199          mpfr_add_si (b, b, m, rnd);
 200          mpfr_set (e, b, rnd);
 201          mpfr_sub_d (b, b, 1, rnd);
 202          mpfr_mul (b, b, e, rnd);
 203          mpfr_div (T, a, b, rnd);
 204        } else {
 205  // d{2m+1} := -x (s+m)(s+t+m)/((s+2m+1)(s+2m))
 206          mpfr_add_si (e, s, m, rnd);
 207          mpfr_add (T, e, t, rnd);
 208          mpfr_mul (a, e, T, rnd);
 209          mpfr_mul (a, a, x, rnd);
 210          mpfr_add_si (b, s, m, rnd);
 211          mpfr_add_si (b, b, m, rnd);
 212          mpfr_set (e, b, rnd);
 213          mpfr_add_d (b, b, 1, rnd);
 214          mpfr_mul (b, b, e, rnd);
 215          mpfr_div (T, a, b, rnd);
 216          mpfr_neg (T, T, rnd);
 217          m++;
 218        }
 219        mpfr_mul (e, T, d, rnd);
 220        mpfr_add_d (d, e, 1, rnd);
 221        mpfr_d_div (d, 1, d, rnd);
 222        mpfr_div (e, T, c, rnd);
 223        mpfr_add_d (c, e, 1, rnd);
 224        mpfr_mul (F, F, c, rnd);
 225        mpfr_mul (F, F, d, rnd);
 226        if (mpfr_cmp (F, W) == 0) {
 227          cont = A68_FALSE;
 228          errno = 0;
 229        } else {
 230          mpfr_set (W, F, rnd);
 231        }
 232      }
 233  // I{x}(s,t)=x^s(1-x)^t / a / B(s,t) F
 234      mpfr_pow (a, x, s, rnd);
 235      mpfr_d_sub (b, 1, x, rnd);
 236      mpfr_pow (b, b, t, rnd);
 237      mpfr_mul (a, a, b, rnd);
 238      mpfr_beta (W, s, t, rnd);
 239      mpfr_sub_d (F, F, 1, rnd);
 240      mpfr_mul (b, F, a, rnd);
 241      mpfr_div (b, b, W, rnd);
 242      mpfr_div (b, b, s, rnd);
 243      mpfr_set (I, b, rnd);
 244      mpfr_clears (a, b, c, d, e, F, T, W, NO_MPFR);
 245    }
 246  }
 247  
 248  //! @brief PROC long long erf = (LONG LONG REAL) LONG LONG REAL
 249  
 250  void genie_mpfr_erf_mp (NODE_T * p)
 251  {
 252    MOID_T *mode = MOID (p);
 253    int digits = DIGITS (mode);
 254    int size = SIZE (mode);
 255    mpfr_t u;
 256    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 257    mpfr_init2 (u, MPFR_MP_BITS);
 258    mp_to_mpfr (p, z, &u, digits);
 259    mpfr_erf (u, u, DEFAULT);
 260    CHECK_MPFR (p, u);
 261    mpfr_to_mp (p, z, &u, digits);
 262    mpfr_clear (u);
 263  }
 264  
 265  //! @brief PROC long long erfc = (LONG LONG REAL) LONG LONG REAL
 266  
 267  void genie_mpfr_erfc_mp (NODE_T * p)
 268  {
 269    MOID_T *mode = MOID (p);
 270    int digits = DIGITS (mode);
 271    int size = SIZE (mode);
 272    mpfr_t u;
 273    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 274    mpfr_init2 (u, MPFR_MP_BITS);
 275    mp_to_mpfr (p, z, &u, digits);
 276    mpfr_erfc (u, u, DEFAULT);
 277    CHECK_MPFR (p, u);
 278    mpfr_to_mp (p, z, &u, digits);
 279    mpfr_clear (u);
 280  }
 281  
 282  //! @brief PROC long long inverf = (LONG LONG REAL) LONG LONG REAL
 283  
 284  void genie_mpfr_inverf_mp (NODE_T * _p_)
 285  {
 286    MOID_T *mode = MOID (_p_);
 287    int digits = DIGITS (mode), size = SIZE (mode);
 288    REAL_T x0;
 289    mpfr_t a, b, u, y;
 290    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 291    A68 (f_entry) = _p_;
 292    mpfr_inits2 (MPFR_MP_BITS, a, b, u, y, NO_MPFR);
 293    mp_to_mpfr (_p_, z, &y, digits);
 294    x0 = a68_inverf (mp_to_real (_p_, z, digits));
 295  //  a = z - 1e-9;
 296  //  b = z + 1e-9;
 297    mpfr_set_d (a, x0 - 1e-9, DEFAULT);
 298    mpfr_set_d (b, x0 + 1e-9, DEFAULT);
 299    zeroin_mpfr (_p_, &u, a, b, y, mpfr_erf);
 300    MATH_RTE (_p_, errno != 0, M_LONG_LONG_REAL, NO_TEXT);
 301    mpfr_to_mp (_p_, z, &u, digits);
 302    mpfr_clears (a, b, u, y, NO_MPFR);
 303  }
 304  
 305  //! @brief PROC long long inverfc = (LONG LONG REAL) LONG LONG REAL
 306  
 307  void genie_mpfr_inverfc_mp (NODE_T * p)
 308  {
 309    MOID_T *mode = MOID (p);
 310    ADDR_T pop_sp = A68_SP;
 311    int digits = DIGITS (mode), size = SIZE (mode);
 312    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 313    one_minus_mp (p, z, z, digits);
 314    A68_SP = pop_sp;
 315    genie_inverf_mp (p);
 316  }
 317  
 318  //! @brief PROC long long gamma = (LONG LONG REAL) LONG LONG REAL
 319  
 320  void genie_gamma_mpfr (NODE_T * p)
 321  {
 322    MOID_T *mode = MOID (p);
 323    int digits = DIGITS (mode);
 324    int size = SIZE (mode);
 325    mpfr_t u;
 326    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 327    mpfr_init2 (u, MPFR_MP_BITS);
 328    mp_to_mpfr (p, z, &u, digits);
 329    mpfr_gamma (u, u, DEFAULT);
 330    CHECK_MPFR (p, u);
 331    mpfr_to_mp (p, z, &u, digits);
 332    mpfr_clear (u);
 333  }
 334  
 335  //! @brief PROC long long ln gamma = (LONG LONG REAL) LONG LONG REAL
 336  
 337  void genie_lngamma_mpfr (NODE_T * p)
 338  {
 339    MOID_T *mode = MOID (p);
 340    int digits = DIGITS (mode);
 341    int size = SIZE (mode);
 342    mpfr_t u;
 343    MP_T *z = (MP_T *) STACK_OFFSET (-size);
 344    mpfr_init2 (u, MPFR_MP_BITS);
 345    mp_to_mpfr (p, z, &u, digits);
 346    mpfr_lngamma (u, u, DEFAULT);
 347    CHECK_MPFR (p, u);
 348    mpfr_to_mp (p, z, &u, digits);
 349    mpfr_clear (u);
 350  }
 351  
 352  void genie_gamma_inc_real_mpfr (NODE_T * p)
 353  {
 354    A68_REAL x, s;
 355    mpfr_t xx, ss;
 356    POP_OBJECT (p, &x, A68_REAL);
 357    POP_OBJECT (p, &s, A68_REAL);
 358    mpfr_inits2 (MPFR_LONG_REAL_BITS, ss, xx, NO_MPFR);
 359    mpfr_set_d (xx, VALUE (&x), DEFAULT);
 360    mpfr_set_d (ss, VALUE (&s), DEFAULT);
 361    mpfr_gamma_inc (ss, ss, xx, DEFAULT);
 362    CHECK_MPFR (p, ss);
 363    PUSH_VALUE (p, mpfr_get_d (ss, DEFAULT), A68_REAL);
 364    mpfr_clears (ss, xx, NO_MPFR);
 365  }
 366  
 367  void genie_gamma_inc_double_real_mpfr (NODE_T * p)
 368  {
 369    A68_LONG_REAL x, s;
 370    mpfr_t xx, ss;
 371    POP_OBJECT (p, &x, A68_LONG_REAL);
 372    POP_OBJECT (p, &s, A68_LONG_REAL);
 373    mpfr_inits2 (MPFR_LONG_REAL_BITS, ss, xx, NO_MPFR);
 374    mpfr_set_float128 (xx, VALUE (&x).f, DEFAULT);
 375    mpfr_set_float128 (ss, VALUE (&s).f, DEFAULT);
 376    mpfr_gamma_inc (ss, ss, xx, DEFAULT);
 377    CHECK_MPFR (p, ss);
 378    PUSH_VALUE (p, dble (mpfr_get_float128 (ss, DEFAULT)), A68_LONG_REAL);
 379    mpfr_clears (ss, xx, NO_MPFR);
 380  }
 381  
 382  void genie_gamma_inc_mpfr (NODE_T * p)
 383  {
 384    int digits = DIGITS (MOID (p)), size = SIZE (MOID (p));
 385    MP_T *x = (MP_T *) STACK_OFFSET (-size);
 386    MP_T *s = (MP_T *) STACK_OFFSET (-2 * size);
 387    mpfr_t xx, ss;
 388    A68_SP -= size;
 389    mpfr_inits2 (MPFR_MP_BITS, ss, xx, NO_MPFR);
 390    mp_to_mpfr (p, x, &xx, digits);
 391    mp_to_mpfr (p, s, &ss, digits);
 392    mpfr_gamma_inc (ss, ss, xx, DEFAULT);
 393    CHECK_MPFR (p, ss);
 394    mpfr_to_mp (p, s, &ss, digits);
 395    mpfr_clears (ss, xx, NO_MPFR);
 396  }
 397  
 398  void genie_beta_mpfr (NODE_T * p)
 399  {
 400    int digits = DIGITS (MOID (p)), size = SIZE (MOID (p));
 401    MP_T *x = (MP_T *) STACK_OFFSET (-size);
 402    MP_T *s = (MP_T *) STACK_OFFSET (-2 * size);
 403    mpfr_t xx, ss;
 404    A68_SP -= size;
 405    mpfr_inits2 (MPFR_MP_BITS, ss, xx, NO_MPFR);
 406    mp_to_mpfr (p, x, &xx, digits);
 407    mp_to_mpfr (p, s, &ss, digits);
 408    mpfr_beta (ss, ss, xx, DEFAULT);
 409    CHECK_MPFR (p, ss);
 410    mpfr_to_mp (p, s, &ss, digits);
 411    mpfr_clears (ss, xx, NO_MPFR);
 412  }
 413  
 414  void genie_ln_beta_mpfr (NODE_T * p)
 415  {
 416    int digits = DIGITS (MOID (p)), size = SIZE (MOID (p));
 417    MP_T *b = (MP_T *) STACK_OFFSET (-size);
 418    MP_T *a = (MP_T *) STACK_OFFSET (-2 * size);
 419    mpfr_t bb, aa, yy, zz;
 420    A68_SP -= size;
 421    mpfr_inits2 (MPFR_MP_BITS, aa, bb, yy, zz, NO_MPFR);
 422    mp_to_mpfr (p, b, &bb, digits);
 423    mp_to_mpfr (p, a, &aa, digits);
 424    mpfr_lngamma (zz, aa, DEFAULT);
 425    mpfr_lngamma (yy, bb, DEFAULT);
 426    mpfr_add (zz, zz, yy, DEFAULT);
 427    mpfr_add (yy, aa, bb, DEFAULT);
 428    mpfr_lngamma (yy, yy, DEFAULT);
 429    mpfr_sub (aa, zz, yy, DEFAULT);
 430    CHECK_MPFR (p, aa);
 431    mpfr_to_mp (p, a, &aa, digits);
 432    mpfr_clears (aa, bb, yy, zz, NO_MPFR);
 433  }
 434  
 435  void genie_beta_inc_mpfr (NODE_T * p)
 436  {
 437    int digits = DIGITS (MOID (p)), size = SIZE (MOID (p));
 438    MP_T *x = (MP_T *) STACK_OFFSET (-size);
 439    MP_T *t = (MP_T *) STACK_OFFSET (-2 * size);
 440    MP_T *s = (MP_T *) STACK_OFFSET (-3 * size);
 441    mpfr_t xx, ss, tt;
 442    A68_SP -= 2 * size;
 443    mpfr_inits2 (MPFR_MP_BITS, ss, tt, xx, NO_MPFR);
 444    mp_to_mpfr (p, x, &xx, digits);
 445    mp_to_mpfr (p, s, &ss, digits);
 446    mp_to_mpfr (p, t, &tt, digits);
 447    mpfr_beta_inc (ss, ss, tt, xx, DEFAULT);
 448    CHECK_MPFR (p, ss);
 449    mpfr_to_mp (p, s, &ss, digits);
 450    mpfr_clears (ss, tt, xx, NO_MPFR);
 451  }
 452  
 453  //! @brief zeroin
 454  
 455  void zeroin_mpfr (NODE_T * _p_, mpfr_t * z, mpfr_t a, mpfr_t b, mpfr_t y, int (*f) (mpfr_t, const mpfr_t, mpfr_rnd_t))
 456  {
 457  // 'zeroin'
 458  // MCA 2310 in 'ALGOL 60 Procedures in Numerical Algebra' by Th.J. Dekker
 459    int sign, its = 5;
 460    BOOL_T go_on = A68_TRUE;
 461    mpfr_t c, fa, fb, fc, tolb, eps, p, q, v, w;
 462    mpfr_inits2 (MPFR_MP_BITS, c, fa, fb, fc, tolb, eps, p, q, v, w, NO_MPFR);
 463    mpfr_set_ui (eps, 10, DEFAULT);
 464    mpfr_pow_si (eps, eps, -(mpfr_digits () - 2), DEFAULT);
 465    f (fa, a, DEFAULT);
 466    mpfr_sub (fa, fa, y, DEFAULT);
 467    f (fb, b, DEFAULT);
 468    mpfr_sub (fb, fb, y, DEFAULT);
 469    mpfr_set (c, a, DEFAULT);
 470    mpfr_set (fc, fa, DEFAULT);
 471    while (go_on && (its--) > 0) {
 472      mpfr_abs (v, fc, DEFAULT);
 473      mpfr_abs (w, fb, DEFAULT);
 474      if (mpfr_cmp (v, w) < 0) {
 475        mpfr_set (a, b, DEFAULT);
 476        mpfr_set (fa, fb, DEFAULT);
 477        mpfr_set (b, c, DEFAULT);
 478        mpfr_set (fb, fc, DEFAULT);
 479        mpfr_set (c, a, DEFAULT);
 480        mpfr_set (fc, fa, DEFAULT);
 481      }
 482      mpfr_abs (tolb, b, DEFAULT);
 483      mpfr_add_ui (tolb, tolb, 1, DEFAULT);
 484      mpfr_mul (tolb, tolb, eps, DEFAULT);
 485      mpfr_add (w, c, b, DEFAULT);
 486      mpfr_div_2ui (w, w, 1, DEFAULT);
 487      mpfr_sub (v, w, b, DEFAULT);
 488      mpfr_abs (v, v, DEFAULT);
 489      go_on = mpfr_cmp (v, tolb) > 0;
 490      if (go_on) {
 491        mpfr_sub (p, b, a, DEFAULT);
 492        mpfr_mul (p, p, fb, DEFAULT);
 493        mpfr_sub (q, fa, fb, DEFAULT);
 494        if (mpfr_cmp_ui (p, 0) < 0) {
 495          mpfr_neg (p, p, DEFAULT);
 496          mpfr_neg (q, q, DEFAULT);
 497        }
 498        mpfr_set (a, b, DEFAULT);
 499        mpfr_set (fa, fb, DEFAULT);
 500        mpfr_abs (v, q, DEFAULT);
 501        mpfr_mul (v, v, tolb, DEFAULT);
 502        if (mpfr_cmp (p, v) <= 0) {
 503          if (mpfr_cmp (c, b) > 0) {
 504            mpfr_add (b, b, tolb, DEFAULT);
 505          } else {
 506            mpfr_sub (b, b, tolb, DEFAULT);
 507          }
 508        } else {
 509          mpfr_sub (v, w, b, DEFAULT);
 510          mpfr_mul (v, v, q, DEFAULT);
 511          if (mpfr_cmp (p, v) < 0) {
 512            mpfr_div (v, p, q, DEFAULT);
 513            mpfr_add (b, v, b, DEFAULT);
 514          } else {
 515            mpfr_set (b, w, DEFAULT);
 516          }
 517        }
 518        f (fb, b, DEFAULT);
 519        mpfr_sub (fb, fb, y, DEFAULT);
 520        sign = mpfr_sgn (fb) + mpfr_sgn (fc);
 521        if (ABS (sign) == 2) {
 522          mpfr_set (c, a, DEFAULT);
 523          mpfr_set (fc, fa, DEFAULT);
 524        }
 525      }
 526    }
 527    CHECK_MPFR (_p_, b);
 528    mpfr_set (*z, b, DEFAULT);
 529    mpfr_clears (c, fa, fb, fc, tolb, eps, p, q, v, w, NO_MPFR);
 530  }
 531  
 532  #endif
 533  
 534  #endif