contrib/mpfr/src/exceptions.c

   1 /* Exception flags and utilities.
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel 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 http://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 unsigned int MPFR_THREAD_ATTR __gmpfr_flags = 0;
  26
  27 mpfr_exp_t MPFR_THREAD_ATTR __gmpfr_emin = MPFR_EMIN_DEFAULT;
  28 mpfr_exp_t MPFR_THREAD_ATTR __gmpfr_emax = MPFR_EMAX_DEFAULT;
  29
  30 #undef mpfr_get_emin
  31
  32 mpfr_exp_t
  33 mpfr_get_emin (void)
  34 {
  35   return __gmpfr_emin;
  36 }
  37
  38 #undef mpfr_set_emin
  39
  40 int
  41 mpfr_set_emin (mpfr_exp_t exponent)
  42 {
  43   if (exponent >= MPFR_EMIN_MIN && exponent <= MPFR_EMIN_MAX)
  44     {
  45       __gmpfr_emin = exponent;
  46       return 0;
  47     }
  48   else
  49     {
  50       return 1;
  51     }
  52 }
  53
  54 mpfr_exp_t
  55 mpfr_get_emin_min (void)
  56 {
  57   return MPFR_EMIN_MIN;
  58 }
  59
  60 mpfr_exp_t
  61 mpfr_get_emin_max (void)
  62 {
  63   return MPFR_EMIN_MAX;
  64 }
  65
  66 #undef mpfr_get_emax
  67
  68 mpfr_exp_t
  69 mpfr_get_emax (void)
  70 {
  71   return __gmpfr_emax;
  72 }
  73
  74 #undef mpfr_set_emax
  75
  76 int
  77 mpfr_set_emax (mpfr_exp_t exponent)
  78 {
  79   if (exponent >= MPFR_EMAX_MIN && exponent <= MPFR_EMAX_MAX)
  80     {
  81       __gmpfr_emax = exponent;
  82       return 0;
  83     }
  84   else
  85     {
  86       return 1;
  87     }
  88 }
  89
  90 mpfr_exp_t
  91 mpfr_get_emax_min (void)
  92 {
  93   return MPFR_EMAX_MIN;
  94 }
  95 mpfr_exp_t
  96 mpfr_get_emax_max (void)
  97 {
  98   return MPFR_EMAX_MAX;
  99 }
 100
 101
 102 #undef mpfr_clear_flags
 103
 104 void
 105 mpfr_clear_flags (void)
 106 {
 107   __gmpfr_flags = 0;
 108 }
 109
 110 #undef mpfr_clear_underflow
 111
 112 void
 113 mpfr_clear_underflow (void)
 114 {
 115   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_UNDERFLOW;
 116 }
 117
 118 #undef mpfr_clear_overflow
 119
 120 void
 121 mpfr_clear_overflow (void)
 122 {
 123   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_OVERFLOW;
 124 }
 125
 126 #undef mpfr_clear_divby0
 127
 128 void
 129 mpfr_clear_divby0 (void)
 130 {
 131   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_DIVBY0;
 132 }
 133
 134 #undef mpfr_clear_nanflag
 135
 136 void
 137 mpfr_clear_nanflag (void)
 138 {
 139   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_NAN;
 140 }
 141
 142 #undef mpfr_clear_inexflag
 143
 144 void
 145 mpfr_clear_inexflag (void)
 146 {
 147   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT;
 148 }
 149
 150 #undef mpfr_clear_erangeflag
 151
 152 void
 153 mpfr_clear_erangeflag (void)
 154 {
 155   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE;
 156 }
 157
 158 #undef mpfr_set_underflow
 159
 160 void
 161 mpfr_set_underflow (void)
 162 {
 163   __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW;
 164 }
 165
 166 #undef mpfr_set_overflow
 167
 168 void
 169 mpfr_set_overflow (void)
 170 {
 171   __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
 172 }
 173
 174 #undef mpfr_set_divby0
 175
 176 void
 177 mpfr_set_divby0 (void)
 178 {
 179   __gmpfr_flags |= MPFR_FLAGS_DIVBY0;
 180 }
 181
 182 #undef mpfr_set_nanflag
 183
 184 void
 185 mpfr_set_nanflag (void)
 186 {
 187   __gmpfr_flags |= MPFR_FLAGS_NAN;
 188 }
 189
 190 #undef mpfr_set_inexflag
 191
 192 void
 193 mpfr_set_inexflag (void)
 194 {
 195   __gmpfr_flags |= MPFR_FLAGS_INEXACT;
 196 }
 197
 198 #undef mpfr_set_erangeflag
 199
 200 void
 201 mpfr_set_erangeflag (void)
 202 {
 203   __gmpfr_flags |= MPFR_FLAGS_ERANGE;
 204 }
 205
 206
 207 #undef mpfr_check_range
 208
 209 int
 210 mpfr_check_range (mpfr_ptr x, int t, mpfr_rnd_t rnd_mode)
 211 {
 212   if (MPFR_LIKELY( MPFR_IS_PURE_FP(x)) )
 213     { /* x is a non-zero FP */
 214       mpfr_exp_t exp = MPFR_EXP (x);  /* Do not use MPFR_GET_EXP */
 215       if (MPFR_UNLIKELY( exp < __gmpfr_emin) )
 216         {
 217           /* The following test is necessary because in the rounding to the
 218            * nearest mode, mpfr_underflow always rounds away from 0. In
 219            * this rounding mode, we need to round to 0 if:
 220            *   _ |x| < 2^(emin-2), or
 221            *   _ |x| = 2^(emin-2) and the absolute value of the exact
 222            *     result is <= 2^(emin-2).
 223            */
 224           if (rnd_mode == MPFR_RNDN &&
 225               (exp + 1 < __gmpfr_emin ||
 226                (mpfr_powerof2_raw(x) &&
 227                 (MPFR_IS_NEG(x) ? t <= 0 : t >= 0))))
 228             rnd_mode = MPFR_RNDZ;
 229           return mpfr_underflow(x, rnd_mode, MPFR_SIGN(x));
 230         }
 231       if (MPFR_UNLIKELY( exp > __gmpfr_emax) )
 232         return mpfr_overflow (x, rnd_mode, MPFR_SIGN(x));
 233     }
 234   else if (MPFR_UNLIKELY (t != 0 && MPFR_IS_INF (x)))
 235     {
 236       /* We need to do the following because most MPFR functions are
 237        * implemented in the following way:
 238        *   Ziv's loop:
 239        *   | Compute an approximation to the result and an error bound.
 240        *   | Possible underflow/overflow detection -> return.
 241        *   | If can_round, break (exit the loop).
 242        *   | Otherwise, increase the working precision and loop.
 243        *   Round the approximation in the target precision.  <== See below
 244        *   Restore the flags (that could have been set due to underflows
 245        *   or overflows during the internal computations).
 246        *   Execute: return mpfr_check_range (...).
 247        * The problem is that an overflow could be generated when rounding the
 248        * approximation (in general, such an overflow could not be detected
 249        * earlier), and the overflow flag is lost when the flags are restored.
 250        * This can occur only when the rounding yields an exponent change
 251        * and the new exponent is larger than the maximum exponent, so that
 252        * an infinity is necessarily obtained.
 253        * So, the simplest solution is to detect this overflow case here in
 254        * mpfr_check_range, which is easy to do since the rounded result is
 255        * necessarily an inexact infinity.
 256        */
 257       __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
 258     }
 259   MPFR_RET (t);  /* propagate inexact ternary value, unlike most functions */
 260 }
 261
 262 #undef mpfr_underflow_p
 263
 264 int
 265 mpfr_underflow_p (void)
 266 {
 267   return __gmpfr_flags & MPFR_FLAGS_UNDERFLOW;
 268 }
 269
 270 #undef mpfr_overflow_p
 271
 272 int
 273 mpfr_overflow_p (void)
 274 {
 275   return __gmpfr_flags & MPFR_FLAGS_OVERFLOW;
 276 }
 277
 278 #undef mpfr_divby0_p
 279
 280 int
 281 mpfr_divby0_p (void)
 282 {
 283   return __gmpfr_flags & MPFR_FLAGS_DIVBY0;
 284 }
 285
 286 #undef mpfr_nanflag_p
 287
 288 int
 289 mpfr_nanflag_p (void)
 290 {
 291   return __gmpfr_flags & MPFR_FLAGS_NAN;
 292 }
 293
 294 #undef mpfr_inexflag_p
 295
 296 int
 297 mpfr_inexflag_p (void)
 298 {
 299   return __gmpfr_flags & MPFR_FLAGS_INEXACT;
 300 }
 301
 302 #undef mpfr_erangeflag_p
 303
 304 int
 305 mpfr_erangeflag_p (void)
 306 {
 307   return __gmpfr_flags & MPFR_FLAGS_ERANGE;
 308 }
 309
 310 /* #undef mpfr_underflow */
 311
 312 /* Note: In the rounding to the nearest mode, mpfr_underflow
 313    always rounds away from 0. In this rounding mode, you must call
 314    mpfr_underflow with rnd_mode = MPFR_RNDZ if the exact result
 315    is <= 2^(emin-2) in absolute value. */
 316
 317 int
 318 mpfr_underflow (mpfr_ptr x, mpfr_rnd_t rnd_mode, int sign)
 319 {
 320   int inex;
 321
 322   MPFR_ASSERT_SIGN (sign);
 323
 324   if (MPFR_IS_LIKE_RNDZ(rnd_mode, sign < 0))
 325     {
 326       MPFR_SET_ZERO(x);
 327       inex = -1;
 328     }
 329   else
 330     {
 331       mpfr_setmin (x, __gmpfr_emin);
 332       inex = 1;
 333     }
 334   MPFR_SET_SIGN(x, sign);
 335   __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
 336   return sign > 0 ? inex : -inex;
 337 }
 338
 339 /* #undef mpfr_overflow */
 340
 341 int
 342 mpfr_overflow (mpfr_ptr x, mpfr_rnd_t rnd_mode, int sign)
 343 {
 344   int inex;
 345
 346   MPFR_ASSERT_SIGN(sign);
 347   if (MPFR_IS_LIKE_RNDZ(rnd_mode, sign < 0))
 348     {
 349       mpfr_setmax (x, __gmpfr_emax);
 350       inex = -1;
 351     }
 352   else
 353     {
 354       MPFR_SET_INF(x);
 355       inex = 1;
 356     }
 357   MPFR_SET_SIGN(x,sign);
 358   __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
 359   return sign > 0 ? inex : -inex;
 360 }