1 /* Exception flags and utilities. Constructors and destructors (debug).
2
3 Copyright 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-impl.h"
24
25 MPFR_THREAD_VAR (mpfr_flags_t, __gmpfr_flags, 0)
26 MPFR_THREAD_VAR (mpfr_exp_t, __gmpfr_emin, MPFR_EMIN_DEFAULT)
27 MPFR_THREAD_VAR (mpfr_exp_t, __gmpfr_emax, MPFR_EMAX_DEFAULT)
28
29 #undef mpfr_get_emin
30
31 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
32 mpfr_get_emin (void)
33 {
34 return __gmpfr_emin;
35 }
36
37 #undef mpfr_set_emin
38
39 int
40 mpfr_set_emin (mpfr_exp_t exponent)
41 {
42 if (MPFR_LIKELY (exponent >= MPFR_EMIN_MIN && exponent <= MPFR_EMIN_MAX))
43 {
44 __gmpfr_emin = exponent;
45 return 0;
46 }
47 else
48 {
49 return 1;
50 }
51 }
52
53 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
54 mpfr_get_emin_min (void)
55 {
56 return MPFR_EMIN_MIN;
57 }
58
59 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
60 mpfr_get_emin_max (void)
61 {
62 return MPFR_EMIN_MAX;
63 }
64
65 #undef mpfr_get_emax
66
67 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
68 mpfr_get_emax (void)
69 {
70 return __gmpfr_emax;
71 }
72
73 #undef mpfr_set_emax
74
75 int
76 mpfr_set_emax (mpfr_exp_t exponent)
77 {
78 if (MPFR_LIKELY (exponent >= MPFR_EMAX_MIN && exponent <= MPFR_EMAX_MAX))
79 {
80 __gmpfr_emax = exponent;
81 return 0;
82 }
83 else
84 {
85 return 1;
86 }
87 }
88
89 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
90 mpfr_get_emax_min (void)
91 {
92 return MPFR_EMAX_MIN;
93 }
94
95 MPFR_COLD_FUNCTION_ATTR mpfr_exp_t
96 mpfr_get_emax_max (void)
97 {
98 return MPFR_EMAX_MAX;
99 }
100
101
102 #undef mpfr_flags_clear
103
104 MPFR_COLD_FUNCTION_ATTR void
105 mpfr_flags_clear (mpfr_flags_t mask)
106 {
107 __gmpfr_flags &= MPFR_FLAGS_ALL ^ mask;
108 }
109
110 #undef mpfr_flags_set
111
112 MPFR_COLD_FUNCTION_ATTR void
113 mpfr_flags_set (mpfr_flags_t mask)
114 {
115 __gmpfr_flags |= mask;
116 }
117
118 #undef mpfr_flags_test
119
120 MPFR_COLD_FUNCTION_ATTR mpfr_flags_t
121 mpfr_flags_test (mpfr_flags_t mask)
122 {
123 return __gmpfr_flags & mask;
124 }
125
126 #undef mpfr_flags_save
127
128 MPFR_COLD_FUNCTION_ATTR mpfr_flags_t
129 mpfr_flags_save (void)
130 {
131 return __gmpfr_flags;
132 }
133
134 #undef mpfr_flags_restore
135
136 MPFR_COLD_FUNCTION_ATTR void
137 mpfr_flags_restore (mpfr_flags_t flags, mpfr_flags_t mask)
138 {
139 __gmpfr_flags =
140 (__gmpfr_flags & (MPFR_FLAGS_ALL ^ mask)) |
141 (flags & mask);
142 }
143
144
145 #undef mpfr_clear_flags
146
147 void
148 mpfr_clear_flags (void)
149 {
150 __gmpfr_flags = 0;
151 }
152
153 #undef mpfr_clear_underflow
154
155 MPFR_COLD_FUNCTION_ATTR void
156 mpfr_clear_underflow (void)
157 {
158 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_UNDERFLOW;
159 }
160
161 #undef mpfr_clear_overflow
162
163 MPFR_COLD_FUNCTION_ATTR void
164 mpfr_clear_overflow (void)
165 {
166 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_OVERFLOW;
167 }
168
169 #undef mpfr_clear_divby0
170
171 MPFR_COLD_FUNCTION_ATTR void
172 mpfr_clear_divby0 (void)
173 {
174 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_DIVBY0;
175 }
176
177 #undef mpfr_clear_nanflag
178
179 MPFR_COLD_FUNCTION_ATTR void
180 mpfr_clear_nanflag (void)
181 {
182 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_NAN;
183 }
184
185 #undef mpfr_clear_inexflag
186
187 MPFR_COLD_FUNCTION_ATTR void
188 mpfr_clear_inexflag (void)
189 {
190 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT;
191 }
192
193 #undef mpfr_clear_erangeflag
194
195 MPFR_COLD_FUNCTION_ATTR void
196 mpfr_clear_erangeflag (void)
197 {
198 __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE;
199 }
200
201 #undef mpfr_set_underflow
202
203 MPFR_COLD_FUNCTION_ATTR void
204 mpfr_set_underflow (void)
205 {
206 __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW;
207 }
208
209 #undef mpfr_set_overflow
210
211 MPFR_COLD_FUNCTION_ATTR void
212 mpfr_set_overflow (void)
213 {
214 __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
215 }
216
217 #undef mpfr_set_divby0
218
219 MPFR_COLD_FUNCTION_ATTR void
220 mpfr_set_divby0 (void)
221 {
222 __gmpfr_flags |= MPFR_FLAGS_DIVBY0;
223 }
224
225 #undef mpfr_set_nanflag
226
227 MPFR_COLD_FUNCTION_ATTR void
228 mpfr_set_nanflag (void)
229 {
230 __gmpfr_flags |= MPFR_FLAGS_NAN;
231 }
232
233 #undef mpfr_set_inexflag
234
235 MPFR_COLD_FUNCTION_ATTR void
236 mpfr_set_inexflag (void)
237 {
238 __gmpfr_flags |= MPFR_FLAGS_INEXACT;
239 }
240
241 #undef mpfr_set_erangeflag
242
243 MPFR_COLD_FUNCTION_ATTR void
244 mpfr_set_erangeflag (void)
245 {
246 __gmpfr_flags |= MPFR_FLAGS_ERANGE;
247 }
248
249
250 #undef mpfr_check_range
251
252 /* Note: It is possible that for pure FP numbers, EXP(x) < MPFR_EMIN_MIN,
253 but the caller must make sure that the difference remains small enough
254 to avoid reaching the special exponent values. */
255 /* This function does not have logging messages. As it is also partly
256 implemented as a macro, if messages are added in the future, the macro
257 may need to be disabled when logging is enabled. */
258 int
259 mpfr_check_range (mpfr_ptr x, int t, mpfr_rnd_t rnd_mode)
260 {
261 if (MPFR_LIKELY (! MPFR_IS_SINGULAR (x)))
262 { /* x is a non-zero FP */
263 mpfr_exp_t exp = MPFR_EXP (x); /* Do not use MPFR_GET_EXP */
264
265 MPFR_ASSERTD (MPFR_IS_NORMALIZED (x));
266 if (MPFR_UNLIKELY (exp < __gmpfr_emin))
267 {
268 /* The following test is necessary because in the rounding to the
269 * nearest mode, mpfr_underflow always rounds away from 0. In
270 * this rounding mode, we need to round to 0 if:
271 * _ |x| < 2^(emin-2), or
272 * _ |x| = 2^(emin-2) and the absolute value of the exact
273 * result is <= 2^(emin-2).
274 */
275 if (rnd_mode == MPFR_RNDN &&
276 (exp + 1 < __gmpfr_emin ||
277 (mpfr_powerof2_raw(x) &&
278 (MPFR_IS_NEG(x) ? t <= 0 : t >= 0))))
279 rnd_mode = MPFR_RNDZ;
280 return mpfr_underflow (x, rnd_mode, MPFR_SIGN(x));
281 }
282 if (MPFR_UNLIKELY (exp > __gmpfr_emax))
283 return mpfr_overflow (x, rnd_mode, MPFR_SIGN(x));
284 }
285 else if (MPFR_UNLIKELY (t != 0 && MPFR_IS_INF (x)))
286 {
287 /* We need to do the following because most MPFR functions are
288 * implemented in the following way:
289 * Ziv's loop:
290 * | Compute an approximation to the result and an error bound.
291 * | Possible underflow/overflow detection -> return.
292 * | If can_round, break (exit the loop).
293 * | Otherwise, increase the working precision and loop.
294 * Round the approximation in the target precision. <== See below
295 * Restore the flags (that could have been set due to underflows
296 * or overflows during the internal computations).
297 * Execute: return mpfr_check_range (...).
298 * The problem is that an overflow could be generated when rounding the
299 * approximation (in general, such an overflow could not be detected
300 * earlier), and the overflow flag is lost when the flags are restored.
301 * This can occur only when the rounding yields an exponent change
302 * and the new exponent is larger than the maximum exponent, so that
303 * an infinity is necessarily obtained.
304 * So, the simplest solution is to detect this overflow case here in
305 * mpfr_check_range, which is easy to do since the rounded result is
306 * necessarily an inexact infinity.
307 */
308 __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
309 }
310 MPFR_RET (t); /* propagate inexact ternary value, unlike most functions */
311 }
312
313
314 #undef mpfr_underflow_p
315
316 MPFR_COLD_FUNCTION_ATTR int
317 mpfr_underflow_p (void)
318 {
319 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_UNDERFLOW <= INT_MAX);
320 return __gmpfr_flags & MPFR_FLAGS_UNDERFLOW;
321 }
322
323 #undef mpfr_overflow_p
324
325 MPFR_COLD_FUNCTION_ATTR int
326 mpfr_overflow_p (void)
327 {
328 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_OVERFLOW <= INT_MAX);
329 return __gmpfr_flags & MPFR_FLAGS_OVERFLOW;
330 }
331
332 #undef mpfr_divby0_p
333
334 MPFR_COLD_FUNCTION_ATTR int
335 mpfr_divby0_p (void)
336 {
337 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_DIVBY0 <= INT_MAX);
338 return __gmpfr_flags & MPFR_FLAGS_DIVBY0;
339 }
340
341 #undef mpfr_nanflag_p
342
343 MPFR_COLD_FUNCTION_ATTR int
344 mpfr_nanflag_p (void)
345 {
346 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_NAN <= INT_MAX);
347 return __gmpfr_flags & MPFR_FLAGS_NAN;
348 }
349
350 #undef mpfr_inexflag_p
351
352 MPFR_COLD_FUNCTION_ATTR int
353 mpfr_inexflag_p (void)
354 {
355 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_INEXACT <= INT_MAX);
356 return __gmpfr_flags & MPFR_FLAGS_INEXACT;
357 }
358
359 #undef mpfr_erangeflag_p
360
361 MPFR_COLD_FUNCTION_ATTR int
362 mpfr_erangeflag_p (void)
363 {
364 MPFR_STAT_STATIC_ASSERT (MPFR_FLAGS_ERANGE <= INT_MAX);
365 return __gmpfr_flags & MPFR_FLAGS_ERANGE;
366 }
367
368
369 /* #undef mpfr_underflow */
370
371 /* Note: In the rounding to the nearest mode, mpfr_underflow
372 always rounds away from 0. In this rounding mode, you must call
373 mpfr_underflow with rnd_mode = MPFR_RNDZ if the exact result
374 is <= 2^(emin-2) in absolute value.
375 We chose the default to round away from zero instead of toward zero
376 because rounding away from zero (MPFR_RNDA) wasn't supported at that
377 time (r1910), so that the caller had no way to change rnd_mode to
378 this mode. */
379
380 MPFR_COLD_FUNCTION_ATTR int
381 mpfr_underflow (mpfr_ptr x, mpfr_rnd_t rnd_mode, int sign)
382 {
383 int inex;
384
385 MPFR_LOG_FUNC
386 (("rnd=%d sign=%d", rnd_mode, sign),
387 ("x[%Pd]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x));
388
389 MPFR_ASSERT_SIGN (sign);
390
391 if (MPFR_IS_LIKE_RNDZ(rnd_mode, sign < 0))
392 {
393 MPFR_SET_ZERO(x);
394 inex = -1;
395 }
396 else
397 {
398 mpfr_setmin (x, __gmpfr_emin);
399 inex = 1;
400 }
401 MPFR_SET_SIGN(x, sign);
402 __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
403 return sign > 0 ? inex : -inex;
404 }
405
406 /* #undef mpfr_overflow */
407
408 MPFR_COLD_FUNCTION_ATTR int
409 mpfr_overflow (mpfr_ptr x, mpfr_rnd_t rnd_mode, int sign)
410 {
411 int inex;
412
413 MPFR_LOG_FUNC
414 (("rnd=%d sign=%d", rnd_mode, sign),
415 ("x[%Pd]=%.*Rg", mpfr_get_prec (x), mpfr_log_prec, x));
416
417 MPFR_ASSERT_SIGN (sign);
418
419 if (MPFR_IS_LIKE_RNDZ(rnd_mode, sign < 0))
420 {
421 mpfr_setmax (x, __gmpfr_emax);
422 inex = -1;
423 }
424 else
425 {
426 MPFR_SET_INF(x);
427 inex = 1;
428 }
429 MPFR_SET_SIGN(x, sign);
430 __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
431 return sign > 0 ? inex : -inex;
432 }
433
434 /**************************************************************************/
435
436 /* Code related to constructors and destructors (for debugging) should
437 be put here. The reason is that such code must be in an object file
438 that will be kept by the linker for symbol resolution, and symbols
439 __gmpfr_emin and __gmpfr_emax from this file will be used by every
440 program calling a MPFR math function (where rounding is involved). */
441
442 #if defined MPFR_DEBUG_PREDICTION
443
444 /* Print prediction statistics at the end of a program.
445 *
446 * Code to debug branch prediction, based on Ulrich Drepper's paper
447 * "What Every Programmer Should Know About Memory":
448 * https://people.freebsd.org/~lstewart/articles/cpumemory.pdf
449 */
450
451 extern long int __start_predict_data;
452 extern long int __stop_predict_data;
453 extern long int __start_predict_line;
454 extern const char *__start_predict_file;
455
456 static void __attribute__ ((destructor))
457 predprint (void)
458 {
459 long int *s = &__start_predict_data;
460 long int *e = &__stop_predict_data;
461 long int *sl = &__start_predict_line;
462 const char **sf = &__start_predict_file;
463
464 while (s < e)
465 {
466 printf("%s:%ld: incorrect=%ld, correct=%ld%s\n",
467 *sf, *sl, s[0], s[1],
468 s[0] > s[1] ? " <==== WARNING" : "");
469 ++sl;
470 ++sf;
471 s += 2;
472 }
473 }
474
475 #endif
476
477 #if MPFR_WANT_ASSERT >= 2
478
479 /* Similar to flags_out in tests/tests.c */
480
481 void
482 flags_fout (FILE *stream, mpfr_flags_t flags)
483 {
484 int none = 1;
485
486 if (flags & MPFR_FLAGS_UNDERFLOW)
487 none = 0, fprintf (stream, " underflow");
488 if (flags & MPFR_FLAGS_OVERFLOW)
489 none = 0, fprintf (stream, " overflow");
490 if (flags & MPFR_FLAGS_NAN)
491 none = 0, fprintf (stream, " nan");
492 if (flags & MPFR_FLAGS_INEXACT)
493 none = 0, fprintf (stream, " inexact");
494 if (flags & MPFR_FLAGS_ERANGE)
495 none = 0, fprintf (stream, " erange");
496 if (none)
497 fprintf (stream, " none");
498 fprintf (stream, " (%u)\n", flags);
499 }
500
501 #endif