summaryrefslogtreecommitdiff
path: root/Build/source/libs/mpfr/mpfr-src/src
diff options
context:
space:
mode:
authorAkira Kakuto <kakuto@fuk.kindai.ac.jp>2016-09-28 08:13:29 +0000
committerAkira Kakuto <kakuto@fuk.kindai.ac.jp>2016-09-28 08:13:29 +0000
commitf7b738cb900737646435336e5932c089a4e3c32f (patch)
tree4ab47d8c518aa6e86216b3451c7f415ed217658a /Build/source/libs/mpfr/mpfr-src/src
parent460968dd6d10a96729c24a3f37e1f54928c39603 (diff)
mpfr 3.1.5
git-svn-id: svn://tug.org/texlive/trunk@42168 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src')
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/Makefile.am2
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/add_ui.c1
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/get_ld.c2
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/mpfr-impl.h54
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/mpfr.h4
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/powerof2.c13
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/round_p.c10
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/round_prec.c308
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sub1.c14
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/sub_ui.c1
-rw-r--r--Build/source/libs/mpfr/mpfr-src/src/version.c2
11 files changed, 314 insertions, 97 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/Makefile.am b/Build/source/libs/mpfr/mpfr-src/src/Makefile.am
index fbd077710b0..1d28bf1ad0e 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/Makefile.am
+++ b/Build/source/libs/mpfr/mpfr-src/src/Makefile.am
@@ -75,7 +75,7 @@ libmpfr_la_LIBADD = @LIBOBJS@
# 2.4.x 3:x:2
# 3.0.x 4:x:0
# 3.1.x 5:x:1
-libmpfr_la_LDFLAGS = $(MPFR_LDFLAGS) $(LIBMPFR_LDFLAGS) -version-info 5:4:1
+libmpfr_la_LDFLAGS = $(MPFR_LDFLAGS) $(LIBMPFR_LDFLAGS) -version-info 5:5:1
# Important note: If for some reason, srcdir is read-only at build time
# (and you use objdir != srcdir), then you need to rebuild get_patches.c
diff --git a/Build/source/libs/mpfr/mpfr-src/src/add_ui.c b/Build/source/libs/mpfr/mpfr-src/src/add_ui.c
index 49c6752a074..700ca6f943c 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/add_ui.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/add_ui.c
@@ -49,6 +49,7 @@ mpfr_add_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode
MPFR_SAVE_EXPO_MARK (expo);
MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
inex = mpfr_add(y, x, uu, rnd_mode);
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range(y, inex, rnd_mode);
}
diff --git a/Build/source/libs/mpfr/mpfr-src/src/get_ld.c b/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
index 1c1efedbc4d..c7c5cd550ef 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/get_ld.c
@@ -122,7 +122,7 @@ mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
inex = mpfr_set (tmp, x, rnd_mode);
- mpfr_set_emin (-16382-63);
+ mpfr_set_emin (-16381-63); /* emin=-16444 */
mpfr_set_emax (16384);
mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
diff --git a/Build/source/libs/mpfr/mpfr-src/src/mpfr-impl.h b/Build/source/libs/mpfr/mpfr-src/src/mpfr-impl.h
index cce39c3824f..0066f90d360 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/mpfr-impl.h
+++ b/Build/source/libs/mpfr/mpfr-src/src/mpfr-impl.h
@@ -251,19 +251,6 @@ typedef struct __gmpfr_cache_s *mpfr_cache_ptr;
# define MPFR_WIN_THREAD_SAFE_DLL 1
#endif
-/* Detect some possible inconsistencies under Unix. */
-#if defined(__unix__)
-# if defined(_WIN32)
-# error "Both __unix__ and _WIN32 are defined"
-# endif
-# if __GMP_LIBGMP_DLL
-# error "__unix__ is defined and __GMP_LIBGMP_DLL is true"
-# endif
-# if defined(MPFR_WIN_THREAD_SAFE_DLL)
-# error "Both __unix__ and MPFR_WIN_THREAD_SAFE_DLL are defined"
-# endif
-#endif
-
#if defined(__MPFR_WITHIN_MPFR) || !defined(MPFR_WIN_THREAD_SAFE_DLL)
extern MPFR_THREAD_ATTR unsigned int __gmpfr_flags;
extern MPFR_THREAD_ATTR mpfr_exp_t __gmpfr_emin;
@@ -861,8 +848,8 @@ typedef intmax_t mpfr_eexp_t;
(MPFR_ASSERTD((s) == MPFR_SIGN_POS || (s) == MPFR_SIGN_NEG))
#define MPFR_SET_SIGN(x, s) \
(MPFR_ASSERT_SIGN(s), MPFR_SIGN(x) = s)
-#define MPFR_IS_POS_SIGN(s1) (s1 > 0)
-#define MPFR_IS_NEG_SIGN(s1) (s1 < 0)
+#define MPFR_IS_POS_SIGN(s1) ((s1) > 0)
+#define MPFR_IS_NEG_SIGN(s1) ((s1) < 0)
#define MPFR_MULT_SIGN(s1, s2) ((s1) * (s2))
/* Transform a sign to 1 or -1 */
#define MPFR_FROM_SIGN_TO_INT(s) (s)
@@ -914,33 +901,37 @@ typedef intmax_t mpfr_eexp_t;
/* We want to test if rnd = Zero, or Away.
'test' is 1 if negative, and 0 if positive. */
#define MPFR_IS_LIKE_RNDZ(rnd, test) \
- ((rnd==MPFR_RNDZ) || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, test))
+ ((rnd) == MPFR_RNDZ || MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, test))
-#define MPFR_IS_LIKE_RNDU(rnd, sign) \
- ((rnd==MPFR_RNDU) || (rnd==MPFR_RNDZ && sign<0) || (rnd==MPFR_RNDA && sign>0))
+#define MPFR_IS_LIKE_RNDU(rnd, sign) \
+ (((rnd) == MPFR_RNDU) || \
+ ((rnd) == MPFR_RNDZ && MPFR_IS_NEG_SIGN (sign)) || \
+ ((rnd) == MPFR_RNDA && MPFR_IS_POS_SIGN (sign)))
-#define MPFR_IS_LIKE_RNDD(rnd, sign) \
- ((rnd==MPFR_RNDD) || (rnd==MPFR_RNDZ && sign>0) || (rnd==MPFR_RNDA && sign<0))
+#define MPFR_IS_LIKE_RNDD(rnd, sign) \
+ (((rnd) == MPFR_RNDD) || \
+ ((rnd) == MPFR_RNDZ && MPFR_IS_POS_SIGN (sign)) || \
+ ((rnd) == MPFR_RNDA && MPFR_IS_NEG_SIGN (sign)))
/* Invert a rounding mode, RNDN, RNDZ and RNDA are unchanged */
-#define MPFR_INVERT_RND(rnd) ((rnd == MPFR_RNDU) ? MPFR_RNDD : \
- ((rnd == MPFR_RNDD) ? MPFR_RNDU : rnd))
+#define MPFR_INVERT_RND(rnd) ((rnd) == MPFR_RNDU ? MPFR_RNDD : \
+ (rnd) == MPFR_RNDD ? MPFR_RNDU : (rnd))
/* Transform RNDU and RNDD to RNDZ according to test */
-#define MPFR_UPDATE_RND_MODE(rnd, test) \
- do { \
- if (MPFR_UNLIKELY(MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test))) \
+#define MPFR_UPDATE_RND_MODE(rnd, test) \
+ do { \
+ if (MPFR_UNLIKELY(MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd, test))) \
rnd = MPFR_RNDZ; \
} while (0)
/* Transform RNDU and RNDD to RNDZ or RNDA according to sign,
leave the other modes unchanged */
-#define MPFR_UPDATE2_RND_MODE(rnd, sign) \
- do { \
- if (rnd == MPFR_RNDU) \
- rnd = (sign > 0) ? MPFR_RNDA : MPFR_RNDZ; \
- else if (rnd == MPFR_RNDD) \
- rnd = (sign < 0) ? MPFR_RNDA : MPFR_RNDZ; \
+#define MPFR_UPDATE2_RND_MODE(rnd, sign) \
+ do { \
+ if (rnd == MPFR_RNDU) \
+ rnd = MPFR_IS_POS_SIGN (sign) ? MPFR_RNDA : MPFR_RNDZ; \
+ else if (rnd == MPFR_RNDD) \
+ rnd = MPFR_IS_NEG_SIGN (sign) ? MPFR_RNDA : MPFR_RNDZ; \
} while (0)
@@ -1897,6 +1888,7 @@ __MPFR_DECLSPEC mpfr_exp_t mpfr_ceil_mul _MPFR_PROTO ((mpfr_exp_t, int, int));
__MPFR_DECLSPEC int mpfr_exp_2 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t));
__MPFR_DECLSPEC int mpfr_exp_3 _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t));
__MPFR_DECLSPEC int mpfr_powerof2_raw _MPFR_PROTO ((mpfr_srcptr));
+__MPFR_DECLSPEC int mpfr_powerof2_raw2 (const mp_limb_t *, mp_size_t);
__MPFR_DECLSPEC int mpfr_pow_general _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,
mpfr_srcptr, mpfr_rnd_t, int, mpfr_save_expo_t *));
diff --git a/Build/source/libs/mpfr/mpfr-src/src/mpfr.h b/Build/source/libs/mpfr/mpfr-src/src/mpfr.h
index e6105c94796..b23b1d2a9c8 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/mpfr.h
+++ b/Build/source/libs/mpfr/mpfr-src/src/mpfr.h
@@ -26,8 +26,8 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
/* Define MPFR version number */
#define MPFR_VERSION_MAJOR 3
#define MPFR_VERSION_MINOR 1
-#define MPFR_VERSION_PATCHLEVEL 4
-#define MPFR_VERSION_STRING "3.1.4"
+#define MPFR_VERSION_PATCHLEVEL 5
+#define MPFR_VERSION_STRING "3.1.5"
/* Macros dealing with MPFR VERSION */
#define MPFR_VERSION_NUM(a,b,c) (((a) << 16L) | ((b) << 8) | (c))
diff --git a/Build/source/libs/mpfr/mpfr-src/src/powerof2.c b/Build/source/libs/mpfr/mpfr-src/src/powerof2.c
index f10ace0ad03..53f0551cd66 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/powerof2.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/powerof2.c
@@ -31,16 +31,17 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
int
mpfr_powerof2_raw (mpfr_srcptr x)
{
- mp_limb_t *xp;
- mp_size_t xn;
-
/* This is an internal function, and we may call it with some
wrong numbers (ie good mantissa but wrong flags or exp)
So we don't want to test if it is a pure FP.
MPFR_ASSERTN(MPFR_IS_PURE_FP(x)); */
- xp = MPFR_MANT(x);
- xn = (MPFR_PREC(x) - 1) / GMP_NUMB_BITS;
- if (xp[xn] != MPFR_LIMB_HIGHBIT)
+ return mpfr_powerof2_raw2 (MPFR_MANT(x), MPFR_LIMB_SIZE(x));
+}
+
+int
+mpfr_powerof2_raw2 (const mp_limb_t *xp, mp_size_t xn)
+{
+ if (xp[--xn] != MPFR_LIMB_HIGHBIT)
return 0;
while (xn > 0)
if (xp[--xn] != 0)
diff --git a/Build/source/libs/mpfr/mpfr-src/src/round_p.c b/Build/source/libs/mpfr/mpfr-src/src/round_p.c
index f19a66433f5..8e2c0bebdaa 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/round_p.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/round_p.c
@@ -35,10 +35,16 @@ mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mpfr_exp_t err0, mpfr_prec_t prec)
i1 = mpfr_round_p_2 (bp, bn, err0, prec);
- /* compare with mpfr_can_round_raw */
+ /* Note: since revision 10747, mpfr_can_round_raw is supposed to be always
+ correct, whereas mpfr_round_p_2 might return 0 in some cases where one
+ could round, for example with err0=67 and prec=54:
+ b = 1111101101010001100011111011100010100011101111011011101111111111
+ thus we cannot compare i1 and i2, we only can check that we don't have
+ i1 <> 0 and i2 = 0.
+ */
i2 = mpfr_can_round_raw (bp, bn, MPFR_SIGN_POS, err0,
MPFR_RNDN, MPFR_RNDZ, prec);
- if (i1 != i2)
+ if (i1 && (i2 == 0))
{
fprintf (stderr, "mpfr_round_p(%d) != mpfr_can_round(%d)!\n"
"bn = %lu, err0 = %ld, prec = %lu\nbp = ", i1, i2,
diff --git a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c b/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
index c99d63e55ef..e5512b439dd 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/round_prec.c
@@ -138,47 +138,173 @@ mpfr_can_round (mpfr_srcptr b, mpfr_exp_t err, mpfr_rnd_t rnd1,
}
int
-mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
+mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec)
{
- mpfr_prec_t err, prec0 = prec;
+ mpfr_prec_t prec2;
mp_size_t k, k1, tn;
int s, s1;
mp_limb_t cc, cc2;
mp_limb_t *tmp;
+ mp_limb_t cy = 0, tmp_hi;
+ int res;
MPFR_TMP_DECL(marker);
- MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
+ /* Since mpfr_can_round is a function in the API, use MPFR_ASSERTN.
+ The specification makes sense only for prec >= 1. */
+ MPFR_ASSERTN (prec >= 1);
- if (MPFR_UNLIKELY(err0 < 0 || (mpfr_uexp_t) err0 <= prec))
- return 0; /* can't round */
+ MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);
MPFR_ASSERT_SIGN(neg);
neg = MPFR_IS_NEG_SIGN(neg);
/* Transform RNDD and RNDU to Zero / Away */
- MPFR_ASSERTD((neg == 0) || (neg == 1));
+ MPFR_ASSERTD (neg == 0 || neg == 1);
if (rnd1 != MPFR_RNDN)
rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA;
if (rnd2 != MPFR_RNDN)
rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA;
+ /* For err < prec (+1 for rnd1=RNDN), we can never round correctly, since
+ the error is at least 2*ulp(b) >= ulp(round(b)).
+ However for err = prec (+1 for rnd1=RNDN), we can round correctly in some
+ rare cases where ulp(b) = 1/2*ulp(U) [see below for the definition of U],
+ which implies rnd1 = RNDZ or RNDN, and rnd2 = RNDA or RNDN. */
+
+ if (MPFR_UNLIKELY (err < prec + (rnd1 == MPFR_RNDN) ||
+ (err == prec + (rnd1 == MPFR_RNDN) &&
+ (rnd1 == MPFR_RNDA ||
+ rnd2 == MPFR_RNDZ))))
+ return 0; /* can't round */
+
+ /* As a consequence... */
+ MPFR_ASSERTD (err >= prec);
+
+ /* The bound c on the error |x-b| is: c = 2^(MPFR_EXP(b)-err) <= b/2.
+ * So, we now know that x and b have the same sign. By symmetry,
+ * assume x > 0 and b > 0. We have: L <= x <= U, where, depending
+ * on rnd1:
+ * MPFR_RNDN: L = b-c, U = b+c
+ * MPFR_RNDZ: L = b, U = b+c
+ * MPFR_RNDA: L = b-c, U = b
+ *
+ * We can round x iff round(L,prec,rnd2) = round(U,prec,rnd2).
+ */
+
if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
- { /* Then prec < PREC(b): we can round:
- (i) in rounding to the nearest iff err0 >= prec + 2
+ { /* Then prec > PREC(b): we can round:
+ (i) in rounding to the nearest as long as err >= prec + 2.
+ When err = prec + 1 and b is not a power
+ of two (so that a change of binade cannot occur), then one
+ can round to nearest thanks to the even rounding rule (in the
+ target precision prec, the significand of b ends with a 0).
+ When err = prec + 1 and b is a power of two, when rnd1 = RNDZ one
+ can round too.
(ii) in directed rounding mode iff rnd1 is compatible with rnd2
- and err0 >= prec + 1, unless b = 2^k and rnd1=rnd2=RNDA in
- which case we need err0 >= prec + 2. */
+ and err >= prec + 1, unless b = 2^k and rnd1 = RNDA or RNDN in
+ which case we need err >= prec + 2.
+ */
+ if ((rnd1 == rnd2 || rnd2 == MPFR_RNDN) && err >= prec + 1)
+ {
+ if (rnd1 != MPFR_RNDZ &&
+ err == prec + 1 &&
+ mpfr_powerof2_raw2 (bp, bn))
+ return 0;
+ else
+ return 1;
+ }
+ return 0;
+ }
+
+ /* now prec <= bn * GMP_NUMB_BITS */
+
+ if (MPFR_UNLIKELY (err > (mpfr_prec_t) bn * GMP_NUMB_BITS))
+ {
+ /* we distinguish the case where b is a power of two:
+ rnd1 rnd2 can round?
+ RNDZ RNDZ ok
+ RNDZ RNDA no
+ RNDZ RNDN ok
+ RNDA RNDZ no
+ RNDA RNDA ok except when err = prec + 1
+ RNDA RNDN ok except when err = prec + 1
+ RNDN RNDZ no
+ RNDN RNDA no
+ RNDN RNDN ok except when err = prec + 1 */
+ if (mpfr_powerof2_raw2 (bp, bn))
+ {
+ if ((rnd2 == MPFR_RNDZ || rnd2 == MPFR_RNDA) && rnd1 != rnd2)
+ return 0;
+ else if (rnd1 == MPFR_RNDZ)
+ return 1; /* RNDZ RNDZ and RNDZ RNDN */
+ else
+ return err > prec + 1;
+ }
+
+ /* now the general case where b is not a power of two:
+ rnd1 rnd2 can round?
+ RNDZ RNDZ ok
+ RNDZ RNDA except when b is representable in precision 'prec'
+ RNDZ RNDN except when b is the middle of two representable numbers in
+ precision 'prec' and b ends with 'xxx0[1]',
+ or b is representable in precision 'prec'
+ and err = prec + 1 and b ends with '1'.
+ RNDA RNDZ except when b is representable in precision 'prec'
+ RNDA RNDA ok
+ RNDA RNDN except when b is the middle of two representable numbers in
+ precision 'prec' and b ends with 'xxx1[1]',
+ or b is representable in precision 'prec'
+ and err = prec + 1 and b ends with '1'.
+ RNDN RNDZ except when b is representable in precision 'prec'
+ RNDN RNDA except when b is representable in precision 'prec'
+ RNDN RNDN except when b is the middle of two representable numbers in
+ precision 'prec', or b is representable in precision 'prec'
+ and err = prec + 1 and b ends with '1'. */
if (rnd2 == MPFR_RNDN)
- return (mpfr_uexp_t) err0 - 2 >= prec;
+ {
+ if (err == prec + 1 && (bp[0] & 1))
+ return 0; /* err == prec + 1 implies prec = bn * GMP_NUMB_BITS */
+ if (prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
+ {
+ k1 = MPFR_PREC2LIMBS (prec + 1);
+ MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
+ if (((bp[bn - k1] >> s1) & 1) &&
+ mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
+ { /* b is the middle of two representable numbers */
+ if (rnd1 == MPFR_RNDN)
+ return 0;
+ k1 = MPFR_PREC2LIMBS (prec);
+ MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
+ return (rnd1 == MPFR_RNDZ) ^
+ (((bp[bn - k1] >> s1) & 1) == 0);
+ }
+ }
+ return 1;
+ }
+ else if (rnd1 == rnd2)
+ {
+ if (rnd1 == MPFR_RNDN && prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
+ {
+ /* then rnd2 = RNDN, and for prec = bn * GMP_NUMB_BITS we cannot
+ have b the middle of two representable numbers */
+ k1 = MPFR_PREC2LIMBS (prec + 1);
+ MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
+ if (((bp[bn - k1] >> s1) & 1) &&
+ mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
+ /* b is representable in precision prec+1 and ends with a 1 */
+ return 0;
+ else
+ return 1;
+ }
+ else
+ return 1;
+ }
else
- return (rnd1 == rnd2) && (mpfr_uexp_t) err0 - 2 >= prec;
+ return mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) != 0;
}
- /* if the error is smaller than ulp(b), then anyway it will propagate
- up to ulp(b) */
- err = ((mpfr_uexp_t) err0 > (mpfr_prec_t) bn * GMP_NUMB_BITS) ?
- (mpfr_prec_t) bn * GMP_NUMB_BITS : (mpfr_prec_t) err0;
+ /* now err <= bn * GMP_NUMB_BITS */
/* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */
k = (err - 1) / GMP_NUMB_BITS;
@@ -196,7 +322,7 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
Warning! The number with updated bn may no longer be normalized. */
k -= k1;
bn -= k1;
- prec -= (mpfr_prec_t) k1 * GMP_NUMB_BITS;
+ prec2 = prec - (mpfr_prec_t) k1 * GMP_NUMB_BITS;
/* We can decide of the correct rounding if rnd2(b-eps) and rnd2(b+eps)
give the same result to the target precision 'prec', i.e., if when
@@ -215,64 +341,150 @@ mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
switch (rnd1)
{
case MPFR_RNDZ:
- /* Round to Zero */
+ /* rnd1 = Round to Zero */
cc = (bp[bn - 1] >> s1) & 1;
/* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
and 0 otherwise */
- cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
+ cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
/* cc is the new value of bit s1 in bp[bn-1] after rounding 'rnd2' */
/* now round b + 2^(MPFR_EXP(b)-err) */
- mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
- /* if there was a carry here, then necessarily bit s1 of bp[bn-1]
- changed, thus we surely cannot round for directed rounding, but this
- will be detected below, with cc2 != cc */
+ cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+ /* propagate carry up to most significant limb */
+ for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
+ cy = ~bp[bn + tn] == 0;
+ if (cy == 0 && err == prec)
+ {
+ res = 0;
+ goto end;
+ }
+ if (MPFR_UNLIKELY(cy))
+ {
+ /* when a carry occurs, we have b < 2^h <= b+c, we can round iff:
+ rnd2 = RNDZ: never, since b and b+c round to different values;
+ rnd2 = RNDA: when b+c is an exact power of two, and err > prec
+ (since for err = prec, b = 2^h - 1/2*ulp(2^h) is
+ exactly representable and thus rounds to itself);
+ rnd2 = RNDN: whenever cc = 0, since err >= prec implies
+ c <= ulp(b) = 1/2*ulp(2^h), thus b+c rounds to 2^h,
+ and b+c >= 2^h implies that bit 'prec' of b is 1,
+ thus cc = 0 means that b is rounded to 2^h too. */
+ res = (rnd2 == MPFR_RNDZ) ? 0
+ : (rnd2 == MPFR_RNDA) ? (err > prec && k == bn && tmp[0] == 0)
+ : cc == 0;
+ goto end;
+ }
break;
case MPFR_RNDN:
- /* Round to nearest */
+ /* rnd1 = Round to nearest */
/* first round b+2^(MPFR_EXP(b)-err) */
- mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
- /* same remark as above in case a carry occurs in mpn_add_1() */
+ cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+ /* propagate carry up to most significant limb */
+ for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
+ cy = ~bp[bn + tn] == 0;
cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
- cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
+ cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2);
/* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
-
+ if (MPFR_UNLIKELY (cy != 0))
+ {
+ /* when a carry occurs, we have b-c < b < 2^h <= b+c, we can round
+ iff:
+ rnd2 = RNDZ: never, since b-c and b+c round to different values;
+ rnd2 = RNDA: when b+c is an exact power of two, and
+ err > prec + 1 (since for err <= prec + 1,
+ b-c <= 2^h - 1/2*ulp(2^h) is exactly representable
+ and thus rounds to itself);
+ rnd2 = RNDN: whenever err > prec + 1, since for err = prec + 1,
+ b+c rounds to 2^h, and b-c rounds to nextbelow(2^h).
+ For err > prec + 1, c <= 1/4*ulp(b) <= 1/8*ulp(2^h),
+ thus
+ 2^h - 1/4*ulp(b) <= b-c < b+c <= 2^h + 1/8*ulp(2^h),
+ therefore both b-c and b+c round to 2^h. */
+ res = (rnd2 == MPFR_RNDZ) ? 0
+ : (rnd2 == MPFR_RNDA) ? (err > prec + 1 && k == bn && tmp[0] == 0)
+ : err > prec + 1;
+ goto end;
+ }
subtract_eps:
- /* now round b-2^(MPFR_EXP(b)-err) */
- cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
+ /* now round b-2^(MPFR_EXP(b)-err), this happens for
+ rnd1 = RNDN or RNDA */
+ MPFR_ASSERTD(rnd1 == MPFR_RNDN || rnd1 == MPFR_RNDA);
+ cy = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
/* propagate the potential borrow up to the most significant limb
(it cannot propagate further since the most significant limb is
- at least MPFR_LIMB_HIGHBIT) */
- for (tn = 0; tn + 1 < k1 && (cc2 != 0); tn ++)
- cc2 = bp[bn + tn] == 0;
- /* We have an exponent decrease when either:
- (i) k1 = 0 and tmp[bn-1] < MPFR_LIMB_HIGHBIT
- (ii) k1 > 0 and cc <> 0 and bp[bn + tn] = MPFR_LIMB_HIGHBIT
- (then necessarily tn = k1-1).
- Then for directed rounding we cannot round,
- and for rounding to nearest we cannot round when err = prec + 1.
+ at least MPFR_LIMB_HIGHBIT).
+ Note: we use the same limb tmp[bn-1] to subtract. */
+ tmp_hi = tmp[bn - 1];
+ for (tn = 0; tn < k1 && cy != 0; tn ++)
+ cy = mpn_sub_1 (&tmp_hi, bp + bn + tn, 1, cy);
+ /* We have an exponent decrease when tn = k1 and
+ tmp[bn-1] < MPFR_LIMB_HIGHBIT:
+ b-c < 2^h <= b (for RNDA) or b+c (for RNDN).
+ Then we surely cannot round when rnd2 = RNDZ, since b or b+c round to
+ a value >= 2^h, and b-c rounds to a value < 2^h.
+ We also surely cannot round when (rnd1,rnd2) = (RNDN,RNDA), since
+ b-c rounds to a value <= 2^h, and b+c > 2^h rounds to a value > 2^h.
+ It thus remains:
+ (rnd1,rnd2) = (RNDA,RNDA), (RNDA,RNDN) and (RNDN,RNDN).
+ For (RNDA,RNDA) we can round only when b-c and b round to 2^h, which
+ implies b = 2^h and err > prec (which is true in that case):
+ a necessary condition is that cc = 0.
+ For (RNDA,RNDN) we can round only when b-c and b round to 2^h, which
+ implies b-c >= 2^h - 1/4*ulp(2^h), and b <= 2^h + 1/2*ulp(2^h);
+ since ulp(2^h) = ulp(b), this implies c <= 3/4*ulp(b), thus
+ err > prec.
+ For (RNDN,RNDN) we can round only when b-c and b+c round to 2^h,
+ which implies b-c >= 2^h - 1/4*ulp(2^h), and
+ b+c <= 2^h + 1/2*ulp(2^h);
+ since ulp(2^h) = ulp(b), this implies 2*c <= 3/4*ulp(b), thus
+ err > prec+1.
*/
- if (((k1 == 0 && tmp[bn - 1] < MPFR_LIMB_HIGHBIT) ||
- (k1 != 0 && cc2 != 0 && bp[bn + tn] == MPFR_LIMB_HIGHBIT)) &&
- (rnd2 != MPFR_RNDN || err0 == prec0 + 1))
+ if (tn == k1 && tmp_hi < MPFR_LIMB_HIGHBIT) /* exponent decrease */
+ {
+ if (rnd2 == MPFR_RNDZ || (rnd1 == MPFR_RNDN && rnd2 == MPFR_RNDA) ||
+ cc != 0 /* b or b+c does not round to 2^h */)
+ {
+ res = 0;
+ goto end;
+ }
+ /* in that case since the most significant bit of tmp is 0, we
+ should consider one more bit; res = 0 when b-c does not round
+ to 2^h. */
+ res = mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2 + 1) != 0;
+ goto end;
+ }
+ if (err == prec + (rnd1 == MPFR_RNDN))
{
- MPFR_TMP_FREE(marker);
- return 0;
+ /* No exponent increase nor decrease, thus we have |U-L| = ulp(b).
+ For rnd2 = RNDZ or RNDA, either [L,U] contains one representable
+ number in the target precision, and then L and U round
+ differently; or both L and U are representable: they round
+ differently too; thus in all cases we cannot round.
+ For rnd2 = RNDN, the only case where we can round is when the
+ middle of [L,U] (i.e. b) is representable, and ends with a 0. */
+ res = (rnd2 == MPFR_RNDN && (((bp[bn - 1] >> s1) & 1) == 0) &&
+ mpfr_round_raw2 (bp, bn, neg, MPFR_RNDZ, prec2) ==
+ mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec2));
+ goto end;
}
break;
default:
- /* Round away */
+ /* rnd1 = Round away */
+ MPFR_ASSERTD (rnd1 == MPFR_RNDA);
cc = (bp[bn - 1] >> s1) & 1;
- cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
+ /* the mpfr_round_raw2() call below returns whether one should add 1 or
+ not for rounding */
+ cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
/* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
goto subtract_eps;
}
cc2 = (tmp[bn - 1] >> s1) & 1;
- cc2 ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
+ res = cc == (cc2 ^ mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2));
+ end:
MPFR_TMP_FREE(marker);
- return cc == cc2;
+ return res;
}
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sub1.c b/Build/source/libs/mpfr/mpfr-src/src/sub1.c
index b00c4cd7cfc..cf361393876 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sub1.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sub1.c
@@ -96,16 +96,15 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* A = S*ABS(B) +/- ulp(a) */
MPFR_SET_EXP (a, MPFR_GET_EXP (b));
MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), MPFR_PREC (b),
- rnd_mode, MPFR_SIGN (a),
- if (MPFR_UNLIKELY ( ++MPFR_EXP (a) > __gmpfr_emax))
- inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a)));
- /* inexact = mpfr_set4 (a, b, rnd_mode, MPFR_SIGN (a)); */
+ rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a));
if (inexact == 0)
{
/* a = b (Exact)
But we know it isn't (Since we have to remove `c')
So if we round to Zero, we have to remove one ulp.
Otherwise the result is correctly rounded. */
+ /* An overflow is not possible. */
+ MPFR_ASSERTD (MPFR_EXP (a) <= __gmpfr_emax);
if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
{
mpfr_nexttozero (a);
@@ -136,9 +135,14 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
i.e. inexact= MPFR_EVEN_INEX */
if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a)))
{
- mpfr_nexttozero (a);
+ if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+ mpfr_setmax (a, __gmpfr_emax);
+ else
+ mpfr_nexttozero (a);
inexact = -MPFR_INT_SIGN (a);
}
+ else if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+ inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
MPFR_RET (inexact);
}
}
diff --git a/Build/source/libs/mpfr/mpfr-src/src/sub_ui.c b/Build/source/libs/mpfr/mpfr-src/src/sub_ui.c
index b366970eb34..718db0a43aa 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/sub_ui.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/sub_ui.c
@@ -52,6 +52,7 @@ mpfr_sub_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode
MPFR_SAVE_EXPO_MARK (expo);
MPFR_SET_EXP (uu, GMP_NUMB_BITS - cnt);
inex = mpfr_sub (y, x, uu, rnd_mode);
+ MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
MPFR_SAVE_EXPO_FREE (expo);
return mpfr_check_range (y, inex, rnd_mode);
}
diff --git a/Build/source/libs/mpfr/mpfr-src/src/version.c b/Build/source/libs/mpfr/mpfr-src/src/version.c
index 130e0f86c1d..0e24a47f022 100644
--- a/Build/source/libs/mpfr/mpfr-src/src/version.c
+++ b/Build/source/libs/mpfr/mpfr-src/src/version.c
@@ -25,5 +25,5 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
const char *
mpfr_get_version (void)
{
- return "3.1.4";
+ return "3.1.5";
}