summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/exp.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/exp.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/exp.c129
1 files changed, 76 insertions, 53 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/exp.c b/Build/source/libs/mpfr/mpfr-src/src/exp.c
index e880cbd990c..aac0a4ffa61 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/exp.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/exp.c
@@ -22,7 +22,15 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#include "mpfr-impl.h"
-/* #define DEBUG */
+/* Cache for emin and emax bounds.
+ Contrary to other caches, it uses a fixed size for the mantissa,
+ so there is no dynamic allocation, and no need to free them. */
+static MPFR_THREAD_ATTR mpfr_exp_t previous_emin;
+static MPFR_THREAD_ATTR mp_limb_t bound_emin_limb[(32 - 1) / GMP_NUMB_BITS + 1];
+static MPFR_THREAD_ATTR mpfr_t bound_emin;
+static MPFR_THREAD_ATTR mpfr_exp_t previous_emax;
+static MPFR_THREAD_ATTR mp_limb_t bound_emax_limb[(32 - 1) / GMP_NUMB_BITS + 1];
+static MPFR_THREAD_ATTR mpfr_t bound_emax;
int
mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
@@ -61,81 +69,96 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
}
/* First, let's detect most overflow and underflow cases. */
- {
- mpfr_t e, bound;
-
- /* We must extended the exponent range and save the flags now. */
- MPFR_SAVE_EXPO_MARK (expo);
-
- mpfr_init2 (e, sizeof (mpfr_exp_t) * CHAR_BIT);
- mpfr_init2 (bound, 32);
-
- inexact = mpfr_set_exp_t (e, expo.saved_emax, MPFR_RNDN);
- MPFR_ASSERTD (inexact == 0);
- mpfr_const_log2 (bound, expo.saved_emax < 0 ? MPFR_RNDD : MPFR_RNDU);
- mpfr_mul (bound, bound, e, MPFR_RNDU);
- if (MPFR_UNLIKELY (mpfr_cmp (x, bound) >= 0))
- {
- /* x > log(2^emax), thus exp(x) > 2^emax */
- mpfr_clears (e, bound, (mpfr_ptr) 0);
- MPFR_SAVE_EXPO_FREE (expo);
- return mpfr_overflow (y, rnd_mode, 1);
- }
-
- inexact = mpfr_set_exp_t (e, expo.saved_emin, MPFR_RNDN);
- MPFR_ASSERTD (inexact == 0);
- inexact = mpfr_sub_ui (e, e, 2, MPFR_RNDN);
- MPFR_ASSERTD (inexact == 0);
- mpfr_const_log2 (bound, expo.saved_emin < 0 ? MPFR_RNDU : MPFR_RNDD);
- mpfr_mul (bound, bound, e, MPFR_RNDD);
- if (MPFR_UNLIKELY (mpfr_cmp (x, bound) <= 0))
- {
- /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
- mpfr_clears (e, bound, (mpfr_ptr) 0);
- MPFR_SAVE_EXPO_FREE (expo);
- return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
- 1);
- }
-
- /* Other overflow/underflow cases must be detected
- by the generic routines. */
- mpfr_clears (e, bound, (mpfr_ptr) 0);
- MPFR_SAVE_EXPO_FREE (expo);
- }
+ /* emax bound is cached. Check if the value in cache is ok */
+ if (MPFR_UNLIKELY (mpfr_get_emax() != previous_emax))
+ {
+ /* Recompute the emax bound */
+ mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
+ mpfr_t e;
+
+ /* We extend the exponent range and save the flags. */
+ MPFR_SAVE_EXPO_MARK (expo);
+
+ MPFR_TMP_INIT1(e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT);
+ MPFR_TMP_INIT1(bound_emax_limb, bound_emax, 32);
+
+ inexact = mpfr_set_exp_t (e, expo.saved_emax, MPFR_RNDN);
+ MPFR_ASSERTD (inexact == 0);
+ mpfr_mul (bound_emax, expo.saved_emax < 0 ?
+ __gmpfr_const_log2_RNDD : __gmpfr_const_log2_RNDU,
+ e, MPFR_RNDU);
+ previous_emax = expo.saved_emax;
+ MPFR_SAVE_EXPO_FREE (expo);
+ }
+
+ /* mpfr_cmp works even in reduced emin,emax range */
+ if (MPFR_UNLIKELY (mpfr_cmp (x, bound_emax) >= 0))
+ {
+ /* x > log(2^emax), thus exp(x) > 2^emax */
+ return mpfr_overflow (y, rnd_mode, 1);
+ }
+
+ /* emin bound is cached. Check if the value in cache is ok */
+ if (MPFR_UNLIKELY (mpfr_get_emin() != previous_emin))
+ {
+ mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
+ mpfr_t e;
+
+ /* We extend the exponent range and save the flags. */
+ MPFR_SAVE_EXPO_MARK (expo);
+
+ /* avoid using a full limb since arithmetic might be slower */
+ MPFR_TMP_INIT1(e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT - 1);
+ MPFR_TMP_INIT1(bound_emin_limb, bound_emin, 32);
+
+ inexact = mpfr_set_exp_t (e, expo.saved_emin, MPFR_RNDN);
+ MPFR_ASSERTD (inexact == 0);
+ inexact = mpfr_sub_ui (e, e, 2, MPFR_RNDN);
+ MPFR_ASSERTD (inexact == 0);
+ mpfr_const_log2 (bound_emin, expo.saved_emin < 0 ? MPFR_RNDU : MPFR_RNDD);
+ mpfr_mul (bound_emin, bound_emin, e, MPFR_RNDD);
+ previous_emin = expo.saved_emin;
+ MPFR_SAVE_EXPO_FREE (expo);
+ }
+
+ if (MPFR_UNLIKELY (mpfr_cmp (x, bound_emin) <= 0))
+ {
+ /* x < log(2^(emin - 2)), thus exp(x) < 2^(emin - 2) */
+ return mpfr_underflow (y, rnd_mode == MPFR_RNDN ? MPFR_RNDZ : rnd_mode,
+ 1);
+ }
expx = MPFR_GET_EXP (x);
precy = MPFR_PREC (y);
- /* if x < 2^(-precy), then exp(x) i.e. gives 1 +/- 1 ulp(1) */
+ /* if x < 2^(-precy), then exp(x) gives 1 +/- 1 ulp(1) */
if (MPFR_UNLIKELY (expx < 0 && (mpfr_uexp_t) (-expx) > precy))
{
mpfr_exp_t emin = __gmpfr_emin;
mpfr_exp_t emax = __gmpfr_emax;
int signx = MPFR_SIGN (x);
+ /* Make sure that the exponent range is large enough:
+ * [0,2] is sufficient in all precisions.
+ */
+ __gmpfr_emin = 0;
+ __gmpfr_emax = 2;
+
MPFR_SET_POS (y);
if (MPFR_IS_NEG_SIGN (signx) && (rnd_mode == MPFR_RNDD ||
rnd_mode == MPFR_RNDZ))
{
- __gmpfr_emin = 0;
- __gmpfr_emax = 0;
mpfr_setmax (y, 0); /* y = 1 - epsilon */
inexact = -1;
}
else
{
- __gmpfr_emin = 1;
- __gmpfr_emax = 1;
mpfr_setmin (y, 1); /* y = 1 */
if (MPFR_IS_POS_SIGN (signx) && (rnd_mode == MPFR_RNDU ||
rnd_mode == MPFR_RNDA))
{
- mp_size_t yn;
- int sh;
-
- yn = MPFR_LIMB_SIZE (y);
- sh = (mpfr_prec_t) yn * GMP_NUMB_BITS - MPFR_PREC(y);
- MPFR_MANT(y)[0] += MPFR_LIMB_ONE << sh;
+ /* Warning: should work for precision 1 bit too! */
+ mpfr_nextabove (y);
inexact = 1;
}
else