summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c96
1 files changed, 63 insertions, 33 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c b/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c
index ce397785234..983dc77a4a6 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/mpn_exp.c
@@ -1,6 +1,6 @@
/* mpfr_mpn_exp -- auxiliary function for mpfr_get_str and mpfr_set_str
-Copyright 1999-2019 Free Software Foundation, Inc.
+Copyright 1999-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
@@ -24,7 +24,43 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#define MPFR_NEED_LONGLONG_H
#include "mpfr-impl.h"
-/* this function computes an approximation of b^e in {a, n}, with exponent
+/* count the number of significant bits of e, i.e.,
+ nbits(mpfr_exp_t) - count_leading_zeros (e) */
+static int
+nbits_mpfr_exp_t (mpfr_exp_t e)
+{
+ int nbits = 0;
+
+ MPFR_ASSERTD (e > 0);
+ while (e >= 0x10000)
+ {
+ e >>= 16;
+ nbits += 16;
+ }
+ MPFR_ASSERTD (e <= 0xffff);
+ if (e >= 0x100)
+ {
+ e >>= 8;
+ nbits += 8;
+ }
+ MPFR_ASSERTD (e <= 0xff);
+ if (e >= 0x10)
+ {
+ e >>= 4;
+ nbits += 4;
+ }
+ MPFR_ASSERTD (e <= 0xf);
+ if (e >= 4)
+ {
+ e >>= 2;
+ nbits += 2;
+ }
+ MPFR_ASSERTD (e <= 3);
+ /* now e = 1, 2, or 3 */
+ return nbits + 1 + (e >= 2);
+}
+
+/* this function computes an approximation to b^e in {a, n}, with exponent
stored in exp_r. The computed value is rounded toward zero (truncated).
It returns an integer f such that the final error is bounded by 2^f ulps,
that is:
@@ -36,23 +72,23 @@ https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
Return -2 if an overflow occurred in the computation of exp_r.
*/
-long
+int
mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
{
mp_limb_t *c, B;
mpfr_exp_t f, h;
int i;
- unsigned long t; /* number of bits in e */
- unsigned long bits;
- size_t n1;
- unsigned int error; /* (number - 1) of loop a^2b inexact */
+ int t; /* number of bits in e */
+ size_t bits, n1;
+ unsigned int error; /* (number - 1) of loops a^2b inexact */
/* error == t means no error */
int err_s_a2 = 0;
int err_s_ab = 0; /* number of error when shift A^2, AB */
MPFR_TMP_DECL(marker);
- MPFR_ASSERTN(e > 0);
- MPFR_ASSERTN((2 <= b) && (b <= 62));
+ MPFR_ASSERTN (n > 0 && n <= ((size_t) -1) / GMP_NUMB_BITS);
+ MPFR_ASSERTN (e > 0);
+ MPFR_ASSERTN (2 <= b && b <= 62);
MPFR_TMP_MARK(marker);
@@ -75,17 +111,16 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
f = h - (n - 1) * GMP_NUMB_BITS;
/* determine number of bits in e */
- count_leading_zeros (t, (mp_limb_t) e);
+ t = nbits_mpfr_exp_t (e);
- t = GMP_NUMB_BITS - t; /* number of bits of exponent e */
-
- error = t; /* error <= GMP_NUMB_BITS */
+ error = t;
+ /* t, error <= bitsize(mpfr_exp_t) */
+ MPFR_ASSERTD (error >= 0);
MPN_ZERO (c, 2 * n);
for (i = t - 2; i >= 0; i--)
{
-
/* determine precision needed */
bits = n * GMP_NUMB_BITS - mpn_scan1 (a, 0);
n1 = (n * GMP_NUMB_BITS - bits) / GMP_NUMB_BITS;
@@ -98,19 +133,19 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
/* set {c+n, 2n1-n} to 0 : {c, n} = {a, n}^2*K^n */
/* check overflow on f */
- if (MPFR_UNLIKELY(f < MPFR_EXP_MIN/2 || f > MPFR_EXP_MAX/2))
+ if (MPFR_UNLIKELY (f < MPFR_EXP_MIN / 2 || f > MPFR_EXP_MAX / 2))
{
overflow:
MPFR_TMP_FREE(marker);
return -2;
}
- /* FIXME: Could f = 2*f + n * GMP_NUMB_BITS be used? */
- f = 2*f;
+ /* FIXME: Could f = 2 * f + n * GMP_NUMB_BITS be used? */
+ f = 2 * f;
MPFR_SADD_OVERFLOW (f, f, n * GMP_NUMB_BITS,
mpfr_exp_t, mpfr_uexp_t,
MPFR_EXP_MIN, MPFR_EXP_MAX,
goto overflow, goto overflow);
- if ((c[2*n - 1] & MPFR_LIMB_HIGHBIT) == 0)
+ if (MPFR_LIMB_MSB (c[2*n - 1]) == 0)
{
/* shift A by one bit to the left */
mpn_lshift (a, c + n, n, 1);
@@ -122,11 +157,11 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
else
MPN_COPY (a, c + n, n);
- if ((error == t) && (2 * n1 <= n) &&
- (mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS))
+ if (error == t && 2 * n1 <= n &&
+ mpn_scan1 (c + 2 * n1, 0) < (n - 2 * n1) * GMP_NUMB_BITS)
error = i;
- if (e & ((mpfr_exp_t) 1 << i))
+ if ((e >> i) & 1)
{
/* multiply A by B */
c[2 * n - 1] = mpn_mul_1 (c + n - 1, a, n, B);
@@ -143,18 +178,22 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
if (error != t)
err_s_ab ++;
}
- if ((error == t) && (c[n - 1] != 0))
+ if (error == t && c[n - 1] != 0)
error = i;
}
}
+ MPFR_ASSERTD (error >= 0);
+ MPFR_ASSERTD (err_s_a2 >= 0);
+ MPFR_ASSERTD (err_s_ab >= 0);
+
MPFR_TMP_FREE(marker);
*exp_r = f;
if (error == t)
return -1; /* result is exact */
- else /* error <= t-2 <= GMP_NUMB_BITS-2
+ else /* error <= t-2 <= bitsize(mpfr_exp_t)-2
err_s_ab, err_s_a2 <= t-1 */
{
/* if there are p loops after the first inexact result, with
@@ -162,15 +201,6 @@ mpfr_mpn_exp (mp_limb_t *a, mpfr_exp_t *exp_r, int b, mpfr_exp_t e, size_t n)
at most 2^(p+ceil((j+1)/2)+l+1)*ulp(res).
This is bounded by 2^(5/2*t-1/2) where t is the number of bits of e.
*/
- error = error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */
-#if 0
- if ((error - 1) >= ((n * GMP_NUMB_BITS - 1) / 2))
- error = n * GMP_NUMB_BITS; /* result is completely wrong:
- this is very unlikely since error is
- at most 5/2*log_2(e), and
- n * GMP_NUMB_BITS is at least
- 3*log_2(e) */
-#endif
- return error;
+ return error + err_s_ab + err_s_a2 / 2 + 3; /* <= 5t/2-1/2 */
}
}