contrib/mpfr/src/exp.c

   1 /* mpfr_exp -- exponential of a floating-point number
   2
   3 Copyright 1999, 2000, 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 /* #define DEBUG */
  26
  27 int
  28 mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
  29 {
  30   mpfr_exp_t expx;
  31   mpfr_prec_t precy;
  32   int inexact;
  33   MPFR_SAVE_EXPO_DECL (expo);
  34
  35   MPFR_LOG_FUNC
  36     (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
  37      ("y[%Pu]=%.*Rg inexact=%d",
  38       mpfr_get_prec (y), mpfr_log_prec, y, inexact));
  39
  40   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
  41     {
  42       if (MPFR_IS_NAN(x))
  43         {
  44           MPFR_SET_NAN(y);
  45           MPFR_RET_NAN;
  46         }
  47       else if (MPFR_IS_INF(x))
  48         {
  49           if (MPFR_IS_POS(x))
  50             MPFR_SET_INF(y);
  51           else
  52             MPFR_SET_ZERO(y);
  53           MPFR_SET_POS(y);
  54           MPFR_RET(0);
  55         }
  56       else
  57         {
  58           MPFR_ASSERTD(MPFR_IS_ZERO(x));
  59           return mpfr_set_ui (y, 1, rnd_mode);
  60         }
  61     }
  62
  63   /* First, let's detect most overflow and underflow cases. */
  64   {
  65     mpfr_t e, bound;
  66
  67     /* We must extended the exponent range and save the flags now. */
  68     MPFR_SAVE_EXPO_MARK (expo);
  69
  70     mpfr_init2 (e, sizeof (mpfr_exp_t) * CHAR_BIT);
  71     mpfr_init2 (bound, 32);
  72
  73     inexact = mpfr_set_exp_t (e, expo.saved_emax, MPFR_RNDN);
  74     MPFR_ASSERTD (inexact == 0);
  75     mpfr_const_log2 (bound, expo.saved_emax < 0 ? MPFR_RNDD : MPFR_RNDU);
  76     mpfr_mul (bound, bound, e, MPFR_RNDU);
  77     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) >= 0))
  78       {
  79         /* x > log(2^emax), thus exp(x) > 2^emax */
  80         mpfr_clears (e, bound, (mpfr_ptr) 0);
  81         MPFR_SAVE_EXPO_FREE (expo);
  82         return mpfr_overflow (y, rnd_mode, 1);
  83       }
  84
  85     inexact = mpfr_set_exp_t (e, expo.saved_emin, MPFR_RNDN);
  86     MPFR_ASSERTD (inexact == 0);
  87     inexact = mpfr_sub_ui (e, e, 2, MPFR_RNDN);
  88     MPFR_ASSERTD (inexact == 0);
  89     mpfr_const_log2 (bound, expo.saved_emin < 0 ? MPFR_RNDU : MPFR_RNDD);
  90     mpfr_mul (bound, bound, e, MPFR_RNDD);
  91     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) <= 0))
  92       {
  93         /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
  94         mpfr_clears (e, bound, (mpfr_ptr) 0);
  95         MPFR_SAVE_EXPO_FREE (expo);
  96         return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
  97                                1);
  98       }
  99
 100     /* Other overflow/underflow cases must be detected
 101        by the generic routines. */
 102     mpfr_clears (e, bound, (mpfr_ptr) 0);
 103     MPFR_SAVE_EXPO_FREE (expo);
 104   }
 105
 106   expx  = MPFR_GET_EXP (x);
 107   precy = MPFR_PREC (y);
 108
 109   /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
 110   if (MPFR_UNLIKELY (expx < 0 && (mpfr_uexp_t) (-expx) > precy))
 111     {
 112       mpfr_exp_t emin = __gmpfr_emin;
 113       mpfr_exp_t emax = __gmpfr_emax;
 114       int signx = MPFR_SIGN (x);
 115
 116       MPFR_SET_POS (y);
 117       if (MPFR_IS_NEG_SIGN (signx) && (rnd_mode == MPFR_RNDD ||
 118                                        rnd_mode == MPFR_RNDZ))
 119         {
 120           __gmpfr_emin = 0;
 121           __gmpfr_emax = 0;
 122           mpfr_setmax (y, 0);  /* y = 1 - epsilon */
 123           inexact = -1;
 124         }
 125       else
 126         {
 127           __gmpfr_emin = 1;
 128           __gmpfr_emax = 1;
 129           mpfr_setmin (y, 1);  /* y = 1 */
 130           if (MPFR_IS_POS_SIGN (signx) && (rnd_mode == MPFR_RNDU ||
 131                                            rnd_mode == MPFR_RNDA))
 132             {
 133               mp_size_t yn;
 134               int sh;
 135
 136               yn = 1 + (MPFR_PREC(y) - 1) / GMP_NUMB_BITS;
 137               sh = (mpfr_prec_t) yn * GMP_NUMB_BITS - MPFR_PREC(y);
 138               MPFR_MANT(y)[0] += MPFR_LIMB_ONE << sh;
 139               inexact = 1;
 140             }
 141           else
 142             inexact = -MPFR_FROM_SIGN_TO_INT(signx);
 143         }
 144
 145       __gmpfr_emin = emin;
 146       __gmpfr_emax = emax;
 147     }
 148   else  /* General case */
 149     {
 150       if (MPFR_UNLIKELY (precy >= MPFR_EXP_THRESHOLD))
 151         /* mpfr_exp_3 saves the exponent range and flags itself, otherwise
 152            the flag changes in mpfr_exp_3 are lost */
 153         inexact = mpfr_exp_3 (y, x, rnd_mode); /* O(M(n) log(n)^2) */
 154       else
 155         {
 156           MPFR_SAVE_EXPO_MARK (expo);
 157           inexact = mpfr_exp_2 (y, x, rnd_mode); /* O(n^(1/3) M(n)) */
 158           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 159           MPFR_SAVE_EXPO_FREE (expo);
 160         }
 161     }
 162
 163   return mpfr_check_range (y, inexact, rnd_mode);
 164 }