1 /* mpfr_exp -- exponential of a floating-point number
3 Copyright 1999-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #include "mpfr-impl.h"
28 mpfr_exp (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
33 MPFR_SAVE_EXPO_DECL (expo
);
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
));
40 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x
) ))
47 else if (MPFR_IS_INF(x
))
58 MPFR_ASSERTD(MPFR_IS_ZERO(x
));
59 return mpfr_set_ui (y
, 1, rnd_mode
);
63 /* First, let's detect most overflow and underflow cases. */
67 /* We must extended the exponent range and save the flags now. */
68 MPFR_SAVE_EXPO_MARK (expo
);
70 mpfr_init2 (e
, sizeof (mpfr_exp_t
) * CHAR_BIT
);
71 mpfr_init2 (bound
, 32);
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))
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);
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))
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
,
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
);
106 expx
= MPFR_GET_EXP (x
);
107 precy
= MPFR_PREC (y
);
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
))
112 mpfr_exp_t emin
= __gmpfr_emin
;
113 mpfr_exp_t emax
= __gmpfr_emax
;
114 int signx
= MPFR_SIGN (x
);
117 if (MPFR_IS_NEG_SIGN (signx
) && (rnd_mode
== MPFR_RNDD
||
118 rnd_mode
== MPFR_RNDZ
))
122 mpfr_setmax (y
, 0); /* y = 1 - epsilon */
129 mpfr_setmin (y
, 1); /* y = 1 */
130 if (MPFR_IS_POS_SIGN (signx
) && (rnd_mode
== MPFR_RNDU
||
131 rnd_mode
== MPFR_RNDA
))
136 yn
= MPFR_LIMB_SIZE (y
);
137 sh
= (mpfr_prec_t
) yn
* GMP_NUMB_BITS
- MPFR_PREC(y
);
138 MPFR_MANT(y
)[0] += MPFR_LIMB_ONE
<< sh
;
142 inexact
= -MPFR_FROM_SIGN_TO_INT(signx
);
148 else /* General case */
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) */
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
);
163 return mpfr_check_range (y
, inexact
, rnd_mode
);