summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/exp_2.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/exp_2.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/exp_2.c177
1 files changed, 106 insertions, 71 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/exp_2.c b/Build/source/libs/mpfr/mpfr-src/src/exp_2.c
index e90d462e654..64a9fd286d0 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/exp_2.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/exp_2.c
@@ -21,7 +21,6 @@ along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
-/* #define DEBUG */
#define MPFR_NEED_LONGLONG_H /* for count_leading_zeros */
#include "mpfr-impl.h"
@@ -44,13 +43,12 @@ mpz_normalize (mpz_t rop, mpz_t z, mpfr_exp_t q)
MPFR_MPZ_SIZEINBASE2 (k, z);
MPFR_ASSERTD (k == (mpfr_uexp_t) k);
- if (q < 0 || (mpfr_uexp_t) k > (mpfr_uexp_t) q)
+ if (MPFR_LIKELY(q < 0 || (mpfr_uexp_t) k > (mpfr_uexp_t) q))
{
mpz_fdiv_q_2exp (rop, z, (unsigned long) ((mpfr_uexp_t) k - q));
return (mpfr_exp_t) k - q;
}
- if (MPFR_UNLIKELY(rop != z))
- mpz_set (rop, z);
+ mpz_set (rop, z);
return 0;
}
@@ -61,7 +59,7 @@ mpz_normalize (mpz_t rop, mpz_t z, mpfr_exp_t q)
static mpfr_exp_t
mpz_normalize2 (mpz_t rop, mpz_t z, mpfr_exp_t expz, mpfr_exp_t target)
{
- if (target > expz)
+ if (MPFR_LIKELY(target > expz))
mpz_fdiv_q_2exp (rop, z, target - expz);
else
mpz_mul_2exp (rop, z, expz - target);
@@ -83,8 +81,9 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_exp_t exps, expx;
mpfr_prec_t q, precy;
int inexact;
- mpfr_t r, s;
+ mpfr_t s, r;
mpz_t ss;
+ MPFR_GROUP_DECL(group);
MPFR_ZIV_DECL (loop);
MPFR_LOG_FUNC
@@ -95,20 +94,47 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
expx = MPFR_GET_EXP (x);
precy = MPFR_PREC(y);
- /* Warning: we cannot use the 'double' type here, since on 64-bit machines
+ /* First perform argument reduction: if x' = x - n*log(2) we have
+ exp(x) = exp(x)*2^n. We should take n near from x/log(2) but it does not
+ need to be exact.
+ Warning: we cannot use the 'double' type here, since on 64-bit machines
x may be as large as 2^62*log(2) without overflow, and then x/log(2)
is about 2^62: not every integer of that size can be represented as a
'double', thus the argument reduction would fail. */
- if (expx <= -2)
+ if (MPFR_UNLIKELY(expx <= -2))
/* |x| <= 0.25, thus n = round(x/log(2)) = 0 */
n = 0;
else
{
- mpfr_init2 (r, sizeof (long) * CHAR_BIT);
- mpfr_const_log2 (r, MPFR_RNDZ);
- mpfr_div (r, x, r, MPFR_RNDN);
+ mp_limb_t r_limb[(sizeof (long) -1) / sizeof(mp_limb_t) + 1];
+ /* Note: we use precision sizeof (long) * CHAR_BIT - 1 here since it is
+ more efficient that full limb precision. */
+ MPFR_TMP_INIT1(r_limb, r, sizeof (long) * CHAR_BIT - 1);
+ mpfr_div (r, x, __gmpfr_const_log2_RNDD, MPFR_RNDN);
+#ifdef MPFR_LONG_WITHIN_LIMB
+ /* The following code assume an unsigned long can fit in a mp_limb_t */
+ {
+ mp_limb_t a;
+ mpfr_exp_t exp;
+ MPFR_STAT_STATIC_ASSERT (MPFR_LIMB_MAX >= ULONG_MAX);
+ /* Read the long directly (faster than using mpfr_get_si
+ since it fits, it is not singular, it can't be zero
+ and there is no conversion to do) */
+ MPFR_ASSERTD (MPFR_NOTZERO (r));
+ exp = MPFR_GET_EXP (r);
+ MPFR_ASSERTD (exp <= GMP_NUMB_BITS);
+ if (exp >= 1)
+ {
+ a = MPFR_MANT(r)[0] >> (GMP_NUMB_BITS - exp);
+ n = MPFR_IS_POS (r) ? a : a <= LONG_MAX ? - (long) a : LONG_MIN;
+ }
+ else
+ n = 0;
+ }
+#else
+ /* Use generic way to get the long */
n = mpfr_get_si (r, MPFR_RNDN);
- mpfr_clear (r);
+#endif
}
/* we have |x| <= (|n|+1)*log(2) */
MPFR_LOG_MSG (("d(x)=%1.30e n=%ld\n", mpfr_get_d1(x), n));
@@ -118,7 +144,8 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
error_r = 0;
else
{
- count_leading_zeros (error_r, (mp_limb_t) SAFE_ABS (unsigned long, n) + 1);
+ count_leading_zeros (error_r,
+ (mp_limb_t) SAFE_ABS (unsigned long, n) + 1);
error_r = GMP_NUMB_BITS - error_r;
/* we have |x| <= 2^error_r * log(2) */
}
@@ -126,21 +153,21 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
/* for the O(n^(1/2)*M(n)) method, the Taylor series computation of
n/K terms costs about n/(2K) multiplications when computed in fixed
point */
- K = (precy < MPFR_EXP_2_THRESHOLD) ? __gmpfr_isqrt ((precy + 1) / 2)
+ K = (precy < MPFR_EXP_2_THRESHOLD) ? __gmpfr_isqrt ((precy + 1) / 2) + 3
: __gmpfr_cuberoot (4*precy);
l = (precy - 1) / K + 1;
err = K + MPFR_INT_CEIL_LOG2 (2 * l + 18);
/* add K extra bits, i.e. failure probability <= 1/2^K = O(1/precy) */
- q = precy + err + K + 8;
+ q = precy + err + K + 10;
/* if |x| >> 1, take into account the cancelled bits */
if (expx > 0)
q += expx;
- /* Note: due to the mpfr_prec_round below, it is not possible to use
- the MPFR_GROUP_* macros here. */
-
- mpfr_init2 (r, q + error_r);
- mpfr_init2 (s, q + error_r);
+ /* Even with to the mpfr_prec_round below, it is possible to use
+ the MPFR_GROUP_* macros here because mpfr_prec_round is only
+ called in a special case. */
+ MPFR_GROUP_INIT_2(group, q + error_r, r, s);
+ mpz_init (ss);
/* the algorithm consists in computing an upper bound of exp(x) using
a precision of q bits, and see if we can round to MPFR_PREC(y) taking
@@ -171,27 +198,31 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_sub (r, x, r, MPFR_RNDU);
- if (MPFR_IS_PURE_FP (r))
- {
- while (MPFR_IS_NEG (r))
- { /* initial approximation n was too large */
- n--;
- mpfr_add (r, r, s, MPFR_RNDU);
- }
+ while (MPFR_IS_PURE_FP(r) && MPFR_IS_NEG (r))
+ { /* initial approximation n was too large */
+ n--;
+ mpfr_add (r, r, s, MPFR_RNDU);
+ }
+ /* if r is 0, we cannot round correctly */
+ if (MPFR_LIKELY(MPFR_IS_PURE_FP (r)))
+ {
/* since there was a cancellation in x - n*log(2), the low error_r
bits from r are zero and thus non significant, thus we can reduce
the working precision */
- if (error_r > 0)
- mpfr_prec_round (r, q, MPFR_RNDU);
+ if (MPFR_LIKELY(error_r > 0))
+ {
+ MPFR_ASSERTD( MPFR_PREC(r) > q);
+ /* since MPFR_PREC(r) > q, there is no reallocation to do,
+ and it is safe to use it with MPFR_GROUP functions */
+ mpfr_prec_round (r, q, MPFR_RNDU);
+ }
/* the error on r is at most 3 ulps (3 ulps if error_r = 0,
and 1 + 3/2 if error_r > 0) */
MPFR_LOG_VAR (r);
MPFR_ASSERTD (MPFR_IS_POS (r));
mpfr_div_2ui (r, r, K, MPFR_RNDU); /* r = (x-n*log(2))/2^K, exact */
- mpz_init (ss);
- exps = mpfr_get_z_2exp (ss, s);
/* s <- 1 + r/1! + r^2/2! + ... + r^l/l! */
MPFR_ASSERTD (MPFR_IS_PURE_FP (r) && MPFR_EXP (r) < 0);
l = (precy < MPFR_EXP_2_THRESHOLD)
@@ -207,10 +238,7 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
exps *= 2;
exps += mpz_normalize (ss, ss, q);
}
- mpfr_set_z (s, ss, MPFR_RNDN);
-
- MPFR_SET_EXP(s, MPFR_GET_EXP (s) + exps);
- mpz_clear (ss);
+ mpfr_set_z_2exp (s, ss, exps, MPFR_RNDN);
/* error is at most 2^K*l, plus 2 to take into account of
the error of 3 ulps on r */
@@ -222,20 +250,17 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
if (MPFR_LIKELY (MPFR_CAN_ROUND (s, q - err, precy, rnd_mode)))
{
- mpfr_clear_flags ();
+ MPFR_CLEAR_FLAGS ();
inexact = mpfr_mul_2si (y, s, n, rnd_mode);
break;
}
}
-
MPFR_ZIV_NEXT (loop, q);
- mpfr_set_prec (r, q + error_r);
- mpfr_set_prec (s, q + error_r);
+ MPFR_GROUP_REPREC_2(group, q+error_r, r, s);
}
MPFR_ZIV_FREE (loop);
-
- mpfr_clear (r);
- mpfr_clear (s);
+ mpz_clear (ss);
+ MPFR_GROUP_CLEAR (group);
return inexact;
}
@@ -246,10 +271,6 @@ mpfr_exp_2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
The absolute error on s is less than 3*l*(l+1)*2^(-q).
Version using fixed-point arithmetic with mpz instead
of mpfr for internal computations.
- NOTE[VL]: the following sentence seems to be obsolete since MY_INIT_MPZ
- is no longer used (r6919); qn was the number of limbs of q.
- s must have at least qn+1 limbs (qn should be enough, but currently fails
- since mpz_mul_2exp(s, s, q-1) reallocates qn+1 limbs)
*/
static unsigned long
mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, mpfr_prec_t q, mpfr_exp_t *exps)
@@ -259,38 +280,52 @@ mpfr_exp2_aux (mpz_t s, mpfr_srcptr r, mpfr_prec_t q, mpfr_exp_t *exps)
mpz_t t, rr;
mp_size_t sbit, tbit;
- MPFR_ASSERTN (MPFR_IS_PURE_FP (r));
+ MPFR_ASSERTD (MPFR_IS_PURE_FP (r));
expt = 0;
*exps = 1 - (mpfr_exp_t) q; /* s = 2^(q-1) */
mpz_init (t);
mpz_init (rr);
- mpz_set_ui(t, 1);
- mpz_set_ui(s, 1);
- mpz_mul_2exp(s, s, q-1);
- expr = mpfr_get_z_2exp(rr, r); /* no error here */
+ mpz_set_ui (t, 1);
+ mpz_set_ui (s, 1);
+ mpz_mul_2exp (s, s, q-1);
+ expr = mpfr_get_z_2exp (rr, r); /* no error here */
l = 0;
- for (;;) {
- l++;
- mpz_mul(t, t, rr);
- expt += expr;
- MPFR_MPZ_SIZEINBASE2 (sbit, s);
- MPFR_MPZ_SIZEINBASE2 (tbit, t);
- dif = *exps + sbit - expt - tbit;
- /* truncates the bits of t which are < ulp(s) = 2^(1-q) */
- expt += mpz_normalize(t, t, (mpfr_exp_t) q-dif); /* error at most 2^(1-q) */
- mpz_fdiv_q_ui (t, t, l); /* error at most 2^(1-q) */
- /* the error wrt t^l/l! is here at most 3*l*ulp(s) */
- MPFR_ASSERTD (expt == *exps);
- if (mpz_sgn (t) == 0)
- break;
- mpz_add(s, s, t); /* no error here: exact */
- /* ensures rr has the same size as t: after several shifts, the error
- on rr is still at most ulp(t)=ulp(s) */
- MPFR_MPZ_SIZEINBASE2 (tbit, t);
- expr += mpz_normalize(rr, rr, tbit);
- }
+ for (;;)
+ {
+ l++;
+ mpz_mul(t, t, rr);
+ expt += expr;
+ MPFR_MPZ_SIZEINBASE2 (sbit, s);
+ MPFR_MPZ_SIZEINBASE2 (tbit, t);
+ dif = *exps + sbit - expt - tbit;
+ /* truncates the bits of t which are < ulp(s) = 2^(1-q) */
+ expt += mpz_normalize (t, t, (mpfr_exp_t) q - dif);
+ /* error at most 2^(1-q) */
+ if (l > 1)
+ {
+ /* GMP doesn't optimize the case of power of 2 */
+ if (IS_POW2(l))
+ {
+ int bits = MPFR_INT_CEIL_LOG2(l);
+ mpz_fdiv_q_2exp (t, t, bits); /* error at most 2^(1-q) */
+ }
+ else
+ {
+ mpz_fdiv_q_ui (t, t, l); /* error at most 2^(1-q) */
+ }
+ /* the error wrt t^l/l! is here at most 3*l*ulp(s) */
+ MPFR_ASSERTD (expt == *exps);
+ }
+ if (mpz_sgn (t) == 0)
+ break;
+ mpz_add (s, s, t); /* no error here: exact */
+ /* ensures rr has the same size as t: after several shifts, the error
+ on rr is still at most ulp(t)=ulp(s) */
+ MPFR_MPZ_SIZEINBASE2 (tbit, t);
+ expr += mpz_normalize (rr, rr, tbit);
+ }
mpz_clear (t);
mpz_clear (rr);