summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile.am2
-rw-r--r--src/add.c20
-rw-r--r--src/add1.c29
-rw-r--r--src/cmp2.c31
-rw-r--r--src/fmma.c139
-rw-r--r--src/mpfr-impl.h80
-rw-r--r--src/mpfr.h2
-rw-r--r--src/print_raw.c15
-rw-r--r--src/sub.c24
-rw-r--r--src/sub1.c119
-rw-r--r--src/ubf.c226
-rw-r--r--src/vasprintf.c18
-rw-r--r--tests/tacosh.c11
-rw-r--r--tests/texp.c31
-rw-r--r--tests/tfmma.c283
-rw-r--r--tests/tpow.c7
-rw-r--r--tests/tpow_z.c30
-rw-r--r--tests/tsprintf.c8
-rw-r--r--tests/tsub1sp.c7
19 files changed, 799 insertions, 283 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index 3c80ea743..d05085a7e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,7 +59,7 @@ buildopt.c digamma.c bernoulli.c isregular.c set_flt.c get_flt.c \
scale2.c set_z_exp.c ai.c gammaonethird.c ieee_floats.h \
grandom.c fpif.c set_float128.c get_float128.c rndna.c nrandom.c \
random_deviate.h random_deviate.c erandom.c mpfr-mini-gmp.c \
-mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c
+mpfr-mini-gmp.h fmma.c log_ui.c gamma_inc.c ubf.c
libmpfr_la_LIBADD = @LIBOBJS@
diff --git a/src/add.c b/src/add.c
index a953952d2..969a64d0a 100644
--- a/src/add.c
+++ b/src/add.c
@@ -31,7 +31,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode),
("a[%Pu]=%.*Rg", mpfr_get_prec (a), mpfr_log_prec, a));
- if (MPFR_ARE_SINGULAR(b,c))
+ if (MPFR_ARE_SINGULAR_OR_UBF (b, c))
{
if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c))
{
@@ -59,7 +59,7 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
MPFR_SET_SAME_SIGN(a, c);
MPFR_RET(0); /* exact */
}
- /* now either b or c is zero */
+ /* now both b and c are finite numbers */
else if (MPFR_IS_ZERO(b))
{
if (MPFR_IS_ZERO(c))
@@ -78,11 +78,23 @@ mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
}
return mpfr_set (a, c, rnd_mode);
}
- else
+ else if (MPFR_IS_ZERO(c))
{
- MPFR_ASSERTD(MPFR_IS_ZERO(c));
return mpfr_set (a, b, rnd_mode);
}
+ else
+ {
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
+ /* mpfr_sub1sp and mpfr_add1sp are not intended to support UBF,
+ for which optimization is less important. */
+ if (MPFR_SIGN(b) != MPFR_SIGN(c))
+ return mpfr_sub1 (a, b, c, rnd_mode);
+ else if (MPFR_EXP_LESS_P (b, c))
+ return mpfr_add1 (a, c, b, rnd_mode);
+ else
+ return mpfr_add1 (a, b, c, rnd_mode);
+ }
}
MPFR_ASSERTD (MPFR_IS_PURE_FP (b));
diff --git a/src/add1.c b/src/add1.c
index 671836899..4e210390e 100644
--- a/src/add1.c
+++ b/src/add1.c
@@ -30,13 +30,21 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mp_limb_t *ap, *bp, *cp;
mpfr_prec_t aq, bq, cq, aq2;
mp_size_t an, bn, cn;
- mpfr_exp_t difw, exp;
+ mpfr_exp_t difw, exp, diff_exp;
int sh, rb, fb, inex;
- mpfr_uexp_t diff_exp;
MPFR_TMP_DECL(marker);
- MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
- MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
+
+ if (MPFR_UNLIKELY (MPFR_IS_UBF (b)))
+ {
+ exp = mpfr_ubf_zexp2exp (MPFR_ZEXP (b));
+ if (exp > __gmpfr_emax)
+ return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));;
+ }
+ else
+ exp = MPFR_GET_EXP (b);
MPFR_TMP_MARK(marker);
@@ -68,13 +76,18 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
MPN_COPY(cp, ap, cn);
}
- exp = MPFR_GET_EXP (b);
MPFR_SET_SAME_SIGN(a, b);
MPFR_UPDATE2_RND_MODE(rnd_mode, MPFR_SIGN(b));
/* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
- /* Note: exponents can be negative, but the unsigned subtraction is
- a modular subtraction, so that one gets the correct result. */
- diff_exp = (mpfr_uexp_t) exp - MPFR_GET_EXP(c);
+ if (MPFR_UNLIKELY (MPFR_IS_UBF (c)))
+ {
+ MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX);
+ diff_exp = mpfr_ubf_diff_exp (b, c);
+ }
+ else
+ diff_exp = exp - MPFR_GET_EXP (c);
+
+ MPFR_ASSERTD (diff_exp >= 0);
/*
* 1. Compute the significant part A', the non-significant bits of A
diff --git a/src/cmp2.c b/src/cmp2.c
index c4eee67f7..9c649b822 100644
--- a/src/cmp2.c
+++ b/src/cmp2.c
@@ -38,6 +38,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
{
mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
mp_size_t bn, cn;
+ mpfr_exp_t sdiff_exp;
mpfr_uexp_t diff_exp;
mpfr_prec_t res = 0;
int sign;
@@ -48,13 +49,21 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
/* the cases b=0 or c=0 are also treated apart in agm and sub
(which calls sub1) */
- MPFR_ASSERTD (MPFR_IS_PURE_FP(b));
- MPFR_ASSERTD (MPFR_IS_PURE_FP(c));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
- if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c))
+ sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ?
+ mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
+
+ /* The returned result is restricted to [MPFR_EXP_MIN,MPFR_EXP_MAX],
+ which is the range of the mpfr_exp_t type. But under the condition
+ below, the value of cancel will not be affected. */
+ MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX);
+
+ if (sdiff_exp >= 0)
{
sign = 1;
- diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
+ diff_exp = sdiff_exp;
bp = MPFR_MANT(b);
cp = MPFR_MANT(c);
@@ -85,7 +94,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
if (MPFR_UNLIKELY (cn < 0))
/* c discards exactly the upper part of b */
{
- unsigned int z;
+ int z;
MPFR_ASSERTD (bn >= 0);
@@ -96,7 +105,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
res += GMP_NUMB_BITS;
}
- count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
+ count_leading_zeros (z, bp[bn]); /* bp[bn] <> 0 */
*cancel = res + z;
return sign;
}
@@ -118,7 +127,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
else /* MPFR_EXP(b) < MPFR_EXP(c) */
{
sign = -1;
- diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b);
+ diff_exp = - (mpfr_uexp_t) sdiff_exp;
bp = MPFR_MANT(c);
cp = MPFR_MANT(b);
@@ -152,6 +161,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
&& (high_dif == 0) && (dif == 1)))
{ /* dif=1 implies diff_exp = 0 or 1 */
+ MPFR_ASSERTD (diff_exp <= 1);
bb = (bn >= 0) ? bp[bn] : 0;
cc = lastc;
if (cn >= 0)
@@ -160,8 +170,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
{
cc += cp[cn];
}
- else /* diff_exp = 1 */
+ else
{
+ MPFR_ASSERTD (diff_exp == 1);
cc += cp[cn] >> 1;
lastc = cp[cn] << (GMP_NUMB_BITS - 1);
}
@@ -188,9 +199,9 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
}
else /* high_dif == 0 */
{
- unsigned int z;
+ int z;
- count_leading_zeros(z, dif); /* dif > 1 here */
+ count_leading_zeros (z, dif); /* dif > 1 here */
res += z;
if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (GMP_NUMB_BITS - z - 1))))
{ /* dif is not a power of two */
diff --git a/src/fmma.c b/src/fmma.c
index 9f0af2faf..fbd94ff6d 100644
--- a/src/fmma.c
+++ b/src/fmma.c
@@ -23,103 +23,43 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#include "mpfr-impl.h"
static int
-mpfr_fmma_slow (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
- mpfr_srcptr d, mpfr_rnd_t rnd)
+mpfr_fmma_aux (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
+ mpfr_srcptr d, mpfr_rnd_t rnd, int neg)
{
- mpfr_t ab, cd;
+ mpfr_ubf_t u, v;
+ mp_size_t un, vn;
+ mpfr_limb_ptr up, vp;
int inex;
+ MPFR_TMP_DECL(marker);
- mpfr_init2 (ab, MPFR_PREC(a) + MPFR_PREC(b));
- mpfr_init2 (cd, MPFR_PREC(c) + MPFR_PREC(d));
- /* FIXME: The following multiplications may overflow or underflow
- (even more often with the fact that the exponent range is not
- extended), in which case the result is not exact. This should
- be solved with future unbounded floats. */
- mpfr_mul (ab, a, b, MPFR_RNDZ); /* exact */
- mpfr_mul (cd, c, d, MPFR_RNDZ); /* exact */
- inex = mpfr_add (z, ab, cd, rnd);
- mpfr_clear (ab);
- mpfr_clear (cd);
- return inex;
-}
+ MPFR_LOG_FUNC
+ (("a[%Pu]=%.*Rg b[%Pu]=%.*Rg c[%Pu]=%.*Rg d[%Pu]=%.*Rg rnd=%d neg=%d",
+ mpfr_get_prec (a), mpfr_log_prec, a,
+ mpfr_get_prec (b), mpfr_log_prec, b,
+ mpfr_get_prec (c), mpfr_log_prec, c,
+ mpfr_get_prec (d), mpfr_log_prec, d, rnd, neg),
+ ("z[%Pu]=%.*Rg inex=%d",
+ mpfr_get_prec (z), mpfr_log_prec, z, inex));
-/* z <- a*b + c*d */
-static int
-mpfr_fmma_fast (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
- mpfr_srcptr d, mpfr_rnd_t rnd)
-{
- /* Assumes that a, b, c, d are finite and non-zero; so any multiplication
- of two of them yielding an infinity is an overflow, and a
- multiplication yielding 0 is an underflow.
- Assumes further that z is distinct from a, b, c, d. */
-
- int inex;
- mpfr_t u, v;
- mp_size_t an, bn, cn, dn;
- mpfr_limb_ptr up, vp;
- MPFR_TMP_DECL(marker);
- MPFR_SAVE_EXPO_DECL (expo);
-
- MPFR_TMP_MARK(marker);
- MPFR_SAVE_EXPO_MARK (expo);
-
- /* u=a*b, v=c*d exactly */
- an = MPFR_LIMB_SIZE(a);
- bn = MPFR_LIMB_SIZE(b);
- cn = MPFR_LIMB_SIZE(c);
- dn = MPFR_LIMB_SIZE(d);
- MPFR_TMP_INIT (up, u, (an + bn) * GMP_NUMB_BITS, an + bn);
- MPFR_TMP_INIT (vp, v, (cn + dn) * GMP_NUMB_BITS, cn + dn);
-
- /* u <- a*b */
- if (an >= bn)
- mpn_mul (up, MPFR_MANT(a), an, MPFR_MANT(b), bn);
- else
- mpn_mul (up, MPFR_MANT(b), bn, MPFR_MANT(a), an);
- if ((up[an + bn - 1] & MPFR_LIMB_HIGHBIT) == 0)
- {
- mpn_lshift (up, up, an + bn, 1);
- /* EXP(a) and EXP(b) are in [1-2^(n-2), 2^(n-2)-1] where
- mpfr_exp_t has is n-bit wide, thus EXP(a)+EXP(b) is in
- [2-2^(n-1), 2^(n-1)-2]. We use the fact here that 2*MPFR_EMIN_MIN-1
- is a valid exponent (see mpfr-impl.h). However we don't use
- MPFR_SET_EXP() which only allows the restricted exponent range
- [1-2^(n-2), 2^(n-2)-1]. */
- MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b) - 1;
- }
- else
- MPFR_EXP(u) = MPFR_EXP(a) + MPFR_EXP(b);
-
- /* v <- c*d */
- if (cn >= dn)
- mpn_mul (vp, MPFR_MANT(c), cn, MPFR_MANT(d), dn);
- else
- mpn_mul (vp, MPFR_MANT(d), dn, MPFR_MANT(c), cn);
- if ((vp[cn + dn - 1] & MPFR_LIMB_HIGHBIT) == 0)
- {
- mpn_lshift (vp, vp, cn + dn, 1);
- MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d) - 1;
- }
- else
- MPFR_EXP(v) = MPFR_EXP(c) + MPFR_EXP(d);
-
- MPFR_PREC(u) = (an + bn) * GMP_NUMB_BITS;
- MPFR_PREC(v) = (cn + dn) * GMP_NUMB_BITS;
- MPFR_SIGN(u) = MPFR_MULT_SIGN(MPFR_SIGN(a), MPFR_SIGN(b));
- MPFR_SIGN(v) = MPFR_MULT_SIGN(MPFR_SIGN(c), MPFR_SIGN(d));
-
- /* tentatively compute z as u+v; here we need z to be distinct
- from a, b, c, d to avoid losing the input values in case we
- need to call mpfr_fmma_slow */
- /* FIXME: The above comment is no longer valid. Anyway, with
- unbounded floats (based on an exact multiplication like above),
- it will no longer be necessary to distinguish fast and slow. */
- inex = mpfr_add (z, u, v, rnd);
-
- MPFR_TMP_FREE(marker);
- MPFR_SAVE_EXPO_FREE (expo);
-
- return mpfr_check_range (z, inex, rnd);
+ MPFR_TMP_MARK (marker);
+
+ un = MPFR_LIMB_SIZE (a) + MPFR_LIMB_SIZE (b);
+ vn = MPFR_LIMB_SIZE (c) + MPFR_LIMB_SIZE (d);
+ MPFR_TMP_INIT (up, u, (mpfr_prec_t) un * GMP_NUMB_BITS, un);
+ MPFR_TMP_INIT (vp, v, (mpfr_prec_t) vn * GMP_NUMB_BITS, vn);
+
+ mpfr_ubf_mul_exact (u, a, b);
+ mpfr_ubf_mul_exact (v, c, d);
+ if (neg)
+ MPFR_CHANGE_SIGN (v);
+ inex = mpfr_add (z, (mpfr_srcptr) u, (mpfr_srcptr) v, rnd);
+
+ MPFR_UBF_CLEAR_EXP (u);
+ MPFR_UBF_CLEAR_EXP (v);
+
+ MPFR_TMP_FREE (marker);
+
+ return inex;
}
/* z <- a*b + c*d */
@@ -127,13 +67,7 @@ int
mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
mpfr_srcptr d, mpfr_rnd_t rnd)
{
- mpfr_limb_ptr zp = MPFR_MANT(z);
-
- return (mpfr_regular_p (a) && mpfr_regular_p (b) && mpfr_regular_p (c) &&
- mpfr_regular_p (d) && zp != MPFR_MANT(a) && zp != MPFR_MANT(b) &&
- zp != MPFR_MANT(c) && zp != MPFR_MANT(d))
- ? mpfr_fmma_fast (z, a, b, c, d, rnd)
- : mpfr_fmma_slow (z, a, b, c, d, rnd);
+ return mpfr_fmma_aux (z, a, b, c, d, rnd, 0);
}
/* z <- a*b - c*d */
@@ -141,8 +75,5 @@ int
mpfr_fmms (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c,
mpfr_srcptr d, mpfr_rnd_t rnd)
{
- mpfr_t minus_c;
-
- MPFR_ALIAS (minus_c, c, -MPFR_SIGN(c), MPFR_EXP(c));
- return mpfr_fmma (z, a, b, minus_c, d, rnd);
+ return mpfr_fmma_aux (z, a, b, c, d, rnd, 1);
}
diff --git a/src/mpfr-impl.h b/src/mpfr-impl.h
index 8e7ccaa60..cc035d346 100644
--- a/src/mpfr-impl.h
+++ b/src/mpfr-impl.h
@@ -840,18 +840,17 @@ typedef intmax_t mpfr_eexp_t;
/* Definition of the exponent limits for MPFR numbers.
* These limits are chosen so that if e is such an exponent, then 2e-1 and
- * 2e+1 are valid exponents. This is useful for intermediate computations,
- * in particular the multiplication. We must have MPFR_EMIN_MIN >= 3-2^(n-2)
- * = 3-MPFR_EXP_INVALID so that 2*MPFR_EMIN_MIN-1 > __MPFR_EXP_INF = 3-2^(n-1).
+ * 2e+1 are representable. This is useful for intermediate computations,
+ * in particular the multiplication.
*/
#undef MPFR_EMIN_MIN
#undef MPFR_EMIN_MAX
#undef MPFR_EMAX_MIN
#undef MPFR_EMAX_MAX
-#define MPFR_EMIN_MIN (3-MPFR_EXP_INVALID)
-#define MPFR_EMIN_MAX (MPFR_EXP_INVALID-3)
-#define MPFR_EMAX_MIN (3-MPFR_EXP_INVALID)
-#define MPFR_EMAX_MAX (MPFR_EXP_INVALID-3)
+#define MPFR_EMIN_MIN (1-MPFR_EXP_INVALID)
+#define MPFR_EMIN_MAX (MPFR_EXP_INVALID-1)
+#define MPFR_EMAX_MIN (1-MPFR_EXP_INVALID)
+#define MPFR_EMAX_MAX (MPFR_EXP_INVALID-1)
/* Use MPFR_GET_EXP and MPFR_SET_EXP instead of MPFR_EXP directly,
unless when the exponent may be out-of-range, for instance when
@@ -880,6 +879,10 @@ typedef intmax_t mpfr_eexp_t;
# define MPFR_SET_INVALID_EXP(x) ((void) 0)
#endif
+#define MPFR_EXP_LESS_P(x,y) \
+ (MPFR_UNLIKELY (MPFR_IS_UBF (x) || MPFR_IS_UBF (y)) ? \
+ mpfr_ubf_exp_less_p (x, y) : MPFR_GET_EXP (x) < MPFR_GET_EXP (y))
+
/******************************************************
********* Singular values (NAN, INF, ZERO) *********
@@ -889,6 +892,7 @@ typedef intmax_t mpfr_eexp_t;
# define MPFR_EXP_ZERO (MPFR_EXP_MIN+1)
# define MPFR_EXP_NAN (MPFR_EXP_MIN+2)
# define MPFR_EXP_INF (MPFR_EXP_MIN+3)
+# define MPFR_EXP_UBF (MPFR_EXP_MIN+4)
#define MPFR_IS_NAN(x) (MPFR_EXP(x) == MPFR_EXP_NAN)
#define MPFR_SET_NAN(x) (MPFR_EXP(x) = MPFR_EXP_NAN)
@@ -897,21 +901,34 @@ typedef intmax_t mpfr_eexp_t;
#define MPFR_IS_ZERO(x) (MPFR_EXP(x) == MPFR_EXP_ZERO)
#define MPFR_SET_ZERO(x) (MPFR_EXP(x) = MPFR_EXP_ZERO)
#define MPFR_NOTZERO(x) (MPFR_EXP(x) != MPFR_EXP_ZERO)
+#define MPFR_IS_UBF(x) (MPFR_EXP(x) == MPFR_EXP_UBF)
+#define MPFR_SET_UBF(x) (MPFR_EXP(x) = MPFR_EXP_UBF)
#define MPFR_IS_NORMALIZED(x) \
(MPFR_LIMB_MSB (MPFR_MANT(x)[MPFR_LAST_LIMB(x)]) != 0)
#define MPFR_IS_FP(x) (!MPFR_IS_NAN(x) && !MPFR_IS_INF(x))
#define MPFR_IS_SINGULAR(x) (MPFR_EXP(x) <= MPFR_EXP_INF)
+#define MPFR_IS_SINGULAR_OR_UBF(x) (MPFR_EXP(x) <= MPFR_EXP_UBF)
#define MPFR_IS_PURE_FP(x) \
(!MPFR_IS_SINGULAR(x) && \
(MPFR_ASSERTD (MPFR_EXP (x) >= MPFR_EMIN_MIN && \
MPFR_EXP (x) <= MPFR_EMAX_MAX && \
MPFR_IS_NORMALIZED (x)), 1))
+#define MPFR_IS_PURE_UBF(x) \
+ (!MPFR_IS_SINGULAR(x) && \
+ (MPFR_ASSERTD ((MPFR_IS_UBF (x) || \
+ (MPFR_EXP (x) >= MPFR_EMIN_MIN && \
+ MPFR_EXP (x) <= MPFR_EMAX_MAX)) && \
+ MPFR_IS_NORMALIZED (x)), 1))
#define MPFR_ARE_SINGULAR(x,y) \
(MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
+#define MPFR_ARE_SINGULAR_OR_UBF(x,y) \
+ (MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(x)) || \
+ MPFR_UNLIKELY(MPFR_IS_SINGULAR_OR_UBF(y)))
+
#define MPFR_IS_POWER_OF_2(x) \
(mpfr_cmp_ui_2exp ((x), 1, MPFR_GET_EXP (x) - 1) == 0)
@@ -2234,4 +2251,53 @@ __MPFR_DECLSPEC extern int __gmpfr_cov_sum_tmd[MPFR_RND_MAX][2][2][3][2][2];
#endif /* MPFR_COV_CHECK */
+/******************************************************
+ ***************** Unbounded Floats *****************
+ ******************************************************/
+
+#if defined (__cplusplus)
+extern "C" {
+#endif
+
+/* An UBF is like a MPFR number, but with an additional mpz_t member,
+ which is assumed to be present (with a value in it) when the usual
+ exponent field has the value MPFR_EXP_UBF. The goal of this compatible
+ representation is to easily be able to support UBF in "normal" code
+ and hopefully avoid aliasing issues at the same time. And code that
+ accepts UBF in input should also accept mpfr_t as a consequence; this
+ makes mpfr_t to UBF conversion unnecessary.
+ When an input of a public function is a UBF, the semantic remains
+ internal to MPFR and can change in the future.
+ Note that functions used for logging need to support UBF (currently
+ done by printing that a number is a UBF, as it may be difficult to
+ do more without significant changes). */
+
+typedef struct {
+ mpfr_prec_t _mpfr_prec;
+ mpfr_sign_t _mpfr_sign;
+ mpfr_exp_t _mpfr_exp;
+ mp_limb_t *_mpfr_d;
+ mpz_t _mpfr_zexp;
+} __mpfr_ubf_struct;
+
+typedef __mpfr_ubf_struct mpfr_ubf_t[1];
+typedef __mpfr_ubf_struct *mpfr_ubf_ptr;
+
+__MPFR_DECLSPEC void mpfr_ubf_mul_exact (mpfr_ubf_ptr,
+ mpfr_srcptr, mpfr_srcptr);
+__MPFR_DECLSPEC int mpfr_ubf_exp_less_p (mpfr_srcptr, mpfr_srcptr);
+__MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_zexp2exp (mpz_ptr);
+__MPFR_DECLSPEC mpfr_exp_t mpfr_ubf_diff_exp (mpfr_srcptr, mpfr_srcptr);
+
+#if defined (__cplusplus)
+}
+#endif
+
+#define MPFR_ZEXP(x) \
+ ((void) (x)->_mpfr_exp /* to check that x has a correct type */, \
+ ((mpfr_ubf_ptr) (x))->_mpfr_zexp)
+
+#define MPFR_UBF_CLEAR_EXP(x) \
+ ((void) (MPFR_IS_UBF (u) && (mpz_clear (MPFR_ZEXP (x)), 0)))
+
#endif /* __MPFR_IMPL_H__ */
diff --git a/src/mpfr.h b/src/mpfr.h
index 0fa453901..9c7a75f3f 100644
--- a/src/mpfr.h
+++ b/src/mpfr.h
@@ -187,7 +187,7 @@ typedef uintmax_t mpfr_uexp_t;
#endif
/* Definition of the standard exponent limits */
-#define MPFR_EMAX_DEFAULT ((mpfr_exp_t) (((mpfr_ulong) 1 << 30) - 3))
+#define MPFR_EMAX_DEFAULT ((mpfr_exp_t) (((mpfr_ulong) 1 << 30) - 1))
#define MPFR_EMIN_DEFAULT (-(MPFR_EMAX_DEFAULT))
/* DON'T USE THIS! (For MPFR-public macros only, see below.)
diff --git a/src/print_raw.c b/src/print_raw.c
index af9c6c1b9..cd19ddee7 100644
--- a/src/print_raw.c
+++ b/src/print_raw.c
@@ -49,25 +49,22 @@ mpfr_fprint_binary (FILE *stream, mpfr_srcptr x)
px = MPFR_PREC (x);
fprintf (stream, "0.");
- for (n = (px - 1) / GMP_NUMB_BITS; ; n--)
+ for (n = (px - 1) / GMP_NUMB_BITS; n >= 0; n--)
{
mp_limb_t wd, t;
- MPFR_ASSERTN (n >= 0);
wd = mx[n];
for (t = MPFR_LIMB_HIGHBIT; t != 0; t >>= 1)
{
putc ((wd & t) == 0 ? '0' : '1', stream);
if (--px == 0)
- {
- mpfr_exp_t ex;
-
- ex = MPFR_EXP (x);
- fprintf (stream, "E%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) ex);
- return;
- }
+ break;
}
}
+ if (MPFR_IS_UBF (x))
+ gmp_fprintf (stream, "E%Zd", MPFR_ZEXP (x));
+ else
+ fprintf (stream, "E%" MPFR_EXP_FSPEC "d", (mpfr_eexp_t) MPFR_EXP (x));
}
}
diff --git a/src/sub.c b/src/sub.c
index 457f67d9d..c95cfdb0a 100644
--- a/src/sub.c
+++ b/src/sub.c
@@ -31,7 +31,7 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpfr_get_prec (c), mpfr_log_prec, c, rnd_mode),
("a[%Pu]=%.*Rg", mpfr_get_prec (a), mpfr_log_prec, a));
- if (MPFR_ARE_SINGULAR (b,c))
+ if (MPFR_ARE_SINGULAR_OR_UBF (b,c))
{
if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
{
@@ -72,11 +72,29 @@ mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
else
return mpfr_neg (a, c, rnd_mode);
}
- else
+ else if (MPFR_IS_ZERO (c))
{
- MPFR_ASSERTD (MPFR_IS_ZERO (c));
return mpfr_set (a, b, rnd_mode);
}
+ else
+ {
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
+ MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
+ /* mpfr_sub1sp and mpfr_add1sp are not intended to support UBF,
+ for which optimization is less important. */
+ if (MPFR_SIGN(b) == MPFR_SIGN(c))
+ return mpfr_sub1 (a, b, c, rnd_mode);
+ else if (MPFR_EXP_LESS_P (b, c))
+ {
+ int inexact;
+ rnd_mode = MPFR_INVERT_RND (rnd_mode);
+ inexact = mpfr_add1 (a, c, b, rnd_mode);
+ MPFR_CHANGE_SIGN (a);
+ return -inexact;
+ }
+ else
+ return mpfr_add1 (a, b, c, rnd_mode);
+ }
}
MPFR_ASSERTD (MPFR_IS_PURE_FP (b));
diff --git a/src/sub1.c b/src/sub1.c
index f062420ab..7d1938af4 100644
--- a/src/sub1.c
+++ b/src/sub1.c
@@ -28,11 +28,13 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
a positive value otherwise.
*/
+/* TODO: check the code in case exp_b == MPFR_EXP_MAX. */
+
int
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
int sign;
- mpfr_exp_t diff_exp;
+ mpfr_exp_t diff_exp, exp_b;
mpfr_prec_t cancel, cancel1;
mp_size_t cancel2, an, bn, cn, cn0;
mp_limb_t *ap, *bp, *cp;
@@ -87,7 +89,17 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
else
MPFR_SET_SAME_SIGN (a,b);
- diff_exp = MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
+ if (MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)))
+ {
+ exp_b = MPFR_IS_UBF (b) ?
+ mpfr_ubf_zexp2exp (MPFR_ZEXP (b)) : MPFR_GET_EXP (b);
+ diff_exp = mpfr_ubf_diff_exp (b, c);
+ }
+ else
+ {
+ exp_b = MPFR_GET_EXP (b);
+ diff_exp = exp_b - MPFR_GET_EXP (c);
+ }
MPFR_ASSERTD (diff_exp >= 0);
aq = MPFR_GET_PREC (a);
@@ -99,62 +111,74 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
but it is more expensive and not very useful */
if (MPFR_UNLIKELY (MAX (aq, bq) + 2 <= diff_exp))
{
+ MPFR_LOG_MSG (("case c small\n", 0));
+
/* Remember, we can't have an exact result! */
/* A.AAAAAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
/* A = S*ABS(B) +/- ulp(a) */
- MPFR_SET_EXP (a, MPFR_GET_EXP (b));
+ MPFR_EXP (a) = exp_b; /* may be up to MPFR_EXP_MAX */
MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq,
- rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a));
+ rnd_mode, MPFR_SIGN (a),
+ if (MPFR_EXP (a) != MPFR_EXP_MAX)
+ ++ MPFR_EXP (a));
+ MPFR_LOG_MSG (("inexact=%d\n", inexact));
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)))
+ /* a = b, but the exact value of b - c is a bit below. Then,
+ except for directed rounding similar to toward zero and
+ before overflow checking: a is the correctly rounded value
+ and since |b| - |c| < |a|, the ternary value value is given
+ by the sign of a. */
+ if (! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
{
- mpfr_nexttozero (a);
- MPFR_RET (- MPFR_INT_SIGN (a));
+ inexact = MPFR_INT_SIGN (a);
+ goto check_overflow;
}
- MPFR_RET (MPFR_INT_SIGN (a));
}
- else
+ else /* inexact != 0 */
{
/* A.AAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
- /* It isn't exact so Prec(b) > Prec(a) and the last
- Prec(b)-Prec(a) bits of `b' are not zeros.
- Which means that removing c from b can't generate a carry
- except in case of even rounding.
- In all other cases the result and the inexact flag should be
- correct (We can't have an exact result).
- In case of EVEN rounding:
+ /* It isn't exact, so PREC(b) > PREC(a) and the last
+ PREC(b)-PREC(a) bits of b are not all zeros.
+ Subtracting c from b will not have an effect on the rounding
+ except in case of a midpoint in the round-to-nearest mode,
+ when the even rounding was done away from zero instead of
+ toward zero.
+ In case of even rounding:
1.BBBBBBBBBBBBBx10
- 1.CCCCCCCCCCCC
- = 1.BBBBBBBBBBBBBx01 Rounded to Prec(b)
- = 1.BBBBBBBBBBBBBx Nearest / Rounded to Prec(a)
+ = 1.BBBBBBBBBBBBBx01 Rounded to PREC(b)
+ = 1.BBBBBBBBBBBBBx Nearest / Rounded to PREC(a)
Set gives:
1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0)
1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
- which means we get a wrong rounded result if x==1,
- i.e. inexact= MPFR_EVEN_INEX */
- if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX * MPFR_INT_SIGN (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);
+ which means we get a wrong rounded result if x == 1,
+ i.e. inexact == MPFR_EVEN_INEX (for positive numbers). */
+ if (MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
+ goto check_overflow;
}
+ /* We need to take the value preceding |a|. We can't use
+ mpfr_nexttozero due to a possible out-of-range exponent.
+ But this will allow us to have more specific code. */
+ MPFR_LOG_MSG (("correcting the value of a\n", 0));
+ sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq;
+ mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
+ if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0))
+ {
+ MPFR_EXP (a) --;
+ /* The following is valid whether an = 1 or an > 1. */
+ ap[an-1] |= MPFR_LIMB_HIGHBIT;
+ }
+ inexact = - MPFR_INT_SIGN (a);
+ check_overflow:
+ if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
+ return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
+ else
+ MPFR_RET (inexact);
}
/* reserve a space to store b aligned with the result, i.e. shifted by
@@ -622,16 +646,21 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
mpfr_exp_t exp_a;
cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
- exp_a = MPFR_GET_EXP (b) - cancel;
- if (MPFR_UNLIKELY(exp_a < __gmpfr_emin))
+ exp_a = exp_b - cancel;
+ if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
{
- MPFR_TMP_FREE(marker);
+ MPFR_TMP_FREE (marker);
if (rnd_mode == MPFR_RNDN &&
(exp_a < __gmpfr_emin - 1 ||
(inexact >= 0 && mpfr_powerof2_raw (a))))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
}
+ if (MPFR_UNLIKELY (exp_a > __gmpfr_emax))
+ {
+ MPFR_TMP_FREE (marker);
+ return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
+ }
MPFR_SET_EXP (a, exp_a);
}
else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
@@ -639,13 +668,11 @@ mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
/* in case cancel = 0, add_exp can still be 1, in case b is just
below a power of two, c is very small, prec(a) < prec(b),
and rnd=away or nearest */
- mpfr_exp_t exp_b;
-
- exp_b = MPFR_GET_EXP (b);
- if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax))
+ MPFR_ASSERTD (add_exp == 0 || add_exp == 1);
+ if (MPFR_UNLIKELY (add_exp && exp_b >= __gmpfr_emax))
{
- MPFR_TMP_FREE(marker);
- return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
+ MPFR_TMP_FREE (marker);
+ return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
}
MPFR_SET_EXP (a, exp_b + add_exp);
}
diff --git a/src/ubf.c b/src/ubf.c
new file mode 100644
index 000000000..e43e0e41d
--- /dev/null
+++ b/src/ubf.c
@@ -0,0 +1,226 @@
+/* Functions to work with unbounded floats (limited low-level interface).
+
+Copyright 2016 Free Software Foundation, Inc.
+Contributed by the AriC and Caramba projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+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. */
+
+#include "mpfr-impl.h"
+
+/* Note: In MPFR math functions, even if UBF code is not called first,
+ we may still need to handle special values NaN and infinities here.
+ Indeed, for MAXR * MAXR + (-inf), even though this is a special case,
+ the computation will also generate an overflow due to MAXR * MAXR,
+ so that UBF code will be called anyway (except if special cases are
+ detected and handled separately, but for polynomial, this should not
+ be needed). */
+
+/* This function does not change the flags. */
+static void
+mpfr_get_zexp (mpz_ptr ez, mpfr_srcptr x)
+{
+ mpz_init (ez);
+
+ if (MPFR_IS_UBF (x))
+ mpz_set (ez, MPFR_ZEXP (x));
+ else
+ {
+ mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
+ mpfr_t e;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ /* TODO: Once this has been tested, optimize based on whether
+ _MPFR_EXP_FORMAT <= 3. */
+ MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT);
+ MPFR_SAVE_EXPO_MARK (expo);
+ MPFR_DBGRES (inex = mpfr_set_exp_t (e, MPFR_GET_EXP (x), MPFR_RNDN));
+ MPFR_ASSERTD (inex == 0);
+ MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN));
+ MPFR_ASSERTD (inex == 0);
+ MPFR_SAVE_EXPO_FREE (expo);
+ }
+}
+
+/* Exact product. The number a is assumed to have enough allocated memory,
+ where the trailing bits are regarded as being part of the input numbers
+ (no reallocation is attempted and no check is performed as MPFR_TMP_INIT
+ could have been used). The arguments b and c may actually be UBF numbers
+ (mpfr_srcptr can be seen a bit like void *, but is stronger).
+ This function does not change the flags, except in case of NaN. */
+void
+mpfr_ubf_mul_exact (mpfr_ubf_ptr a, mpfr_srcptr b, mpfr_srcptr c)
+{
+ MPFR_LOG_FUNC
+ (("b[%Pu]=%.*Rg c[%Pu]=%.*Rg",
+ mpfr_get_prec (b), mpfr_log_prec, b,
+ mpfr_get_prec (c), mpfr_log_prec, c),
+ ("a[%Pu]=%.*Rg",
+ mpfr_get_prec (a), mpfr_log_prec, a));
+
+ MPFR_ASSERTD ((mpfr_ptr) a != b);
+ MPFR_ASSERTD ((mpfr_ptr) a != c);
+ MPFR_SIGN (a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c));
+
+ if (MPFR_ARE_SINGULAR (b, c))
+ {
+ if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
+ MPFR_SET_NAN (a);
+ else if (MPFR_IS_INF (b))
+ {
+ if (MPFR_NOTZERO (c))
+ MPFR_SET_INF (a);
+ else
+ MPFR_SET_NAN (a);
+ }
+ else if (MPFR_IS_INF (c))
+ {
+ if (!MPFR_IS_ZERO (b))
+ MPFR_SET_INF (a);
+ else
+ MPFR_SET_NAN (a);
+ }
+ else
+ {
+ MPFR_ASSERTD (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));
+ MPFR_SET_ZERO (a);
+ }
+ }
+ else
+ {
+ mpfr_exp_t e;
+ mp_size_t bn, cn;
+ mpfr_limb_ptr ap;
+ mp_limb_t u, v;
+ int m;
+
+ /* Note about the code below: For the choice of the precision of
+ * the result a, one could choose PREC(b) + PREC(c), instead of
+ * taking whole limbs into account, but in most cases where one
+ * would gain one limb, one would need to copy the significand
+ * instead of a no-op (see the mul.c code).
+ * But in the case MPFR_LIMB_MSB (u) == 0, if the result fits in
+ * an-1 limbs, one could actually do
+ * mpn_rshift (ap, ap, k, GMP_NUMB_BITS - 1)
+ * instead of
+ * mpn_lshift (ap, ap, k, 1)
+ * to gain one limb (and reduce the precision), replacing a shift
+ * by another one. Would this be interesting?
+ */
+
+ bn = MPFR_LIMB_SIZE (b);
+ cn = MPFR_LIMB_SIZE (c);
+
+ ap = MPFR_MANT (a);
+
+ u = (bn >= cn) ?
+ mpn_mul (ap, MPFR_MANT (b), bn, MPFR_MANT (c), cn) :
+ mpn_mul (ap, MPFR_MANT (c), cn, MPFR_MANT (b), bn);
+ if (MPFR_UNLIKELY (MPFR_LIMB_MSB (u) == 0))
+ {
+ m = 1;
+ MPFR_DBGRES (v = mpn_lshift (ap, ap, bn + cn, 1));
+ MPFR_ASSERTD (v == 0);
+ }
+ else
+ m = 0;
+
+ if (! MPFR_IS_UBF (b) && ! MPFR_IS_UBF (c) &&
+ (e = MPFR_GET_EXP (b) + MPFR_GET_EXP (c) - m,
+ MPFR_EXP_IN_RANGE (e)))
+ {
+ MPFR_SET_EXP (a, e);
+ }
+ else
+ {
+ mpz_t be, ce;
+
+ mpz_init (MPFR_ZEXP (a));
+
+ /* This may involve copies of mpz_t, but exponents should not be
+ very large integers anyway. */
+ mpfr_get_zexp (be, b);
+ mpfr_get_zexp (ce, c);
+ mpz_add (MPFR_ZEXP (a), be, ce);
+ mpz_clear (be);
+ mpz_clear (ce);
+ mpz_sub_ui (MPFR_ZEXP (a), MPFR_ZEXP (a), m);
+ MPFR_SET_UBF (a);
+ }
+ }
+}
+
+int
+mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y)
+{
+ mpz_t xe, ye;
+ int c;
+
+ mpfr_get_zexp (xe, x);
+ mpfr_get_zexp (ye, y);
+ c = mpz_cmp (xe, ye) < 0;
+ mpz_clear (xe);
+ mpz_clear (ye);
+ return c;
+}
+
+/* Convert an mpz_t to an mpfr_exp_t, restricted to
+ the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
+mpfr_exp_t
+mpfr_ubf_zexp2exp (mpz_ptr ez)
+{
+ mp_size_t n;
+ mpfr_eexp_t e;
+ mpfr_t d;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ n = ABSIZ (ez); /* limb size of ez */
+ if (n == 0)
+ return 0;
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (d, n * GMP_NUMB_BITS);
+ MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN));
+ MPFR_ASSERTD (inex == 0);
+ e = mpfr_get_exp_t (d, MPFR_RNDZ);
+ mpfr_clear (d);
+ MPFR_SAVE_EXPO_FREE (expo);
+ if (MPFR_UNLIKELY (e < MPFR_EXP_MIN))
+ return MPFR_EXP_MIN;
+ if (MPFR_UNLIKELY (e > MPFR_EXP_MAX))
+ return MPFR_EXP_MAX;
+ return e;
+}
+
+/* Return the difference of the exponents of x and y, restricted to
+ the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
+mpfr_exp_t
+mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y)
+{
+ mpz_t xe, ye;
+ mpfr_exp_t e;
+
+ mpfr_get_zexp (xe, x);
+ mpfr_get_zexp (ye, y);
+ mpz_sub (xe, xe, ye);
+ mpz_clear (ye);
+ e = mpfr_ubf_zexp2exp (xe);
+ mpz_clear (xe);
+ return e;
+}
diff --git a/src/vasprintf.c b/src/vasprintf.c
index 634b5b189..3e5fb0672 100644
--- a/src/vasprintf.c
+++ b/src/vasprintf.c
@@ -1552,6 +1552,24 @@ partition_number (struct number_parts *np, mpfr_srcptr p,
}
}
}
+ else if (MPFR_UNLIKELY (MPFR_IS_UBF (p)))
+ {
+ /* mpfr_get_str does not support UBF, so that UBF numbers are regarded
+ as special cases here. This is not much a problem since UBF numbers
+ are internal to MPFR and here, they only for logging. */
+ if (np->pad_type == LEADING_ZEROS)
+ /* change to right justification padding with left spaces */
+ np->pad_type = LEFT;
+
+ if (MPFR_IS_NEG (p))
+ np->sign = '-';
+
+ np->ip_size = 3;
+ str = (char *) (*__gmp_allocate_func) (1 + np->ip_size);
+ strcpy (str, uppercase ? "UBF" : "ubf");
+ np->ip_ptr = register_string (np->sl, str);
+ /* TODO: output more information (e.g. the exponent) if need be. */
+ }
else
{
MPFR_ASSERTD (MPFR_IS_PURE_FP (p));
diff --git a/tests/tacosh.c b/tests/tacosh.c
index d10dbe4f2..fbabc3098 100644
--- a/tests/tacosh.c
+++ b/tests/tacosh.c
@@ -176,17 +176,7 @@ huge (void)
{
mpfr_t x, y, z;
int inex;
- mpfr_exp_t emax;
- if (mpfr_get_emax_max () < 1073741823)
- return;
-
- emax = mpfr_get_emax ();
- mpfr_set_emax (1073741823);
-
- /* FIXME: The main purpose of this test was on 32-bit ABI,
- but it is no longer run there. Solution: implement the
- TODO below. */
/* TODO: extend the exponent range and use mpfr_get_emax (). */
mpfr_inits2 (32, x, y, z, (mpfr_ptr) 0);
mpfr_set_ui_2exp (x, 1, 1073741822, MPFR_RNDN);
@@ -203,7 +193,6 @@ huge (void)
MPFR_ASSERTN (inex < 0);
mpfr_clears (x, y, z, (mpfr_ptr) 0);
- mpfr_set_emax (emax);
}
int
diff --git a/tests/texp.c b/tests/texp.c
index 0e85589b4..e0c9e26f2 100644
--- a/tests/texp.c
+++ b/tests/texp.c
@@ -273,26 +273,21 @@ check_special (void)
}
/* Check overflow. Corner case of mpfr_exp_2 */
- /* FIXME: The main purpose of the test below was on 32-bit ABI,
- but it is no longer run there. */
mpfr_set_prec (x, 64);
- if (mpfr_set_emax (1073741823) == 0)
- { /* 1073741823 is in the allowed exponent range */
- mpfr_set_emin (MPFR_EMIN_DEFAULT);
- mpfr_set_str (x,
- "0.1011000101110010000101111111010100001100000001110001100111001101E30",
- 2, MPFR_RNDN);
- mpfr_exp (x, x, MPFR_RNDD);
- if (mpfr_cmp_str (x,
- ".1111111111111111111111111111111111111111111111111111111111111111E1073741823",
- 2, MPFR_RNDN) != 0)
- {
- printf ("Wrong overflow detection in mpfr_exp\n");
- mpfr_dump (x);
- exit (1);
- }
+ mpfr_set_emax (MPFR_EMAX_DEFAULT);
+ mpfr_set_emin (MPFR_EMIN_DEFAULT);
+ mpfr_set_str (x,
+ "0.1011000101110010000101111111010100001100000001110001100111001101E30",
+ 2, MPFR_RNDN);
+ mpfr_exp (x, x, MPFR_RNDD);
+ if (mpfr_cmp_str (x,
+".1111111111111111111111111111111111111111111111111111111111111111E1073741823",
+ 2, MPFR_RNDN) != 0)
+ {
+ printf ("Wrong overflow detection in mpfr_exp\n");
+ mpfr_dump (x);
+ exit (1);
}
-
/* Check underflow. Corner case of mpfr_exp_2 */
mpfr_set_str (x,
"-0.1011000101110010000101111111011111010001110011110111100110101100E30",
diff --git a/tests/tfmma.c b/tests/tfmma.c
index b825de646..c45995a57 100644
--- a/tests/tfmma.c
+++ b/tests/tfmma.c
@@ -22,8 +22,6 @@ http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
#include "mpfr-test.h"
-/* TODO: add more tests, with special values and exception checking. */
-
/* check both mpfr_fmma and mpfr_fmms */
static void
random_test (mpfr_t a, mpfr_t b, mpfr_t c, mpfr_t d, mpfr_rnd_t rnd)
@@ -123,7 +121,7 @@ zero_tests (void)
set_emin (MPFR_EMIN_MIN);
set_emax (MPFR_EMAX_MAX);
- mpfr_inits2 (64, a, b, c, d, res, (mpfr_ptr) 0);
+ mpfr_inits2 (GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0);
for (i = 0; i <= 4; i++)
{
switch (i)
@@ -143,38 +141,278 @@ zero_tests (void)
}
RND_LOOP (r)
{
- int inex;
+ int j, inex;
mpfr_flags_t flags;
mpfr_set (b, a, MPFR_RNDN);
mpfr_set (c, a, MPFR_RNDN);
mpfr_neg (d, a, MPFR_RNDN);
+ /* We also want to test cases where the precision of the
+ result is twice the precision of each input, in case
+ add1sp.c and/or sub1sp.c could be involved. */
+ for (j = 1; j <= 2; j++)
+ {
+ mpfr_init2 (res, GMP_NUMB_BITS * j);
+ mpfr_clear_flags ();
+ inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r);
+ flags = __gmpfr_flags;
+ if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) ||
+ (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res)))
+ {
+ printf ("Error in zero_tests on i = %d, %s\n",
+ i, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
+ printf ("Expected %c0, inex = 0\n",
+ r == MPFR_RNDD ? '-' : '+');
+ printf ("Got ");
+ if (MPFR_IS_POS (res))
+ printf ("+");
+ mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN);
+ printf (", inex = %d\n", inex);
+ printf ("Expected flags:");
+ flags_out (0);
+ printf ("Got flags: ");
+ flags_out (flags);
+ exit (1);
+ }
+ mpfr_clear (res);
+ } /* j */
+ } /* r */
+ } /* i */
+ mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
+
+ set_emin (emin);
+ set_emax (emax);
+}
+
+static void
+max_tests (void)
+{
+ mpfr_exp_t emin, emax;
+ mpfr_t x, y1, y2;
+ int r;
+ int i, inex1, inex2;
+ mpfr_flags_t flags1, flags2;
+
+ emin = mpfr_get_emin ();
+ emax = mpfr_get_emax ();
+ set_emin (MPFR_EMIN_MIN);
+ set_emax (MPFR_EMAX_MAX);
+
+ mpfr_init2 (x, GMP_NUMB_BITS);
+ mpfr_setmax (x, MPFR_EMAX_MAX);
+ flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
+ RND_LOOP (r)
+ {
+ /* We also want to test cases where the precision of the
+ result is twice the precision of each input, in case
+ add1sp.c and/or sub1sp.c could be involved. */
+ for (i = 1; i <= 2; i++)
+ {
+ mpfr_inits2 (GMP_NUMB_BITS * i, y1, y2, (mpfr_ptr) 0);
+ inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r);
mpfr_clear_flags ();
- inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r);
- flags = __gmpfr_flags;
- if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) ||
- (r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res)))
+ inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r);
+ flags2 = __gmpfr_flags;
+ if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+ mpfr_equal_p (y1, y2)))
{
- printf ("Error in zero_tests on i = %d, %s\n",
- i, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
- printf ("Expected %c0, inex = 0\n", r == MPFR_RNDD ? '-' : '+');
+ printf ("Error in max_tests for %s\n",
+ mpfr_print_rnd_mode ((mpfr_rnd_t) r));
+ printf ("Expected ");
+ mpfr_dump (y1);
+ printf (" with inex = %d, flags =", inex1);
+ flags_out (flags1);
printf ("Got ");
- if (MPFR_IS_POS (res))
- printf ("+");
- mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN);
- printf (", inex = %d\n", inex);
- printf ("Expected flags:");
- flags_out (0);
- printf ("Got flags: ");
- flags_out (flags);
+ mpfr_dump (y2);
+ printf (" with inex = %d, flags =", inex2);
+ flags_out (flags2);
exit (1);
}
+ mpfr_clears (y1, y2, (mpfr_ptr) 0);
+ } /* i */
+ } /* r */
+ mpfr_clear (x);
+
+ set_emin (emin);
+ set_emax (emax);
+}
+
+/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't.
+ * a^2 + cd where a^2 overflows and cd doesn't affect the overflow
+ * (and cd may even underflow).
+ */
+static void
+overflow_tests (void)
+{
+ mpfr_exp_t old_emax;
+ int i;
+
+ old_emax = mpfr_get_emax ();
+
+ for (i = 0; i < 40; i++)
+ {
+ mpfr_exp_t emax, exp_a;
+ mpfr_t a, k, c, d, z1, z2;
+ mpfr_rnd_t rnd;
+ mpfr_prec_t prec_a, prec_z;
+ int add = i & 1, inex, inex1, inex2;
+ mpfr_flags_t flags1, flags2;
+
+ /* In most cases, we do the test with the maximum exponent. */
+ emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + (randlimb () % 1);
+ set_emax (emax);
+ exp_a = emax/2 + 32;
+
+ rnd = RND_RAND ();
+ prec_a = 64 + randlimb () % 100;
+ prec_z = MPFR_PREC_MIN + randlimb () % 160;
+
+ mpfr_init2 (a, prec_a);
+ mpfr_urandom (a, RANDS, MPFR_RNDU);
+ mpfr_set_exp (a, exp_a);
+
+ mpfr_init2 (k, prec_a - 32);
+ mpfr_urandom (k, RANDS, MPFR_RNDU);
+ mpfr_set_exp (k, exp_a - 32);
+
+ mpfr_init2 (c, prec_a + 1);
+ inex = mpfr_add (c, a, k, MPFR_RNDN);
+ MPFR_ASSERTN (inex == 0);
+
+ mpfr_init2 (d, prec_a);
+ inex = mpfr_sub (d, a, k, MPFR_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ if (add)
+ mpfr_neg (d, d, MPFR_RNDN);
+
+ mpfr_init2 (z1, prec_z);
+ mpfr_clear_flags ();
+ inex1 = mpfr_sqr (z1, k, rnd);
+ flags1 = __gmpfr_flags;
+
+ mpfr_init2 (z2, prec_z);
+ mpfr_clear_flags ();
+ inex2 = add ?
+ mpfr_fmma (z2, a, a, c, d, rnd) :
+ mpfr_fmms (z2, a, a, c, d, rnd);
+ flags2 = __gmpfr_flags;
+
+ if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+ mpfr_equal_p (z1, z2)))
+ {
+ printf ("Error 1 in overflow_tests for %s\n",
+ mpfr_print_rnd_mode (rnd));
+ printf ("Expected ");
+ mpfr_dump (z1);
+ printf (" with inex = %d, flags =", inex1);
+ flags_out (flags1);
+ printf ("Got ");
+ mpfr_dump (z2);
+ printf (" with inex = %d, flags =", inex2);
+ flags_out (flags2);
+ exit (1);
+ }
+
+ /* c and d such that a^2 +/- cd ~= a^2 (overflow) */
+ mpfr_urandom (c, RANDS, MPFR_RNDU);
+ mpfr_set_exp (c, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin);
+ if (randlimb () % 1)
+ mpfr_neg (c, c, MPFR_RNDN);
+ mpfr_urandom (d, RANDS, MPFR_RNDU);
+ mpfr_set_exp (d, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin);
+ if (randlimb () % 1)
+ mpfr_neg (d, d, MPFR_RNDN);
+
+ mpfr_clear_flags ();
+ inex1 = mpfr_sqr (z1, a, rnd);
+ flags1 = __gmpfr_flags;
+ MPFR_ASSERTN (flags1 == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT));
+
+ mpfr_clear_flags ();
+ inex2 = add ?
+ mpfr_fmma (z2, a, a, c, d, rnd) :
+ mpfr_fmms (z2, a, a, c, d, rnd);
+ flags2 = __gmpfr_flags;
+
+ if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
+ mpfr_equal_p (z1, z2)))
+ {
+ printf ("Error 2 in overflow_tests for %s\n",
+ mpfr_print_rnd_mode (rnd));
+ printf ("Expected ");
+ mpfr_dump (z1);
+ printf (" with inex = %d, flags =", inex1);
+ flags_out (flags1);
+ printf ("Got ");
+ mpfr_dump (z2);
+ printf (" with inex = %d, flags =", inex2);
+ flags_out (flags2);
+ exit (1);
+ }
+
+ mpfr_clears (a, k, c, d, z1, z2, (mpfr_ptr) 0);
+ }
+
+ set_emax (old_emax);
+}
+
+/* (1/2)x + (1/2)x = x tested near underflow */
+static void
+half_plus_half (void)
+{
+ mpfr_exp_t emin;
+ mpfr_t h, x1, x2, y;
+ int neg, r, i, inex;
+ mpfr_flags_t flags;
+
+ emin = mpfr_get_emin ();
+ set_emin (MPFR_EMIN_MIN);
+ mpfr_inits2 (4, h, x1, (mpfr_ptr) 0);
+ mpfr_init2 (x2, GMP_NUMB_BITS);
+ mpfr_set_ui_2exp (h, 1, -1, MPFR_RNDN);
+
+ for (mpfr_setmin (x1, __gmpfr_emin);
+ MPFR_GET_EXP (x1) < __gmpfr_emin + 2;
+ mpfr_nextabove (x1))
+ {
+ inex = mpfr_set (x2, x1, MPFR_RNDN);
+ MPFR_ASSERTN (inex == 0);
+ for (neg = 0; neg < 2; neg++)
+ {
+ RND_LOOP (r)
+ {
+ /* We also want to test cases where the precision of the
+ result is twice the precision of each input, in case
+ add1sp.c and/or sub1sp.c could be involved. */
+ for (i = 1; i <= 2; i++)
+ {
+ mpfr_init2 (y, GMP_NUMB_BITS * i);
+ mpfr_clear_flags ();
+ inex = mpfr_fmma (y, h, x2, h, x2, (mpfr_rnd_t) r);
+ flags = __gmpfr_flags;
+ if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2)))
+ {
+ printf ("Error in half_plus_half for %s\n",
+ mpfr_print_rnd_mode ((mpfr_rnd_t) r));
+ printf ("Expected ");
+ mpfr_dump (x2);
+ printf (" with inex = 0, flags =");
+ flags_out (0);
+ printf ("Got ");
+ mpfr_dump (y);
+ printf (" with inex = %d, flags =", inex);
+ flags_out (flags);
+ exit (1);
+ }
+ mpfr_clear (y);
+ }
+ }
+ mpfr_neg (x2, x2, MPFR_RNDN);
}
}
- mpfr_clears (a, b, c, d, res, (mpfr_ptr) 0);
+ mpfr_clears (h, x1, x2, (mpfr_ptr) 0);
set_emin (emin);
- set_emax (emax);
}
int
@@ -184,6 +422,9 @@ main (int argc, char *argv[])
random_tests ();
zero_tests ();
+ max_tests ();
+ overflow_tests ();
+ half_plus_half ();
tests_end_mpfr ();
return 0;
diff --git a/tests/tpow.c b/tests/tpow.c
index 3bcf7a9a7..13819ff07 100644
--- a/tests/tpow.c
+++ b/tests/tpow.c
@@ -1406,13 +1406,7 @@ bug20071218 (void)
{
mpfr_t x, y, z, t;
int tern;
- mpfr_exp_t emin;
-
- if (mpfr_get_emin_min () > -1073741823)
- return;
- emin = mpfr_get_emin ();
- mpfr_set_emin (-1073741823);
mpfr_inits2 (64, x, y, z, t, (mpfr_ptr) 0);
mpfr_set_str (x, "0x.80000000000002P-1023", 0, MPFR_RNDN);
mpfr_set_str (y, "100000.000000002", 16, MPFR_RNDN);
@@ -1448,7 +1442,6 @@ bug20071218 (void)
exit (1);
}
mpfr_clears (x, y, z, t, (mpfr_ptr) 0);
- mpfr_set_emin (emin);
}
/* With revision 5429, this gives:
diff --git a/tests/tpow_z.c b/tests/tpow_z.c
index 1514bd041..004435473 100644
--- a/tests/tpow_z.c
+++ b/tests/tpow_z.c
@@ -302,28 +302,22 @@ static void
bug20080223 (void)
{
mpfr_t a, exp, answer;
- mpfr_exp_t emin = mpfr_get_emin ();
- mpfr_set_emin (mpfr_get_emin_min ());
- if (mpfr_get_emin () <= -1073741823)
- {
- mpfr_init2 (a, 53);
- mpfr_init2 (exp, 53);
- mpfr_init2 (answer, 53);
+ mpfr_init2 (a, 53);
+ mpfr_init2 (exp, 53);
+ mpfr_init2 (answer, 53);
- mpfr_set_si (exp, -1073741824, MPFR_RNDN);
+ mpfr_set_si (exp, -1073741824, MPFR_RNDN);
- mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN);
- /* a = 562949953139837/2^48 */
- mpfr_pow (answer, a, exp, MPFR_RNDN);
- mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823");
- MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0);
+ mpfr_set_str (a, "1.999999999", 10, MPFR_RNDN);
+ /* a = 562949953139837/2^48 */
+ mpfr_pow (answer, a, exp, MPFR_RNDN);
+ mpfr_set_str_binary (a, "0.110110101111011001110000111111100011101000111011101E-1073741823");
+ MPFR_ASSERTN(mpfr_cmp0 (answer, a) == 0);
- mpfr_clear (a);
- mpfr_clear (exp);
- mpfr_clear (answer);
- }
- mpfr_set_emin (emin);
+ mpfr_clear (a);
+ mpfr_clear (exp);
+ mpfr_clear (answer);
}
static void
diff --git a/tests/tsprintf.c b/tests/tsprintf.c
index fba32511e..b33355c88 100644
--- a/tests/tsprintf.c
+++ b/tests/tsprintf.c
@@ -1128,7 +1128,6 @@ bug20111102 (void)
* for %Ra and %Rb is not done on the MPFR number itself (as it
* would overflow). Note: it has been reported on comp.std.c that
* some C libraries behave differently on %a, but this is a bug.
- * FIXME: this function assumes e = 3 mod 4.
*/
static void
check_emax_aux (mpfr_exp_t e)
@@ -1138,9 +1137,6 @@ check_emax_aux (mpfr_exp_t e)
int i;
mpfr_exp_t emax;
- if ((e % 4) != 3)
- return;
-
MPFR_ASSERTN (e <= LONG_MAX);
emax = mpfr_get_emax ();
set_emax (e);
@@ -1199,7 +1195,6 @@ check_emax (void)
check_emax_aux (MPFR_EMAX_MAX);
}
-/* FIXME: this function assumes e = 1 mod 4 */
static void
check_emin_aux (mpfr_exp_t e)
{
@@ -1209,9 +1204,6 @@ check_emin_aux (mpfr_exp_t e)
mpfr_exp_t emin;
mpz_t ee;
- if ((e % 4) != 1)
- return;
-
MPFR_ASSERTN (e >= LONG_MIN);
emin = mpfr_get_emin ();
set_emin (e);
diff --git a/tests/tsub1sp.c b/tests/tsub1sp.c
index 460a5dd6f..b191e350f 100644
--- a/tests/tsub1sp.c
+++ b/tests/tsub1sp.c
@@ -351,19 +351,12 @@ check_special (void)
MPFR_SET_NAN(x);
MPFR_SET_NAN(x2);
-#if 0
mpfr_set_str_binary (y,
"0.1000000000000000000000000000000000000000000000000000000000000000"
"E-1073741823");
mpfr_set_str_binary (z,
"0.1100000000000000000000000000000000000000000000000000000000000000"
"E-1073741823");
-#else
- mpfr_set_ui_2exp (y, 1, MPFR_EMIN_DEFAULT-1, MPFR_RNDN);
- /* y = 0.5*2^EMIN_DEFAULT */
- mpfr_set_ui_2exp (z, 3, MPFR_EMIN_DEFAULT-2, MPFR_RNDN);
- /* z = 0.75*2^EMIN_DEFAULT */
-#endif
inexact1 = mpfr_sub1(x2, y, z, (mpfr_rnd_t) r);
inexact2 = mpfr_sub1sp(x, y, z, (mpfr_rnd_t) r);
if (mpfr_cmp(x, x2))
View Lorry Controller Status