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