contrib/mpfr/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 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 /* #define DEBUG */
  26
  27 int
  28 mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
  29 {
  30   mp_exp_t expx;
  31   mp_prec_t precy;
  32   int inexact;
  33   MPFR_SAVE_EXPO_DECL (expo);
  34
  35   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode),
  36                  ("y[%#R]=%R inexact=%d", y, y, inexact));
  37
  38   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
  39     {
  40       if (MPFR_IS_NAN(x))
  41         {
  42           MPFR_SET_NAN(y);
  43           MPFR_RET_NAN;
  44         }
  45       else if (MPFR_IS_INF(x))
  46         {
  47           if (MPFR_IS_POS(x))
  48             MPFR_SET_INF(y);
  49           else
  50             MPFR_SET_ZERO(y);
  51           MPFR_SET_POS(y);
  52           MPFR_RET(0);
  53         }
  54       else
  55         {
  56           MPFR_ASSERTD(MPFR_IS_ZERO(x));
  57           return mpfr_set_ui (y, 1, rnd_mode);
  58         }
  59     }
  60
  61   /* First, let's detect most overflow and underflow cases. */
  62   {
  63     mpfr_t e, bound;
  64
  65     /* We must extended the exponent range and save the flags now. */
  66     MPFR_SAVE_EXPO_MARK (expo);
  67
  68     mpfr_init2 (e, sizeof (mp_exp_t) * CHAR_BIT);
  69     mpfr_init2 (bound, 32);
  70
  71     inexact = mpfr_set_exp_t (e, expo.saved_emax, GMP_RNDN);
  72     MPFR_ASSERTD (inexact == 0);
  73     mpfr_const_log2 (bound, expo.saved_emax < 0 ? GMP_RNDD : GMP_RNDU);
  74     mpfr_mul (bound, bound, e, GMP_RNDU);
  75     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) >= 0))
  76       {
  77         /* x > log(2^emax), thus exp(x) > 2^emax */
  78         mpfr_clears (e, bound, (mpfr_ptr) 0);
  79         MPFR_SAVE_EXPO_FREE (expo);
  80         return mpfr_overflow (y, rnd_mode, 1);
  81       }
  82
  83     inexact = mpfr_set_exp_t (e, expo.saved_emin, GMP_RNDN);
  84     MPFR_ASSERTD (inexact == 0);
  85     inexact = mpfr_sub_ui (e, e, 2, GMP_RNDN);
  86     MPFR_ASSERTD (inexact == 0);
  87     mpfr_const_log2 (bound, expo.saved_emin < 0 ? GMP_RNDU : GMP_RNDD);
  88     mpfr_mul (bound, bound, e, GMP_RNDD);
  89     if (MPFR_UNLIKELY (mpfr_cmp (x, bound) <= 0))
  90       {
  91         /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
  92         mpfr_clears (e, bound, (mpfr_ptr) 0);
  93         MPFR_SAVE_EXPO_FREE (expo);
  94         return mpfr_underflow (y, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
  95                                1);
  96       }
  97
  98     /* Other overflow/underflow cases must be detected
  99        by the generic routines. */
 100     mpfr_clears (e, bound, (mpfr_ptr) 0);
 101     MPFR_SAVE_EXPO_FREE (expo);
 102   }
 103
 104   MPFR_CLEAR_FLAGS (y);
 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       mp_exp_t emin = __gmpfr_emin;
 113       mp_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 == GMP_RNDD ||
 118                                        rnd_mode == GMP_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 == GMP_RNDU)
 131             {
 132               mp_size_t yn;
 133               int sh;
 134
 135               yn = 1 + (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB;
 136               sh = (mp_prec_t) yn * BITS_PER_MP_LIMB - MPFR_PREC(y);
 137               MPFR_MANT(y)[0] += MPFR_LIMB_ONE << sh;
 138               inexact = 1;
 139             }
 140           else
 141             inexact = -MPFR_FROM_SIGN_TO_INT(signx);
 142         }
 143
 144       __gmpfr_emin = emin;
 145       __gmpfr_emax = emax;
 146     }
 147   else  /* General case */
 148     {
 149       if (MPFR_UNLIKELY (precy > MPFR_EXP_THRESHOLD))
 150         /* mpfr_exp_3 saves the exponent range and flags itself, otherwise
 151            the flag changes in mpfr_exp_3 are lost */
 152         inexact = mpfr_exp_3 (y, x, rnd_mode); /* O(M(n) log(n)^2) */
 153       else
 154         {
 155           MPFR_SAVE_EXPO_MARK (expo);
 156           inexact = mpfr_exp_2 (y, x, rnd_mode); /* O(n^(1/3) M(n)) */
 157           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 158           MPFR_SAVE_EXPO_FREE (expo);
 159         }
 160     }
 161
 162   return mpfr_check_range (y, inexact, rnd_mode);
 163 }