contrib/mpfr/exceptions.c

   1 /* Exception flags and utilities.
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #include "mpfr-impl.h"
  24
  25 unsigned int MPFR_THREAD_ATTR __gmpfr_flags = 0;
  26
  27 mp_exp_t MPFR_THREAD_ATTR __gmpfr_emin = MPFR_EMIN_DEFAULT;
  28 mp_exp_t MPFR_THREAD_ATTR __gmpfr_emax = MPFR_EMAX_DEFAULT;
  29
  30 #undef mpfr_get_emin
  31
  32 mp_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 (mp_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 mp_exp_t
  55 mpfr_get_emin_min (void)
  56 {
  57   return MPFR_EMIN_MIN;
  58 }
  59
  60 mp_exp_t
  61 mpfr_get_emin_max (void)
  62 {
  63   return MPFR_EMIN_MAX;
  64 }
  65
  66 #undef mpfr_get_emax
  67
  68 mp_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 (mp_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 mp_exp_t
  91 mpfr_get_emax_min (void)
  92 {
  93   return MPFR_EMAX_MIN;
  94 }
  95 mp_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_nanflag
 127
 128 void
 129 mpfr_clear_nanflag (void)
 130 {
 131   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_NAN;
 132 }
 133
 134 #undef mpfr_clear_inexflag
 135
 136 void
 137 mpfr_clear_inexflag (void)
 138 {
 139   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_INEXACT;
 140 }
 141
 142 #undef mpfr_clear_erangeflag
 143
 144 void
 145 mpfr_clear_erangeflag (void)
 146 {
 147   __gmpfr_flags &= MPFR_FLAGS_ALL ^ MPFR_FLAGS_ERANGE;
 148 }
 149
 150 #undef mpfr_clear_underflow
 151
 152 void
 153 mpfr_set_underflow (void)
 154 {
 155   __gmpfr_flags |= MPFR_FLAGS_UNDERFLOW;
 156 }
 157
 158 #undef mpfr_clear_overflow
 159
 160 void
 161 mpfr_set_overflow (void)
 162 {
 163   __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
 164 }
 165
 166 #undef mpfr_clear_nanflag
 167
 168 void
 169 mpfr_set_nanflag (void)
 170 {
 171   __gmpfr_flags |= MPFR_FLAGS_NAN;
 172 }
 173
 174 #undef mpfr_clear_inexflag
 175
 176 void
 177 mpfr_set_inexflag (void)
 178 {
 179   __gmpfr_flags |= MPFR_FLAGS_INEXACT;
 180 }
 181
 182 #undef mpfr_clear_erangeflag
 183
 184 void
 185 mpfr_set_erangeflag (void)
 186 {
 187   __gmpfr_flags |= MPFR_FLAGS_ERANGE;
 188 }
 189
 190
 191 #undef mpfr_check_range
 192
 193 int
 194 mpfr_check_range (mpfr_ptr x, int t, mp_rnd_t rnd_mode)
 195 {
 196   if (MPFR_LIKELY( MPFR_IS_PURE_FP(x)) )
 197     { /* x is a non-zero FP */
 198       mp_exp_t exp = MPFR_EXP (x);  /* Do not use MPFR_GET_EXP */
 199       if (MPFR_UNLIKELY( exp < __gmpfr_emin) )
 200         {
 201           /* The following test is necessary because in the rounding to the
 202            * nearest mode, mpfr_underflow always rounds away from 0. In
 203            * this rounding mode, we need to round to 0 if:
 204            *   _ |x| < 2^(emin-2), or
 205            *   _ |x| = 2^(emin-2) and the absolute value of the exact
 206            *     result is <= 2^(emin-2).
 207            */
 208           if (rnd_mode == GMP_RNDN &&
 209               (exp + 1 < __gmpfr_emin ||
 210                (mpfr_powerof2_raw(x) &&
 211                 (MPFR_IS_NEG(x) ? t <= 0 : t >= 0))))
 212             rnd_mode = GMP_RNDZ;
 213           return mpfr_underflow(x, rnd_mode, MPFR_SIGN(x));
 214         }
 215       if (MPFR_UNLIKELY( exp > __gmpfr_emax) )
 216         return mpfr_overflow(x, rnd_mode, MPFR_SIGN(x));
 217     }
 218   else if (MPFR_UNLIKELY (t != 0 && MPFR_IS_INF (x)))
 219     {
 220       /* We need to do the following because most MPFR functions are
 221        * implemented in the following way:
 222        *   Ziv's loop:
 223        *   | Compute an approximation to the result and an error bound.
 224        *   | Possible underflow/overflow detection -> return.
 225        *   | If can_round, break (exit the loop).
 226        *   | Otherwise, increase the working precision and loop.
 227        *   Round the approximation in the target precision.
 228        *   Restore the flags (that could have been set due to underflows
 229        *   or overflows during the internal computations).
 230        *   Execute: return mpfr_check_range (...).
 231        * The problem is that an overflow could be generated when rounding the
 232        * approximation (in general, such an overflow could not be detected
 233        * earlier), and the overflow flag is lost when the flags are restored.
 234        * So, the simplest solution is to detect this overflow case here in
 235        * mpfr_check_range, which is easy to do since the rounded result is
 236        * necessarily an inexact infinity.
 237        */
 238       __gmpfr_flags |= MPFR_FLAGS_OVERFLOW;
 239     }
 240   MPFR_RET (t);  /* propagate inexact ternary value, unlike most functions */
 241 }
 242
 243 #undef mpfr_underflow_p
 244
 245 int
 246 mpfr_underflow_p (void)
 247 {
 248   return __gmpfr_flags & MPFR_FLAGS_UNDERFLOW;
 249 }
 250
 251 #undef mpfr_overflow_p
 252
 253 int
 254 mpfr_overflow_p (void)
 255 {
 256   return __gmpfr_flags & MPFR_FLAGS_OVERFLOW;
 257 }
 258
 259 #undef mpfr_nanflag_p
 260
 261 int
 262 mpfr_nanflag_p (void)
 263 {
 264   return __gmpfr_flags & MPFR_FLAGS_NAN;
 265 }
 266
 267 #undef mpfr_inexflag_p
 268
 269 int
 270 mpfr_inexflag_p (void)
 271 {
 272   return __gmpfr_flags & MPFR_FLAGS_INEXACT;
 273 }
 274
 275 #undef mpfr_erangeflag_p
 276
 277 int
 278 mpfr_erangeflag_p (void)
 279 {
 280   return __gmpfr_flags & MPFR_FLAGS_ERANGE;
 281 }
 282
 283 /* #undef mpfr_underflow */
 284
 285 /* Note: In the rounding to the nearest mode, mpfr_underflow
 286    always rounds away from 0. In this rounding mode, you must call
 287    mpfr_underflow with rnd_mode = GMP_RNDZ if the exact result
 288    is <= 2^(emin-2) in absolute value. */
 289
 290 int
 291 mpfr_underflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
 292 {
 293   int inex;
 294
 295   MPFR_ASSERT_SIGN (sign);
 296
 297   if (rnd_mode == GMP_RNDN
 298       || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_POS_SIGN (sign)))
 299     {
 300       mpfr_setmin (x, __gmpfr_emin);
 301       inex = 1;
 302     }
 303   else
 304     {
 305       MPFR_SET_ZERO(x);
 306       inex = -1;
 307     }
 308   MPFR_SET_SIGN(x, sign);
 309   __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_UNDERFLOW;
 310   return sign > 0 ? inex : -inex;
 311 }
 312
 313 /* #undef mpfr_overflow */
 314
 315 int
 316 mpfr_overflow (mpfr_ptr x, mp_rnd_t rnd_mode, int sign)
 317 {
 318   int inex;
 319
 320   MPFR_ASSERT_SIGN(sign);
 321   MPFR_CLEAR_FLAGS(x);
 322   if (rnd_mode == GMP_RNDN
 323       || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, sign > 0))
 324     {
 325       MPFR_SET_INF(x);
 326       inex = 1;
 327     }
 328   else
 329     {
 330       mpfr_setmax (x, __gmpfr_emax);
 331       inex = -1;
 332     }
 333   MPFR_SET_SIGN(x,sign);
 334   __gmpfr_flags |= MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
 335   return sign > 0 ? inex : -inex;
 336 }