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