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
