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