diff options
author | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2003-10-28 16:31:13 +0000 |
---|---|---|
committer | pelissip <pelissip@280ebfd0-de03-0410-8827-d642c229c3f4> | 2003-10-28 16:31:13 +0000 |
commit | 117edf14c822a22cdd9c25689aeadc904a1a30d1 (patch) | |
tree | e39bd61cefc24cc6cfbc5b2b956e4fb36015d111 | |
parent | 734c0a144b04e2cae4e67b394010e3f6e3cadecc (diff) | |
download | mpfr-117edf14c822a22cdd9c25689aeadc904a1a30d1.tar.gz |
Use of MPFR_UNLIKELY and MPFR_IS_SINGULAR for fast detection of special values (Nan, Inf or Zero).
Start to encapsulate the sign to be independant of the reprensation (Must be 1 or -1).
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2525 280ebfd0-de03-0410-8827-d642c229c3f4
-rw-r--r-- | acos.c | 31 | ||||
-rw-r--r-- | acosh.c | 36 | ||||
-rw-r--r-- | add.c | 95 | ||||
-rw-r--r-- | add1.c | 4 | ||||
-rw-r--r-- | add_one_ulp.c | 15 | ||||
-rw-r--r-- | agm.c | 65 | ||||
-rw-r--r-- | asin.c | 28 | ||||
-rw-r--r-- | asinh.c | 44 | ||||
-rw-r--r-- | atan.c | 77 | ||||
-rw-r--r-- | atanh.c | 66 | ||||
-rw-r--r-- | cbrt.c | 54 | ||||
-rw-r--r-- | cmp.c | 38 | ||||
-rw-r--r-- | cmp2.c | 4 | ||||
-rw-r--r-- | cmp_abs.c | 29 | ||||
-rw-r--r-- | cos.c | 19 | ||||
-rw-r--r-- | cosh.c | 36 | ||||
-rw-r--r-- | div.c | 113 | ||||
-rw-r--r-- | eq.c | 39 | ||||
-rw-r--r-- | erf.c | 27 | ||||
-rw-r--r-- | exp.c | 41 | ||||
-rw-r--r-- | exp2.c | 48 | ||||
-rw-r--r-- | expm1.c | 53 | ||||
-rw-r--r-- | fma.c | 126 | ||||
-rw-r--r-- | frac.c | 6 | ||||
-rw-r--r-- | gamma.c | 56 | ||||
-rw-r--r-- | gammaPiAGMformula.c | 31 | ||||
-rw-r--r-- | hypot.c | 40 | ||||
-rw-r--r-- | log.c | 59 | ||||
-rw-r--r-- | log10.c | 64 | ||||
-rw-r--r-- | log1p.c | 77 | ||||
-rw-r--r-- | log2.c | 68 | ||||
-rw-r--r-- | minmax.c | 73 | ||||
-rw-r--r-- | mpfr-impl.h | 24 | ||||
-rw-r--r-- | mpfr.h | 12 | ||||
-rw-r--r-- | mul.c | 95 | ||||
-rw-r--r-- | mul_ui.c | 63 | ||||
-rw-r--r-- | next.c | 27 | ||||
-rw-r--r-- | pow.c | 192 | ||||
-rw-r--r-- | pow_si.c | 67 | ||||
-rw-r--r-- | pow_ui.c | 71 | ||||
-rw-r--r-- | round_prec.c | 79 | ||||
-rw-r--r-- | set.c | 87 | ||||
-rw-r--r-- | set_si.c | 15 | ||||
-rw-r--r-- | set_ui.c | 8 | ||||
-rw-r--r-- | sin.c | 29 | ||||
-rw-r--r-- | sinh.c | 42 | ||||
-rw-r--r-- | sqrt.c | 80 | ||||
-rw-r--r-- | sub.c | 101 | ||||
-rw-r--r-- | sub_one_ulp.c | 16 | ||||
-rw-r--r-- | tan.c | 26 | ||||
-rw-r--r-- | tanh.c | 53 | ||||
-rw-r--r-- | tests/tpow.c | 14 | ||||
-rw-r--r-- | ui_div.c | 62 | ||||
-rw-r--r-- | ui_sub.c | 35 | ||||
-rw-r--r-- | zeta.c | 45 |
55 files changed, 1427 insertions, 1378 deletions
@@ -40,11 +40,22 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) int compared; int inexact = 0; - /* Trivial cases */ - if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + /* Special cases */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(acos); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + { + MPFR_SET_NAN(acos); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(x)) + { + /* acos(0)=Pi/2 */ + mpfr_const_pi (acos, rnd_mode); + MPFR_SET_EXP (acos, MPFR_GET_EXP (acos) - 1); + return 1; /* inexact */ + } + MPFR_ASSERTN(1); } /* Set x_p=|x| */ @@ -64,7 +75,7 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) if (compared == 0) { mpfr_clear (xp); - if (signe > 0) /* acos(+1) = 0 */ + if (MPFR_IS_POS_SIGN(signe)) /* acos(+1) = 0 */ return mpfr_set_ui (acos, 0, rnd_mode); else /* acos(-1) = Pi */ { @@ -73,18 +84,10 @@ mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) } } - if (MPFR_IS_ZERO(x)) /* acos(0)=Pi/2 */ - { - mpfr_clear (xp); - mpfr_const_pi (acos, rnd_mode); - MPFR_SET_EXP (acos, MPFR_GET_EXP (acos) - 1); - return 1; /* inexact */ - } - prec_acos = MPFR_PREC(acos); mpfr_ui_sub (xp, 1, xp, GMP_RNDD); - if (signe > 0) + if (MPFR_IS_POS_SIGN(signe)) supplement = 2 - 2 * MPFR_GET_EXP (xp); else supplement = 2 - MPFR_GET_EXP(xp); @@ -36,29 +36,35 @@ mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) int inexact = 0; int comp; - if (MPFR_IS_NAN(x) || (comp = mpfr_cmp_ui (x, 1)) < 0) + /* Deal with special cases */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) + { + /* Nan, or zero or -Inf */ + if (MPFR_IS_NAN(x) || MPFR_IS_ZERO(x) || MPFR_IS_NEG(x) ) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(x)) + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + MPFR_RET(0); + } + MPFR_ASSERTN(1); + } + else if (MPFR_UNLIKELY( (comp = mpfr_cmp_ui (x, 1)) < 0)) { MPFR_SET_NAN(y); MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - if (comp == 0) + } + MPFR_CLEAR_FLAGS(y); + if (MPFR_UNLIKELY( comp == 0 )) { MPFR_SET_ZERO(y); /* acosh(1) = 0 */ MPFR_SET_POS(y); MPFR_RET(0); } - - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); /* General case */ { @@ -28,59 +28,60 @@ MA 02111-1307, USA. */ int mpfr_add (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { - if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + if (MPFR_ARE_SINGULAR(b,c)) { - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(a); - - if (MPFR_IS_INF(b)) - { - if (!MPFR_IS_INF(c) || MPFR_SIGN(b) == MPFR_SIGN(c)) - { - MPFR_SET_INF(a); - MPFR_SET_SAME_SIGN(a, b); - MPFR_RET(0); /* exact */ - } + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(a); + if (MPFR_IS_INF(b)) + { + if (!MPFR_IS_INF(c) || MPFR_SIGN(b) == MPFR_SIGN(c)) + { + MPFR_SET_INF(a); + MPFR_SET_SAME_SIGN(a, b); + MPFR_RET(0); /* exact */ + } + else + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + } else - { - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - } - else - if (MPFR_IS_INF(c)) - { - MPFR_SET_INF(a); - MPFR_SET_SAME_SIGN(a, c); - MPFR_RET(0); /* exact */ - } - - MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c)); - - if (MPFR_IS_ZERO(b)) - { + if (MPFR_IS_INF(c)) + { + MPFR_SET_INF(a); + MPFR_SET_SAME_SIGN(a, c); + MPFR_RET(0); /* exact */ + } + if (MPFR_IS_ZERO(b)) + { + if (MPFR_IS_ZERO(c)) + { + MPFR_SET_SIGN(a, + (rnd_mode != GMP_RNDD ? + ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) < 0) ? -1 : 1) : + ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) > 0) ? 1 : -1))); + MPFR_CLEAR_INF(a); + MPFR_SET_ZERO(a); + MPFR_RET(0); /* 0 + 0 is exact */ + } + return mpfr_set(a, c, rnd_mode); + } if (MPFR_IS_ZERO(c)) - { - MPFR_SET_SIGN(a, - (rnd_mode != GMP_RNDD ? - ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) < 0) ? -1 : 1) : - ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) > 0) ? 1 : -1))); - MPFR_CLEAR_INF(a); - MPFR_SET_ZERO(a); - MPFR_RET(0); /* 0 + 0 is exact */ - } - return mpfr_set(a, c, rnd_mode); + { + return mpfr_set(a, b, rnd_mode); + } + /* Should never reach here */ + MPFR_ASSERTN(1); } - if (MPFR_IS_ZERO(c)) - { - return mpfr_set(a, b, rnd_mode); - } + MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); - MPFR_CLEAR_INF(a); /* clear Inf flag */ + MPFR_CLEAR_FLAGS(a); /* clear flags */ if (MPFR_SIGN(b) != MPFR_SIGN(c)) { /* signs differ, it's a subtraction */ @@ -40,8 +40,8 @@ mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, int sh, rb, fb, inex; TMP_DECL(marker); - MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_NOTZERO(b)); - MPFR_ASSERTN(MPFR_IS_FP(c) && MPFR_NOTZERO(c)); + MPFR_ASSERTD(MPFR_IS_PURE_FP(b)); + MPFR_ASSERTD(MPFR_IS_PURE_FP(c)); TMP_MARK(marker); ap = MPFR_MANT(a); diff --git a/add_one_ulp.c b/add_one_ulp.c index 5ae30196b..2ffe731e9 100644 --- a/add_one_ulp.c +++ b/add_one_ulp.c @@ -32,11 +32,14 @@ mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode) int sh; mp_limb_t *xp; - if (MPFR_IS_NAN(x)) - MPFR_RET_NAN; - - if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) - return 0; + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) + { + if (MPFR_IS_NAN(x)) + MPFR_RET_NAN; + if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) + return 0; + MPFR_ASSERTN(1); + } xn = 1 + (MPFR_PREC(x) - 1) / BITS_PER_MP_LIMB; sh = (mp_prec_t) xn * BITS_PER_MP_LIMB - MPFR_PREC(x); @@ -44,7 +47,7 @@ mpfr_add_one_ulp (mpfr_ptr x, mp_rnd_t rnd_mode) if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ { mp_exp_t exp = MPFR_EXP (x); - if (exp == __gmpfr_emax) + if (MPFR_UNLIKELY(exp == __gmpfr_emax)) return mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x)); else { @@ -34,39 +34,40 @@ mpfr_agm (mpfr_ptr r, mpfr_srcptr op2, mpfr_srcptr op1, mp_rnd_t rnd_mode) mpfr_t u, v, tmp, tmpu, tmpv, a, b; TMP_DECL(marker); - /* If a or b is NaN, the result is NaN */ - if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) + /* Deal with special values */ + if (MPFR_ARE_SINGULAR(op1, op2)) { - MPFR_SET_NAN(r); - __gmpfr_flags |= MPFR_FLAGS_NAN; - MPFR_RET_NAN; - } - - /* If a or b is negative (including -Infinity), the result is NaN */ - if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0)) - { - MPFR_SET_NAN(r); - __gmpfr_flags |= MPFR_FLAGS_NAN; - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(r); - - /* If a or b is +Infinity, the result is +Infinity */ - if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2)) - { - MPFR_SET_INF(r); - MPFR_SET_SAME_SIGN(r, op1); - MPFR_RET(0); /* exact */ - } - - MPFR_CLEAR_INF(r); - - /* If a or b is 0, the result is 0 */ - if ((MPFR_NOTZERO(op1) && MPFR_NOTZERO(op2)) == 0) - { - MPFR_SET_ZERO(r); - MPFR_RET(0); /* exact */ + /* If a or b is NaN, the result is NaN */ + if (MPFR_IS_NAN(op1) || MPFR_IS_NAN(op2)) + { + MPFR_SET_NAN(r); + __gmpfr_flags |= MPFR_FLAGS_NAN; + MPFR_RET_NAN; + } + /* If a or b is negative (including -Infinity), the result is NaN */ + if ((MPFR_SIGN(op1) < 0) || (MPFR_SIGN(op2) < 0)) + { + MPFR_SET_NAN(r); + __gmpfr_flags |= MPFR_FLAGS_NAN; + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(r); + /* If a or b is +Infinity, the result is +Infinity */ + if (MPFR_IS_INF(op1) || MPFR_IS_INF(op2)) + { + MPFR_SET_INF(r); + MPFR_SET_SAME_SIGN(r, op1); + MPFR_RET(0); /* exact */ + } + MPFR_CLEAR_INF(r); + /* If a or b is 0, the result is 0 */ + if (MPFR_IS_ZERO(op1) || MPFR_IS_ZERO(op2)) + { + MPFR_SET_POS(r); + MPFR_SET_ZERO(r); + MPFR_RET(0); /* exact */ + } + MPFR_ASSERTN(1); } /* precision of the following calculus */ @@ -41,19 +41,26 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) int compared; int inexact; - /* Trivial cases */ - if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + /* Special cases */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(asin); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + { + MPFR_SET_NAN(asin); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(x)) /* x = 0 */ + { + mpfr_set_ui (asin, 0, GMP_RNDN); + MPFR_RET(0); /* exact result */ + } + MPFR_ASSERTN(1); } /* Set x_p=|x| */ signe = MPFR_SIGN(x); mpfr_init2 (xp, MPFR_PREC(x)); - mpfr_set (xp, x, rnd_mode); - if (signe == -1) - MPFR_CHANGE_SIGN(xp); + mpfr_abs (xp, x, rnd_mode); compared = mpfr_cmp_ui (xp, 1); @@ -82,13 +89,6 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode) return inexact; } - if (MPFR_IS_ZERO(x)) /* x = 0 */ - { - mpfr_set_ui (asin, 0, GMP_RNDN); - mpfr_clear (xp); - return 0; /* exact result */ - } - prec_asin = MPFR_PREC(asin); mpfr_ui_sub (xp, 1, xp, GMP_RNDD); @@ -37,34 +37,34 @@ mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_t t, te, ti; /* auxiliary variables */ long int err; - if (MPFR_IS_NAN(x)) - { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_ZERO(y); /* asinh(0) = 0 */ - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + if (MPFR_IS_INF(x)) + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + MPFR_CLEAR_INF(y); + if (MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); /* asinh(0) = 0 */ + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + MPFR_ASSERTN(1); } Nx = MPFR_PREC(x); /* Precision of input variable */ Ny = MPFR_PREC(y); /* Precision of output variable */ - neg = MPFR_SIGN(x) < 0; + neg = MPFR_IS_NEG(x); /* General case */ @@ -82,45 +82,48 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) int logn; /* Trivial cases */ - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(arctangent); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(x)) - { - MPFR_CLEAR_FLAGS(arctangent); - if (MPFR_SIGN(x) > 0) /* arctan(+inf) = Pi/2 */ - inexact = mpfr_const_pi (arctangent, rnd_mode); - else /* arctan(-inf) = -Pi/2 */ + if (MPFR_IS_NAN(x)) { - if (rnd_mode == GMP_RNDU) - rnd_mode = GMP_RNDD; - else if (rnd_mode == GMP_RNDD) - rnd_mode = GMP_RNDU; - inexact = -mpfr_const_pi (arctangent, rnd_mode); - MPFR_CHANGE_SIGN (arctangent); + MPFR_SET_NAN(arctangent); + MPFR_RET_NAN; } - MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 1); - return inexact; + else if (MPFR_IS_INF(x)) + { + MPFR_CLEAR_FLAGS(arctangent); + if (MPFR_IS_POS(x)) + /* arctan(+inf) = Pi/2 */ + inexact = mpfr_const_pi (arctangent, rnd_mode); + else + /* arctan(-inf) = -Pi/2 */ + { + if (rnd_mode == GMP_RNDU) + rnd_mode = GMP_RNDD; + else if (rnd_mode == GMP_RNDD) + rnd_mode = GMP_RNDU; + inexact = -mpfr_const_pi (arctangent, rnd_mode); + MPFR_CHANGE_SIGN (arctangent); + } + MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 1); + return inexact; + } + else if (MPFR_IS_ZERO(x)) + { + mpfr_set_ui (arctangent, 0, GMP_RNDN); + return 0; /* exact result */ + } + MPFR_ASSERTN(1); } MPFR_CLEAR_FLAGS(arctangent); - if (MPFR_IS_ZERO(x)) - { - mpfr_set_ui (arctangent, 0, GMP_RNDN); - return 0; /* exact result */ - } signe = MPFR_SIGN(x); prec_arctan = MPFR_PREC(arctangent); /* Set x_p=|x| */ mpfr_init2 (xp, MPFR_PREC(x)); - mpfr_set (xp, x, rnd_mode); - if (signe == -1) - MPFR_CHANGE_SIGN(xp); + mpfr_abs (xp, x, rnd_mode); /* Other simple case arctang(-+1)=-+pi/4 */ comparaison = mpfr_cmp_ui (xp, 1); @@ -128,7 +131,7 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) { inexact = mpfr_const_pi (arctangent, rnd_mode); MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 2); - if (signe == -1) + if (MPFR_IS_NEG_SIGN( signe )) { inexact = -inexact; MPFR_CHANGE_SIGN(arctangent); @@ -144,7 +147,8 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) prec_x = __gmpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB); logn = __gmpfr_ceil_log2 ((double) prec_x); - if (logn < 2) logn = 2; + if (logn < 2) + logn = 2; realprec = prec_arctan + __gmpfr_ceil_log2((double) prec_arctan) + 4; mpz_init (ukz); mpz_init (square); @@ -191,7 +195,7 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) /* Calculation of arctan(Ak) */ mpz_mul(square, ukz, ukz); - mpz_neg(square, square); + mpz_neg(square, square); mpfr_atan_aux(t_arctan, square, 2*twopoweri, N0 - i); mpfr_set_z(Ak, ukz, GMP_RNDN); mpfr_div_2ui(Ak, Ak, twopoweri, GMP_RNDN); @@ -210,17 +214,10 @@ mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode) } if (comparaison > 0) - { - mpfr_sub(arctgt, Pisur2, tmp_arctan, GMP_RNDN); - if (signe == -1) - MPFR_CHANGE_SIGN(arctgt); - } + mpfr_sub(arctgt, Pisur2, tmp_arctan, GMP_RNDN); else - { - mpfr_set(arctgt, tmp_arctan, GMP_RNDN); - if (signe == -1) - MPFR_CHANGE_SIGN(arctgt); - } + mpfr_set(arctgt, tmp_arctan, GMP_RNDN); + MPFR_SET_POS(arctgt); if (mpfr_can_round (arctgt, realprec, GMP_RNDN, GMP_RNDZ, MPFR_PREC (arctangent) + (rnd_mode == GMP_RNDN))) @@ -34,49 +34,36 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) { int inexact =0; mpfr_t x; - int flag_neg=0; - mp_prec_t Nx=MPFR_PREC(xt); /* Precision of input variable */ - mpfr_init2(x,Nx); - mpfr_set(x,xt,GMP_RNDN); - - if (MPFR_SIGN(x) < 0) - { - MPFR_CHANGE_SIGN(x); - flag_neg=1; - } - if (MPFR_IS_NAN(x)) - { - MPFR_SET_NAN(y); - mpfr_clear(x); - MPFR_RET_NAN; - } - MPFR_CLEAR_NAN(y); - - - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y,x); - if(flag_neg) - MPFR_CHANGE_SIGN(y); - mpfr_clear(x); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(x)) + /* Special cases */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(xt) )) { - MPFR_SET_ZERO(y); /* atanh(0) = 0 */ - MPFR_SET_SAME_SIGN(y,x); - if(flag_neg) - MPFR_CHANGE_SIGN(y); - mpfr_clear(x); - MPFR_RET(0); + if (MPFR_IS_NAN(xt)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + if (MPFR_IS_INF(xt)) + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, xt); + MPFR_RET(0); + } + MPFR_CLEAR_INF(y); + if (MPFR_IS_ZERO(xt)) + { + MPFR_SET_ZERO(y); /* atanh(0) = 0 */ + MPFR_SET_SAME_SIGN(y,xt); + MPFR_RET(0); + } + MPFR_ASSERTN(1); } + mpfr_init2(x,Nx); + mpfr_abs(x, xt, GMP_RNDN); + /* General case */ { /* Declaration of the intermediary variable */ @@ -126,7 +113,7 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) Ny + (rnd_mode == GMP_RNDN)) || MPFR_IS_ZERO(t))); - if(flag_neg) + if (MPFR_IS_NEG(xt)) MPFR_CHANGE_SIGN(t); inexact = mpfr_set (y, t, rnd_mode); @@ -138,3 +125,4 @@ mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mpfr_clear(x); return inexact; } + @@ -52,40 +52,38 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpz_t m; mp_exp_t e, r, sh; mp_prec_t n, size_m; - int inexact, sign_x; + int inexact, x_neg; /* special values */ - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + if (MPFR_IS_INF(x)) + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN (y, x); + return 0; + } + MPFR_CLEAR_INF(y); + /* case 0: cbrt(+/- 0) = +/- 0 */ + if (MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN (y, x); + return 0; + } + MPFR_ASSERTN(1); } - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN (y, x); - return 0; - } - - MPFR_CLEAR_INF(y); - - /* case 0: cbrt(+/- 0) = +/- 0 */ - if (MPFR_IS_ZERO(x)) - { - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN (y, x); - return 0; - } - - sign_x = MPFR_SIGN(x); - mpz_init (m); e = mpfr_get_z_exp (m, x); /* x = m * 2^e */ - if (sign_x < 0) + if ((x_neg = MPFR_IS_NEG(x))) mpz_neg (m, m); r = e % 3; if (r < 0) @@ -136,9 +134,9 @@ mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) inexact += mpfr_set_z (y, m, GMP_RNDN); MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3); - if (sign_x < 0) + if (x_neg) { - mpfr_neg (y, y, GMP_RNDN); + MPFR_CHANGE_SIGN(y); inexact = -inexact; } @@ -36,25 +36,29 @@ mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, int s) mp_size_t bn, cn; mp_limb_t *bp, *cp; - MPFR_ASSERTN(!MPFR_IS_NAN(b)); - MPFR_ASSERTN(!MPFR_IS_NAN(c)); - s *= MPFR_SIGN(c); - - if (MPFR_IS_INF(b)) - return (MPFR_IS_INF(c) && s * MPFR_SIGN(b) > 0) ? 0 : MPFR_SIGN(b); - - if (MPFR_IS_INF(c)) - return -s; + s = MPFR_MULT_SIGN( s , MPFR_SIGN(c) ); + if (MPFR_ARE_SINGULAR(b,c)) + { + if (MPFR_IS_INF(b)) + { + if (MPFR_IS_INF(c) && s == MPFR_SIGN(b) ) + return 0; + else + return MPFR_SIGN(b); + } + else if (MPFR_IS_INF(c)) + return -s; + else if (MPFR_IS_ZERO(b)) + return MPFR_IS_ZERO(c) ? 0 : -s; + else if (MPFR_IS_ZERO(c)) + return MPFR_SIGN(b); + else + /* Should never reach here */ + MPFR_ASSERTN(1); + } /* b and c are real numbers */ - - if (MPFR_IS_ZERO(b)) - return MPFR_IS_ZERO(c) ? 0 : -s; - - if (MPFR_IS_ZERO(c)) - return MPFR_SIGN(b); - - if (s * MPFR_SIGN(b) < 0) + if (s != MPFR_SIGN(b)) return MPFR_SIGN(b); /* now signs are equal */ @@ -43,8 +43,8 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel) mp_prec_t res = 0; int sign; - MPFR_ASSERTN(MPFR_IS_FP(b)); - MPFR_ASSERTN(MPFR_IS_FP(c)); + MPFR_ASSERTD(MPFR_IS_FP(b)); + MPFR_ASSERTD(MPFR_IS_FP(c)); /* Optimized case x - x */ if (b == c) @@ -34,18 +34,19 @@ mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c) mp_size_t bn, cn; mp_limb_t *bp, *cp; - MPFR_ASSERTN (! MPFR_IS_NAN (b)); - MPFR_ASSERTN (! MPFR_IS_NAN (c)); - - if (MPFR_IS_INF (b)) - return ! MPFR_IS_INF (c); - if (MPFR_IS_INF (c)) - return -1; - - if (MPFR_IS_ZERO (c)) - return ! MPFR_IS_ZERO (b); - if (MPFR_IS_ZERO (b)) - return -1; + if (MPFR_ARE_SINGULAR (b, c)) + { + if (MPFR_IS_INF (b)) + return ! MPFR_IS_INF (c); + else if (MPFR_IS_INF (c)) + return -1; + else if (MPFR_IS_ZERO (c)) + return ! MPFR_IS_ZERO (b); + else if (MPFR_IS_ZERO (b)) + return -1; + else + MPFR_ASSERTN(1); + } be = MPFR_GET_EXP (b); ce = MPFR_GET_EXP (c); @@ -56,8 +57,8 @@ mpfr_cmpabs (mpfr_srcptr b, mpfr_srcptr c) /* exponents are equal */ - bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB; - cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB; + bn = MPFR_LIMB_SIZE(b)-1; + cn = MPFR_LIMB_SIZE(c)-1; bp = MPFR_MANT(b); cp = MPFR_MANT(c); @@ -33,16 +33,17 @@ mpfr_cos (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int K0, K, precy, m, k, l, inexact; mpfr_t r, s; - if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - if (MPFR_IS_ZERO(x)) - { - mpfr_set_ui (y, 1, GMP_RNDN); - return 0; + if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(x)) + return mpfr_set_ui (y, 1, GMP_RNDN); + /* Never reach this */ + MPFR_ASSERTN(1); } precy = MPFR_PREC(y); @@ -39,28 +39,27 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mp_prec_t Nxt = MPFR_PREC(xt); int inexact =0; - if (MPFR_IS_NAN(xt)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (MPFR_IS_NAN(xt)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(xt)) + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + MPFR_RET(0); + } + else if (MPFR_IS_ZERO(xt)) + return mpfr_set_ui(y, 1, rnd_mode); /* cosh(0) = 1 */ + /* Should never reach this code */ + MPFR_ASSERTN(1); } - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(xt)) - { - MPFR_SET_INF(y); - if (MPFR_SIGN(y) < 0) - MPFR_CHANGE_SIGN(y); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(xt)) - return mpfr_set_ui(y,1,rnd_mode); /* cosh(0) = 1 */ mpfr_init2(x,Nxt); - mpfr_set4(x, xt, GMP_RNDN, 1); + mpfr_abs(x, xt, GMP_RNDN); /* General case */ { @@ -85,7 +84,6 @@ mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) mpfr_init (te); mpfr_init (ti); - /* First computation of cosh */ do { @@ -51,61 +51,62 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) * * **************************************************************************/ - if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) + if (MPFR_ARE_SINGULAR(u,v)) { - MPFR_SET_NAN(q); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(q); - - sign_quotient = MPFR_SIGN(u) * MPFR_SIGN(v); - if (MPFR_SIGN(q) != sign_quotient) - MPFR_CHANGE_SIGN(q); - - if (MPFR_IS_INF(u)) - { - if (MPFR_IS_INF(v)) + if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) { MPFR_SET_NAN(q); - MPFR_RET_NAN; + MPFR_RET_NAN; } - else + MPFR_CLEAR_NAN(q); + sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); + MPFR_SET_SIGN(q, sign_quotient); + if (MPFR_IS_INF(u)) { - MPFR_SET_INF(q); - MPFR_RET(0); + if (MPFR_IS_INF(v)) + { + MPFR_SET_NAN(q); + MPFR_RET_NAN; + } + else + { + MPFR_SET_INF(q); + MPFR_RET(0); + } } - } - else - if (MPFR_IS_INF(v)) - { - MPFR_CLEAR_INF(q); - MPFR_SET_ZERO(q); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(q); /* clear Inf flag */ - - if (MPFR_IS_ZERO(v)) - { - if (MPFR_IS_ZERO(u)) + else + if (MPFR_IS_INF(v)) + { + MPFR_CLEAR_INF(q); + MPFR_SET_ZERO(q); + MPFR_RET(0); + } + MPFR_CLEAR_INF(q); /* clear Inf flag */ + if (MPFR_IS_ZERO(v)) { - MPFR_SET_NAN(q); - MPFR_RET_NAN; + if (MPFR_IS_ZERO(u)) + { + MPFR_SET_NAN(q); + MPFR_RET_NAN; + } + else + { + MPFR_SET_INF(q); + MPFR_RET(0); + } } - else + if (MPFR_IS_ZERO(u)) { - MPFR_SET_INF(q); - MPFR_RET(0); + MPFR_SET_ZERO(q); + MPFR_RET(0); } + /* Never reach this !*/ + MPFR_ASSERTN(1); } - - if (MPFR_IS_ZERO(u)) - { - MPFR_SET_ZERO(q); - MPFR_RET(0); - } - + + sign_quotient = MPFR_MULT_SIGN( MPFR_SIGN(u) , MPFR_SIGN(v) ); + MPFR_SET_SIGN(q, sign_quotient); + /************************************************************************** * * * End of the part concerning special values. * @@ -208,9 +209,12 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) { case GMP_RNDU : rnd_mode2 = GMP_RNDD; break; case GMP_RNDD : rnd_mode2 = GMP_RNDU; break; - case GMP_RNDZ : rnd_mode2 = sign_quotient == 1 ? GMP_RNDU : GMP_RNDD; + case GMP_RNDZ : + rnd_mode2 = MPFR_IS_POS_SIGN(sign_quotient) ? GMP_RNDU : GMP_RNDD; + break; + default: + rnd_mode2 = GMP_RNDZ; break; - default : rnd_mode2 = GMP_RNDZ; } can_round2 = mpfr_can_round_raw(qp, qsize + 1, sign_quotient, err + sh + @@ -340,7 +344,8 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) /* Hack : qp[qsize] is 0, 1 or 2, hence if not 0, = 2^(qp[qsize] - 1). */ { near = mpn_rshift(qp, qp, qsize, qp[qsize]); - qp[qsize - 1] |= MPFR_LIMB_HIGHBIT; qexp += qp[qsize]; + qp[qsize - 1] |= MPFR_LIMB_HIGHBIT; + qexp += qp[qsize]; } else { @@ -352,7 +357,8 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) } } - cc = mpfr_round_raw_generic(qp, qp, err, (sign_quotient == -1 ? 1 : 0), + cc = mpfr_round_raw_generic(qp, qp, err, + (MPFR_IS_NEG_SIGN(sign_quotient) ? 1 : 0), MPFR_PREC(q), rnd_mode, &inex, 1); qp += qsize - MPFR_LIMB_SIZE(q); /* 0 or 1 */ @@ -384,15 +390,15 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) if (k >= 0) /* Remainder is nonzero. */ { - if ((rnd_mode == GMP_RNDD && sign_quotient == -1) - || (rnd_mode == GMP_RNDU && sign_quotient == 1)) + if ((rnd_mode == GMP_RNDD && MPFR_IS_NEG_SIGN(sign_quotient)) + || (rnd_mode == GMP_RNDU && MPFR_IS_POS_SIGN(sign_quotient))) /* Rounding to infinity. */ { - inex = sign_quotient; + inex = MPFR_FROM_SIGN_TO_INT( sign_quotient ); cc = 1; } /* rounding to zero. */ - else inex = -sign_quotient; + else inex = -MPFR_FROM_SIGN_TO_INT( sign_quotient ); } } else /* We might have to correct an even rounding if remainder @@ -414,7 +420,8 @@ mpfr_div (mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mp_rnd_t rnd_mode) if (k >= 0) /* In fact the quotient is larger than expected */ { - inex = sign_quotient; /* To infinity, finally. */ + inex = MPFR_FROM_SIGN_TO_INT( sign_quotient ); + /* To infinity, finally. */ cc = 1; } } @@ -34,35 +34,24 @@ mpfr_eq (mpfr_srcptr u, mpfr_srcptr v, unsigned long int n_bits) mp_srcptr up, vp; mp_size_t usize, vsize, size, i; mp_exp_t uexp, vexp; - int usign, k; + int k; - if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) - return 0; /* non equal */ - - usign = MPFR_SIGN(u); - - if (MPFR_IS_INF(u)) - return MPFR_IS_INF(v) && usign == MPFR_SIGN(v); /* +Inf = +Inf */ - else if (MPFR_IS_INF(v)) - return 0; /* +Inf != -Inf */ - - /* 1. Are the signs different? */ - if (usign == MPFR_SIGN(v)) + if (MPFR_ARE_SINGULAR(u, v)) { - /* U and V are both non-negative or both negative. */ - if (MPFR_IS_ZERO(u)) - return MPFR_IS_ZERO(v); /* 0 = 0 */ - if (MPFR_IS_ZERO(v)) - return MPFR_IS_ZERO(u); /* 0 = 0 */ - - /* Fall out. */ - } - else - { - /* Either U or V is negative, but not both. */ - return MPFR_IS_ZERO(u) && MPFR_IS_ZERO(v); + if (MPFR_IS_NAN(u) || MPFR_IS_NAN(v)) + return 0; /* non equal */ + else if (MPFR_IS_INF(u) && MPFR_IS_INF(v)) + return (MPFR_SIGN(u) == MPFR_SIGN(v)); + else if (MPFR_IS_ZERO(u) && MPFR_IS_ZERO(v)) + return 1; + else + return 0; } + /* 1. Are the signs different? */ + if (MPFR_SIGN(u) != MPFR_SIGN(v)) + return 0; + uexp = MPFR_GET_EXP (u); vexp = MPFR_GET_EXP (v); @@ -44,26 +44,27 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) double n = (double) MPFR_PREC(y); int inex; - if (MPFR_IS_NAN(x)) + sign_x = MPFR_SIGN (x); + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + if (MPFR_IS_INF(x)) /* erf(+inf) = +1, erf(-inf) = -1 */ + return mpfr_set_si (y, MPFR_FROM_SIGN_TO_INT(sign_x), GMP_RNDN); + if (MPFR_IS_ZERO(x)) /* erf(+0) = +0, erf(-0) = -0 */ + return mpfr_set (y, x, GMP_RNDN); /* should keep the sign of x */ + MPFR_ASSERTN(1); } - sign_x = MPFR_SIGN (x); - - if (MPFR_IS_INF(x)) /* erf(+inf) = +1, erf(-inf) = -1 */ - return mpfr_set_si (y, sign_x, GMP_RNDN); - - if (MPFR_IS_ZERO(x)) /* erf(+0) = +0, erf(-0) = -0 */ - return mpfr_set (y, x, GMP_RNDN); /* should keep the sign of x */ - /* now x is neither NaN, Inf nor 0 */ xf = mpfr_get_d (x, GMP_RNDN); xf = xf * xf; /* xf ~ x^2 */ - if (sign_x > 0) + if (MPFR_IS_POS_SIGN(sign_x)) rnd2 = rnd_mode; else { @@ -96,7 +97,7 @@ mpfr_erf (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) inex = mpfr_erf_0 (y, x, rnd2); } - if (sign_x < 0) + if (MPFR_IS_NEG_SIGN(sign_x)) { MPFR_CHANGE_SIGN (y); return - inex; @@ -38,34 +38,33 @@ mpfr_exp (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int expx, precy; double d; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) > 0) + if (MPFR_IS_NAN(x)) { - MPFR_SET_INF(y); + MPFR_SET_NAN(y); + MPFR_RET_NAN; } - else + MPFR_CLEAR_NAN(y); + if (MPFR_IS_INF(x)) { - MPFR_CLEAR_INF(y); - MPFR_SET_ZERO(y); + if (MPFR_IS_POS(x)) + { + MPFR_SET_INF(y); + } + else + { + MPFR_CLEAR_INF(y); + MPFR_SET_ZERO(y); + } + MPFR_SET_POS(y); + MPFR_RET(0); } - MPFR_SET_POS(y); - MPFR_RET(0); + MPFR_CLEAR_INF(y); + if (MPFR_IS_ZERO(x)) + return mpfr_set_ui (y, 1, GMP_RNDN); } - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(x)) - return mpfr_set_ui (y, 1, GMP_RNDN); - expx = MPFR_GET_EXP (x); precy = MPFR_PREC(y); @@ -36,33 +36,33 @@ mpfr_exp2 (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int inexact; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + if (MPFR_IS_INF(x)) + { + if (MPFR_IS_POS(x)) + { + MPFR_SET_INF(y); + } + else + { + MPFR_CLEAR_INF(y); + MPFR_SET_ZERO(y); + } + MPFR_SET_POS(y); + MPFR_RET(0); + } + /* 2^0 = 1 */ + if (MPFR_IS_ZERO(x)) + return mpfr_set_ui (y, 1, rnd_mode); } - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) > 0) - { - MPFR_SET_INF(y); - } - else - { - MPFR_CLEAR_INF(y); - MPFR_SET_ZERO(y); - } - MPFR_SET_POS(y); - MPFR_RET(0); - } - - /* 2^0 = 1 */ - if (MPFR_IS_ZERO(x)) - return mpfr_set_ui (y, 1, rnd_mode); - /* since the smallest representable non-zero float is 1/2*2^__gmpfr_emin, if x < __gmpfr_emin - 1, the result is either 1/2*2^__gmpfr_emin or 0 */ MPFR_ASSERTN(MPFR_EMIN_MIN - 2 >= LONG_MIN); @@ -34,34 +34,33 @@ mpfr_expm1 (mpfr_ptr y, mpfr_srcptr x , mp_rnd_t rnd_mode) { int inexact = 0; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - /* check for inf or -inf (expm1(-inf)=-1) */ - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) > 0) - { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - return 0; - } - else - return mpfr_set_si(y, -1, rnd_mode); - } - - MPFR_CLEAR_INF(y); - - if(MPFR_IS_ZERO(x)) - { - MPFR_SET_ZERO(y); /* expm1(+/- 0) = +/- 0 */ - MPFR_SET_SAME_SIGN(y,x); - MPFR_RET(0); + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + /* check for inf or -inf (expm1(-inf)=-1) */ + if (MPFR_IS_INF(x)) + { + if (MPFR_IS_POS(x)) + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + return 0; + } + else + return mpfr_set_si(y, -1, rnd_mode); + } + MPFR_CLEAR_INF(y); + if(MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); /* expm1(+/- 0) = +/- 0 */ + MPFR_SET_SAME_SIGN(y,x); + MPFR_RET(0); + } } /* General case */ @@ -40,76 +40,72 @@ mpfr_fma (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z, mpfr_t u; /* particular cases */ - if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) || + MPFR_IS_SINGULAR(y) || + MPFR_IS_SINGULAR(z) )) { - MPFR_SET_NAN(s); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) - { - /* cases Inf*0+z, 0*Inf+z, Inf-Inf */ - if ((MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) || - (MPFR_IS_FP(x) && MPFR_IS_ZERO(x)) || - (MPFR_IS_INF(z) && ((MPFR_SIGN(x) * MPFR_SIGN(y)) != MPFR_SIGN(z)))) - { - MPFR_SET_NAN(s); - MPFR_RET_NAN; - } - + if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z)) + { + MPFR_SET_NAN(s); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) + { + /* cases Inf*0+z, 0*Inf+z, Inf-Inf */ + if ((MPFR_IS_ZERO(y)) || + (MPFR_IS_ZERO(x)) || + (MPFR_IS_INF(z) && + ((MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y))) != MPFR_SIGN(z)))) + { + MPFR_SET_NAN(s); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(s); + if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */ + { + MPFR_SET_INF(s); + MPFR_SET_SAME_SIGN(s, z); + MPFR_RET(0); + } + else /* z is finite */ + { + MPFR_SET_INF(s); + MPFR_SET_SIGN(s, MPFR_MULT_SIGN(MPFR_SIGN(x) , MPFR_SIGN(y))); + MPFR_RET(0); + } + } MPFR_CLEAR_NAN(s); - - if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */ - { - MPFR_SET_INF(s); - MPFR_SET_SAME_SIGN(s, z); - MPFR_RET(0); - } - else /* z is finite */ - { - MPFR_SET_INF(s); - if (MPFR_SIGN(s) != (MPFR_SIGN(x) * MPFR_SIGN(y))) - MPFR_CHANGE_SIGN(s); - MPFR_RET(0); - } - } - - MPFR_CLEAR_NAN(s); - - /* now x and y are finite */ - if (MPFR_IS_INF(z)) - { - MPFR_SET_INF(s); - MPFR_SET_SAME_SIGN(s, z); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(s); - - if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y)) - { + /* now x and y are finite */ + if (MPFR_IS_INF(z)) + { + MPFR_SET_INF(s); + MPFR_SET_SAME_SIGN(s, z); + MPFR_RET(0); + } + MPFR_CLEAR_INF(s); + if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y)) + { + if (MPFR_IS_ZERO(z)) + { + int sign_p, sign_z; + sign_p = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) ); + sign_z = MPFR_SIGN(z); + MPFR_SET_SIGN(s,(rnd_mode != GMP_RNDD ? + ((MPFR_IS_NEG_SIGN(sign_p) && MPFR_IS_NEG_SIGN(sign_z)) + ? -1 : 1) : + ((MPFR_IS_POS_SIGN(sign_p) && MPFR_IS_POS_SIGN(sign_z)) + ? 1 : -1))); + MPFR_SET_ZERO(s); + MPFR_RET(0); + } + else + return mpfr_set (s, z, rnd_mode); + } if (MPFR_IS_ZERO(z)) - { - int sign_p, sign_z; - sign_p = MPFR_SIGN(x) * MPFR_SIGN(y); - sign_z = MPFR_SIGN(z); - if (MPFR_SIGN(s) != - (rnd_mode != GMP_RNDD ? - ((sign_p < 0 && sign_z < 0) ? -1 : 1) : - ((sign_p > 0 && sign_z > 0) ? 1 : -1))) - { - MPFR_CHANGE_SIGN(s); - } - MPFR_SET_ZERO(s); - MPFR_RET(0); - } - else - return mpfr_set (s, z, rnd_mode); + return mpfr_mul (s, x, y, rnd_mode); + MPFR_ASSERTN(1); } - if (MPFR_IS_ZERO(z)) - return mpfr_mul (s, x, y, rnd_mode); - /* if we take prec(u) >= prec(x) + prec(y), the product u <- x*y is always exact */ mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y)); @@ -36,13 +36,13 @@ mpfr_frac (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) mpfr_t tmp; mpfr_ptr t; - if (MPFR_IS_NAN(u)) + /* Special cases */ + if (MPFR_UNLIKELY(MPFR_IS_NAN(u))) { MPFR_SET_NAN(r); MPFR_RET_NAN; } - - if (MPFR_IS_INF(u) || mpfr_integer_p (u)) + else if (MPFR_UNLIKELY(MPFR_IS_INF(u) || mpfr_integer_p (u))) { MPFR_CLEAR_FLAGS(r); MPFR_SET_SAME_SIGN(r, u); @@ -60,34 +60,36 @@ mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode) int inex; /* Trivial cases */ - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(gamma); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) < 0) - { - MPFR_SET_NAN(gamma); - MPFR_RET_NAN; - } - else - { - MPFR_CLEAR_NAN(gamma); - MPFR_SET_INF(gamma); - MPFR_SET_POS(gamma); - return 0; /* exact */ - } - } - - if (MPFR_IS_ZERO(x)) - { - MPFR_CLEAR_NAN(gamma); - MPFR_SET_INF(gamma); - MPFR_SET_SAME_SIGN(gamma, x); - return 0; /* exact */ + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(gamma); + MPFR_RET_NAN; + } + if (MPFR_IS_INF(x)) + { + if (MPFR_IS_NEG(x)) + { + MPFR_SET_NAN(gamma); + MPFR_RET_NAN; + } + else + { + MPFR_CLEAR_NAN(gamma); + MPFR_SET_INF(gamma); + MPFR_SET_POS(gamma); + return 0; /* exact */ + } + } + if (MPFR_IS_ZERO(x)) + { + MPFR_CLEAR_NAN(gamma); + MPFR_SET_INF(gamma); + MPFR_SET_SAME_SIGN(gamma, x); + return 0; /* exact */ + } + MPFR_ASSERTN(1); } /* Set x_p=x if x> 1 else set x_p=2-x */ diff --git a/gammaPiAGMformula.c b/gammaPiAGMformula.c index 82f590c6e..acd467d50 100644 --- a/gammaPiAGMformula.c +++ b/gammaPiAGMformula.c @@ -72,21 +72,26 @@ mpfr_gamma (gamma, x, rnd_mode) int k; /* Trivial cases */ - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(gamma); - return 1; - } - if (!MPFR_NOTZERO(x)) - { - MPFR_SET_INF(gamma); - return 1; - } - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(gamma); - return 1; + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(gamma); + return 1; + } + if (MPFR_ISZERO(x)) + { + MPFR_SET_INF(gamma); + return 1; + } + if (MPFR_IS_INF(x)) + { + MPFR_SET_INF(gamma); + return 1; + } + MPFR_ASSERTN(1); } + /* Set x_p=x if x> 1 else set x_p=2-x */ prec_gamma = MPFR_PREC(gamma); compared = mpfr_cmp_ui(x, 1); @@ -44,30 +44,28 @@ mpfr_hypot (mpfr_ptr z, mpfr_srcptr x , mpfr_srcptr y , mp_rnd_t rnd_mode) mp_exp_unsigned_t diff_exp; /* particular cases */ - - if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y)) - { - MPFR_SET_NAN(z); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(z); - - if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) + if (MPFR_ARE_SINGULAR(x,y)) { - MPFR_SET_INF(z); - MPFR_SET_POS(z); - MPFR_RET(0); + if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y)) + { + MPFR_SET_NAN(z); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(z); + if (MPFR_IS_INF(x) || MPFR_IS_INF(y)) + { + MPFR_SET_INF(z); + MPFR_SET_POS(z); + MPFR_RET(0); + } + MPFR_CLEAR_INF(z); + if (MPFR_IS_ZERO(x)) + return mpfr_abs (z, y, rnd_mode); + if (MPFR_IS_ZERO(y)) + return mpfr_abs (z, x, rnd_mode); + MPFR_ASSERTN(1); } - MPFR_CLEAR_INF(z); - - if (MPFR_IS_ZERO(x)) - return mpfr_abs (z, y, rnd_mode); - - if (MPFR_IS_ZERO(y)) - return mpfr_abs (z, x, rnd_mode); - if (mpfr_cmpabs (x, y) < 0) { mpfr_srcptr t; @@ -51,50 +51,51 @@ mpfr_log (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) double ref; TMP_DECL(marker); - /* If a is NaN, the result is NaN */ - if (MPFR_IS_NAN(a)) + /* Special cases */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) )) { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(r); - - /* check for infinity before zero */ - if (MPFR_IS_INF(a)) - { - if (MPFR_SIGN(a) < 0) /* log(-Inf) = NaN */ + /* If a is NaN, the result is NaN */ + if (MPFR_IS_NAN(a)) { MPFR_SET_NAN(r); - MPFR_RET_NAN; + MPFR_RET_NAN; } - else /* log(+Inf) = +Inf */ + MPFR_CLEAR_NAN(r); + /* check for infinity before zero */ + if (MPFR_IS_INF(a)) + { + if (MPFR_IS_NEG(a)) + /* log(-Inf) = NaN */ + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + else /* log(+Inf) = +Inf */ + { + MPFR_SET_INF(r); + MPFR_SET_POS(r); + MPFR_RET(0); + } + } + /* Now we can clear the flags without damage even if r == a */ + MPFR_CLEAR_INF(r); + if (MPFR_IS_ZERO(a)) { MPFR_SET_INF(r); - MPFR_SET_POS(r); - MPFR_RET(0); + MPFR_SET_NEG(r); + MPFR_RET(0); /* log(0) is an exact -infinity */ } } - - /* Now we can clear the flags without damage even if r == a */ - MPFR_CLEAR_INF(r); - - if (MPFR_IS_ZERO(a)) - { - MPFR_SET_INF(r); - MPFR_SET_NEG(r); - MPFR_RET(0); /* log(0) is an exact -infinity */ - } - + /* If a is negative, the result is NaN */ - if (MPFR_SIGN(a) < 0) + if (MPFR_UNLIKELY( MPFR_IS_NEG(a) )) { MPFR_SET_NAN(r); MPFR_RET_NAN; } /* If a is 1, the result is 0 */ - if (mpfr_cmp_ui (a, 1) == 0) + if (MPFR_UNLIKELY( mpfr_cmp_ui (a, 1) == 0) ) { MPFR_SET_ZERO(r); MPFR_SET_POS(r); @@ -36,42 +36,42 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) int inexact = 0; /* If a is NaN, the result is NaN */ - if (MPFR_IS_NAN(a)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) )) { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(r); - - /* check for infinity before zero */ - if (MPFR_IS_INF(a)) - { - if (MPFR_SIGN(a) < 0) /* log10(-Inf) = NaN */ - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - else /* log10(+Inf) = +Inf */ - { - MPFR_SET_INF(r); - MPFR_SET_POS(r); - MPFR_RET(0); /* exact */ - } - } - - /* Now we can clear the flags without damage even if r == a */ - MPFR_CLEAR_INF(r); - - if (MPFR_IS_ZERO(a)) - { - MPFR_SET_INF(r); - MPFR_SET_NEG(r); - MPFR_RET(0); /* log10(0) is an exact -infinity */ + if (MPFR_IS_NAN(a)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(r); + /* check for infinity before zero */ + if (MPFR_IS_INF(a)) + { + if (MPFR_IS_NEG(a)) + /* log10(-Inf) = NaN */ + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + else /* log10(+Inf) = +Inf */ + { + MPFR_SET_INF(r); + MPFR_SET_POS(r); + MPFR_RET(0); /* exact */ + } + } + /* Now we can clear the flags without damage even if r == a */ + MPFR_CLEAR_INF(r); + if (MPFR_IS_ZERO(a)) + { + MPFR_SET_INF(r); + MPFR_SET_NEG(r); + MPFR_RET(0); /* log10(0) is an exact -infinity */ + } } /* If a is negative, the result is NaN */ - if (MPFR_SIGN(a) < 0) + if (MPFR_UNLIKELY( MPFR_IS_NEG(a) )) { MPFR_SET_NAN(r); MPFR_RET_NAN; @@ -34,53 +34,54 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) { int comp, inexact = 0; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - /* check for inf or -inf (result is not defined) */ - if (MPFR_IS_INF(x)) - { - if (MPFR_SIGN(x) > 0) - { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - MPFR_RET(0); - } - else - { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(y); + /* check for inf or -inf (result is not defined) */ + if (MPFR_IS_INF(x)) + { + if (MPFR_IS_POS(x)) + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + MPFR_RET(0); + } + else + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + } + if (MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); /* log1p(+/- 0) = +/- 0 */ + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + MPFR_ASSERTN(1); } - + comp = mpfr_cmp_si(x,-1); /* x<-1 undefined */ - if (comp < 0) + if (MPFR_UNLIKELY(comp <= 0)) { + if (comp == 0) + /* x=0: log1p(-1)=-inf (division by zero) */ + { + MPFR_SET_INF(y); + MPFR_SET_POS(y); + MPFR_RET_NAN; + } MPFR_SET_NAN(y); MPFR_RET_NAN; } - /* x=0: log1p(-1)=-inf (division by zero) */ - if (comp == 0) - { - MPFR_SET_INF(y); - MPFR_SET_POS(y); - MPFR_RET_NAN; - } - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(x)) - { - MPFR_SET_ZERO(y); /* log1p(+/- 0) = +/- 0 */ - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); - } + MPFR_CLEAR_FLAGS(y); /* General case */ { @@ -35,43 +35,43 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode) { int inexact = 0; - /* If a is NaN, the result is NaN */ - if (MPFR_IS_NAN(a)) - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(r); - - /* check for infinity before zero */ - if (MPFR_IS_INF(a)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) )) { - if (MPFR_SIGN(a) < 0) /* log(-Inf) = NaN */ - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - else /* log(+Inf) = +Inf */ - { - MPFR_SET_INF(r); - MPFR_SET_POS(r); - MPFR_RET(0); - } + /* If a is NaN, the result is NaN */ + if (MPFR_IS_NAN(a)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(r); + /* check for infinity before zero */ + if (MPFR_IS_INF(a)) + { + if (MPFR_IS_NEG(a)) + /* log(-Inf) = NaN */ + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + else /* log(+Inf) = +Inf */ + { + MPFR_SET_INF(r); + MPFR_SET_POS(r); + MPFR_RET(0); + } + } + /* Now we can clear the flags without damage even if r == a */ + MPFR_CLEAR_INF(r); + if (MPFR_IS_ZERO(a)) + { + MPFR_SET_INF(r); + MPFR_SET_NEG(r); + MPFR_RET(0); /* log2(0) is an exact -infinity */ + } } - - /* Now we can clear the flags without damage even if r == a */ - MPFR_CLEAR_INF(r); - - if (MPFR_IS_ZERO(a)) - { - MPFR_SET_INF(r); - MPFR_SET_NEG(r); - MPFR_RET(0); /* log2(0) is an exact -infinity */ - } - + /* If a is negative, the result is NaN */ - if (MPFR_SIGN(a) < 0) + if (MPFR_UNLIKELY(MPFR_IS_NEG(a))) { MPFR_SET_NAN(r); MPFR_RET_NAN; @@ -33,27 +33,26 @@ MA 02111-1307, USA. */ int mpfr_min (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) { - if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) ) - { - MPFR_SET_NAN(z); - MPFR_RET_NAN; - } - MPFR_CLEAR_NAN(z); - - if (MPFR_IS_NAN(x)) - return mpfr_set(z, y, rnd_mode); - - if (MPFR_IS_NAN(y)) - return mpfr_set(z, x, rnd_mode); - - if (MPFR_IS_FP(x) && MPFR_IS_ZERO(x) && MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) + if (MPFR_ARE_SINGULAR(x,y)) { - if (MPFR_SIGN(x) < 0) - return mpfr_set(z, x, rnd_mode); - else - return mpfr_set(z, y, rnd_mode); + if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) ) + { + MPFR_SET_NAN(z); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(z); + if (MPFR_IS_NAN(x)) + return mpfr_set(z, y, rnd_mode); + if (MPFR_IS_NAN(y)) + return mpfr_set(z, x, rnd_mode); + if (MPFR_IS_ZERO(x) && MPFR_IS_ZERO(y)) + { + if (MPFR_IS_NEG(x)) + return mpfr_set(z, x, rnd_mode); + else + return mpfr_set(z, y, rnd_mode); + } } - if (mpfr_cmp(x,y) <= 0) return mpfr_set(z, x, rnd_mode); else @@ -69,25 +68,25 @@ mpfr_min (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) int mpfr_max (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) { - if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) ) - { - MPFR_SET_NAN(z); - MPFR_RET_NAN; - } - MPFR_CLEAR_NAN(z); - - if (MPFR_IS_NAN(x)) - return mpfr_set(z, y, rnd_mode); - - if (MPFR_IS_NAN(y)) - return mpfr_set(z, x, rnd_mode); - - if (MPFR_IS_FP(x) && MPFR_IS_ZERO(x) && MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) + if (MPFR_ARE_SINGULAR(x,y)) { - if (MPFR_SIGN(x) < 0) - return mpfr_set(z, y, rnd_mode); - else - return mpfr_set(z, x, rnd_mode); + if (MPFR_IS_NAN(x) && MPFR_IS_NAN(y) ) + { + MPFR_SET_NAN(z); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(z); + if (MPFR_IS_NAN(x)) + return mpfr_set(z, y, rnd_mode); + if (MPFR_IS_NAN(y)) + return mpfr_set(z, x, rnd_mode); + if (MPFR_IS_ZERO(x) && MPFR_IS_ZERO(y)) + { + if (MPFR_IS_NEG(x)) + return mpfr_set(z, y, rnd_mode); + else + return mpfr_set(z, x, rnd_mode); + } } if (mpfr_cmp(x,y) <= 0) diff --git a/mpfr-impl.h b/mpfr-impl.h index 9f1e8380f..5fab45fa5 100644 --- a/mpfr-impl.h +++ b/mpfr-impl.h @@ -93,7 +93,8 @@ typedef unsigned long int mp_size_unsigned_t; /* MPFR_ASSERTN(expr): assertions that should always be checked */ /* #define MPFR_ASSERTN(expr) ASSERT_ALWAYS(expr) */ -#define MPFR_ASSERTN(expr) ((void) ((expr) || (ASSERT_FAIL (expr), 0))) +#define MPFR_ASSERTN(expr) \ + ((void) ((MPFR_UNLIKELY(expr)) || (ASSERT_FAIL (expr), 0))) /* MPFR_ASSERTD(expr): assertions that should be checked when testing */ /* #define MPFR_ASSERTD(expr) ASSERT(expr) */ @@ -240,6 +241,7 @@ long double __gmpfr_longdouble_volatile __GMP_PROTO ((long double)) ATTRIBUTE_CO /* Definition of the special values of the exponent */ /* Clear flags macros are still defined and should be still used */ /* since functions shouldn't rely on a specific format */ +/* NOTE: Well, until a real spec of how to use them isn't created, it is unusable */ #define MPFR_PREC(x) ((x)->_mpfr_prec) #define MPFR_EXP(x) ((x)->_mpfr_exp) @@ -269,23 +271,37 @@ long double __gmpfr_longdouble_volatile __GMP_PROTO ((long double)) ATTRIBUTE_CO #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_REAL_FP(x) (!MPFR_IS_SINGULAR(x)) +#define MPFR_IS_PURE_FP(x) (!MPFR_IS_SINGULAR(x)) -/* FIXME: NOTZERO real usefull ? */ +#define MPFR_ARE_SINGULAR(x,y) \ + (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x)) || MPFR_UNLIKELY(MPFR_IS_SINGULAR(y))) + +/* Without zeros */ #define MPFR_ISNONNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) > 0) #define MPFR_ISNEG(x) (MPFR_NOTZERO((x)) && MPFR_SIGN(x) < 0) +/* With signed zeros */ +#define MPFR_IS_NEG(x) (MPFR_SIGN(x) < 0) +#define MPFR_IS_POS(x) (MPFR_SIGN(x) > 0) + #define MPFR_SET_POS(x) (MPFR_SIGN(x) = MPFR_SIGN_POS) #define MPFR_SET_NEG(x) (MPFR_SIGN(x) = MPFR_SIGN_NEG) #define MPFR_CHANGE_SIGN(x) (MPFR_SIGN(x) = -MPFR_SIGN(x)) #define MPFR_SET_SAME_SIGN(x, y) (MPFR_SIGN(x) = MPFR_SIGN(y)) +#define MPFR_SET_OPPOSITE_SIGN(x, y) (MPFR_SIGN(x) = - MPFR_SIGN(y)) #define MPFR_CHECK_SIGN(s) \ (MPFR_ASSERTD((s) == MPFR_SIGN_POS || (s) == MPFR_SIGN_NEG)) #define MPFR_SET_SIGN(x, s) \ (MPFR_CHECK_SIGN(s), MPFR_SIGN(x) = s) -#define MPFR_MULT_SIGN(x, s) \ +#define MPFR_IS_POS_SIGN(s1) (s1 > 0) +#define MPFR_IS_NEG_SIGN(s1) (s1 < 0) +#define MPFR_MULT_SIGN(s1, s2) \ + ((s1) * (s2)) +#define MPFR_SET_MULT_SIGN(x, s) \ (MPFR_CHECK_SIGN(s), MPFR_SIGN(x) *= s) +/* Transform a sign to 1 or -1 */ +#define MPFR_FROM_SIGN_TO_INT(s) (s) /* Special inexact value */ #define MPFR_EVEN_INEX 2 @@ -34,10 +34,10 @@ MA 02111-1307, USA. */ /* Definition of rounding modes */ typedef enum { - MPFR_RNDN=0, - MPFR_RNDZ=1, - MPFR_RNDU=2, - MPFR_RNDD=3 + GMP_RNDN=0, + GMP_RNDZ=1, + GMP_RNDU=2, + GMP_RNDD=3 } mpfr_rnd_t; /* Flags of __gmpfr_flags */ @@ -84,10 +84,6 @@ typedef struct { } __mpfr_struct; /* Compatibility with previous versions of MPFR */ -#define GMP_RNDN MPFR_RNDN -#define GMP_RNDZ MPFR_RNDZ -#define GMP_RNDU MPFR_RNDU -#define GMP_RNDD MPFR_RNDD #define mp_rnd_t mpfr_rnd_t #define mp_prec_t mpfr_prec_t @@ -35,59 +35,57 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) mp_size_t an, bn, cn, tn, k; TMP_DECL(marker); - /* deal with NaN and zero */ - if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + /* deal with special cases */ + if (MPFR_ARE_SINGULAR(b,c)) { - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(a); - - sign_product = MPFR_SIGN(b) * MPFR_SIGN(c); - - if (MPFR_IS_INF(b)) - { - if (MPFR_IS_INF(c) || MPFR_NOTZERO(c)) - { - if (MPFR_SIGN(a) != sign_product) - MPFR_CHANGE_SIGN(a); - MPFR_SET_INF(a); - MPFR_RET(0); /* exact */ - } - else + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) { MPFR_SET_NAN(a); - MPFR_RET_NAN; + MPFR_RET_NAN; } - } - else if (MPFR_IS_INF(c)) - { - if (MPFR_NOTZERO(b)) - { - if (MPFR_SIGN(a) != sign_product) - MPFR_CHANGE_SIGN(a); - MPFR_SET_INF(a); - MPFR_RET(0); /* exact */ - } - else + MPFR_CLEAR_NAN(a); + sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) ); + if (MPFR_IS_INF(b)) { - MPFR_SET_NAN(a); - MPFR_RET_NAN; + if (MPFR_IS_INF(c) || MPFR_NOTZERO(c)) + { + MPFR_SET_SIGN(a,sign_product); + MPFR_SET_INF(a); + MPFR_RET(0); /* exact */ + } + else + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } } + else if (MPFR_IS_INF(c)) + { + if (MPFR_NOTZERO(b)) + { + MPFR_SET_SIGN(a, sign_product); + MPFR_SET_INF(a); + MPFR_RET(0); /* exact */ + } + else + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + } + MPFR_CLEAR_INF(a); /* clear Inf flag */ + if (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c)) + { + MPFR_SET_SIGN(a, sign_product); + MPFR_SET_ZERO(a); + MPFR_RET(0); /* 0 * 0 is exact */ + } + /* Should never reach here */ + MPFR_ASSERTN(1); } - - MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c)); - MPFR_CLEAR_INF(a); /* clear Inf flag */ - - if (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c)) - { - if (MPFR_SIGN(a) != sign_product) - MPFR_CHANGE_SIGN(a); - MPFR_SET_ZERO(a); - MPFR_RET(0); /* 0 * 0 is exact */ - } - + MPFR_CLEAR_FLAGS(a); + sign_product = MPFR_MULT_SIGN( MPFR_SIGN(b) , MPFR_SIGN(c) ); + bx = MPFR_GET_EXP (b); cx = MPFR_GET_EXP (c); /* Note: the exponent of the exact result will be e = bx + cx + ec with @@ -133,7 +131,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) tmp += k - tn; if (b1 == 0) mpn_lshift (tmp, tmp, tn, 1); - cc = mpfr_round_raw (ap, tmp, bq + cq, sign_product < 0, aq, + cc = mpfr_round_raw (ap, tmp, bq + cq, MPFR_IS_NEG_SIGN(sign_product), aq, rnd_mode, &inexact); if (cc) /* cc = 1 ==> result is a power of two */ ap[an-1] = MPFR_LIMB_HIGHBIT; @@ -157,8 +155,7 @@ mpfr_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) } MPFR_SET_EXP (a, ax); - if (MPFR_SIGN(a) != sign_product) - MPFR_CHANGE_SIGN(a); + MPFR_SET_SIGN(a, sign_product); return inexact; } @@ -33,40 +33,49 @@ mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) int cnt, c, inexact; TMP_DECL(marker); - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(x)) - { - if (u != 0) + if (MPFR_IS_NAN(x)) { - MPFR_CLEAR_FLAGS(y); - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); /* infinity is exact */ + MPFR_SET_NAN(y); + MPFR_RET_NAN; } - else /* 0 * infinity */ + else if (MPFR_IS_INF(x)) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (u != 0) + { + MPFR_CLEAR_FLAGS(y); + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); /* infinity is exact */ + } + else /* 0 * infinity */ + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + } + else if (MPFR_IS_ZERO(x)) + { + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); /* zero is exact */ } + else + MPFR_ASSERTN(1); } - - MPFR_CLEAR_FLAGS(y); - - if (u == 0 || MPFR_IS_ZERO(x)) + else if (MPFR_UNLIKELY(u <= 1)) { - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); /* zero is exact */ + if (u == 0) + { + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); /* zero is exact */ + } + else + return mpfr_set (y, x, rnd_mode); } - if (u == 1) - return mpfr_set (y, x, rnd_mode); - TMP_MARK(marker); yp = MPFR_MANT(y); yn = (MPFR_PREC(y) - 1) / BITS_PER_MP_LIMB + 1; @@ -112,8 +121,8 @@ mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mp_rnd_t rnd_mode) TMP_FREE(marker); - if (__gmpfr_emax < MPFR_EMAX_MIN + cnt || - MPFR_GET_EXP (x) > __gmpfr_emax - cnt) + if (MPFR_UNLIKELY(__gmpfr_emax < MPFR_EMAX_MIN + cnt || + MPFR_GET_EXP (x) > __gmpfr_emax - cnt)) return mpfr_set_overflow(y, rnd_mode, MPFR_SIGN(x)); MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt); @@ -29,14 +29,13 @@ MA 02111-1307, USA. */ static void mpfr_nexttozero (mpfr_ptr x) { - if (MPFR_IS_INF(x)) + if (MPFR_UNLIKELY(MPFR_IS_INF(x))) { MPFR_CLEAR_FLAGS(x); mpfr_setmax (x, __gmpfr_emax); return; } - - if (MPFR_IS_ZERO(x)) + else if (MPFR_UNLIKELY( MPFR_IS_ZERO(x) )) { MPFR_CHANGE_SIGN(x); mpfr_setmin (x, __gmpfr_emin); @@ -54,7 +53,7 @@ mpfr_nexttozero (mpfr_ptr x) if (xp[xn-1] >> (BITS_PER_MP_LIMB - 1) == 0) { /* was an exact power of two: not normalized any more */ mp_exp_t exp = MPFR_EXP (x); - if (exp == __gmpfr_emin) + if (MPFR_UNLIKELY(exp == __gmpfr_emin)) MPFR_SET_ZERO(x); else { @@ -71,10 +70,9 @@ mpfr_nexttozero (mpfr_ptr x) static void mpfr_nexttoinf (mpfr_ptr x) { - if (MPFR_IS_INF(x)) + if (MPFR_UNLIKELY(MPFR_IS_INF(x))) return; - - if (MPFR_IS_ZERO(x)) + else if (MPFR_UNLIKELY(MPFR_IS_ZERO(x))) mpfr_setmin (x, __gmpfr_emin); else { @@ -88,7 +86,7 @@ mpfr_nexttoinf (mpfr_ptr x) if (mpn_add_1 (xp, xp, xn, MP_LIMB_T_ONE << sh)) /* got 1.0000... */ { mp_exp_t exp = MPFR_EXP (x); - if (exp == __gmpfr_emax) + if (MPFR_UNLIKELY(exp == __gmpfr_emax)) MPFR_SET_INF(x); else { @@ -102,13 +100,12 @@ mpfr_nexttoinf (mpfr_ptr x) void mpfr_nextabove (mpfr_ptr x) { - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY(MPFR_IS_NAN(x))) { __gmpfr_flags |= MPFR_FLAGS_NAN; return; } - - if (MPFR_SIGN(x) < 0) + if (MPFR_IS_NEG(x)) mpfr_nexttozero (x); else mpfr_nexttoinf (x); @@ -117,13 +114,13 @@ mpfr_nextabove (mpfr_ptr x) void mpfr_nextbelow (mpfr_ptr x) { - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY(MPFR_IS_NAN(x))) { __gmpfr_flags |= MPFR_FLAGS_NAN; return; } - if (MPFR_SIGN(x) < 0) + if (MPFR_IS_NEG(x)) mpfr_nexttoinf (x); else mpfr_nexttozero (x); @@ -134,7 +131,7 @@ mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y) { int s; - if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y)) + if (MPFR_UNLIKELY(MPFR_IS_NAN(x) || MPFR_IS_NAN(y))) { __gmpfr_flags |= MPFR_FLAGS_NAN; return; @@ -143,7 +140,7 @@ mpfr_nexttoward (mpfr_ptr x, mpfr_srcptr y) s = mpfr_cmp (x, y); if (s == 0) return; - if (s < 0) + else if (s < 0) mpfr_nextabove (x); else mpfr_nextbelow (x); @@ -150,109 +150,105 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) int inexact = 0; /* pow(x, ±0) returns 1 for any x, even a NaN. */ - if (MPFR_IS_FP(y) && MPFR_IS_ZERO(y)) + if (MPFR_UNLIKELY( MPFR_IS_ZERO(y) )) return mpfr_set_ui (z, 1, GMP_RNDN); - if (MPFR_IS_NAN(x)) - { - MPFR_SET_NAN(z); - MPFR_RET_NAN; - } - - if (MPFR_IS_NAN(y)) - { - /* pow(+1, NaN) returns 1. */ - if (MPFR_IS_FP(x) && mpfr_cmp_ui(x, 1) == 0) - return mpfr_set_ui (z, 1, GMP_RNDN); - - MPFR_SET_NAN(z); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(y)) + if (MPFR_ARE_SINGULAR(x,y)) { + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(z); + MPFR_RET_NAN; + } + else if (MPFR_IS_NAN(y)) + { + /* pow(+1, NaN) returns 1. */ + if (mpfr_cmp_ui(x, 1) == 0) + return mpfr_set_ui (z, 1, GMP_RNDN); + MPFR_SET_NAN(z); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(y)) + { + if (MPFR_IS_INF(x)) + { + MPFR_CLEAR_FLAGS(z); + if (MPFR_IS_POS(y)) + MPFR_SET_INF(z); + else + MPFR_SET_ZERO(z); + MPFR_SET_POS(z); + MPFR_RET(0); + } + else + { + mpfr_t one; + int cmp; + mpfr_init2(one, BITS_PER_MP_LIMB); + mpfr_set_ui(one, 1, GMP_RNDN); + cmp = mpfr_cmpabs(x, one) * MPFR_FROM_SIGN_TO_INT(MPFR_SIGN(y)); + mpfr_clear(one); + MPFR_CLEAR_FLAGS(z); + MPFR_SET_POS(z); + if (cmp > 0) + { + /* Return +inf. */ + MPFR_SET_INF(z); + MPFR_RET(0); + } + else if (cmp < 0) + { + /* Return +0. */ + MPFR_SET_ZERO(z); + MPFR_RET(0); + } + else + { + /* Return 1. */ + return mpfr_set_ui (z, 1, GMP_RNDN); + } + } + } + MPFR_ASSERTD(MPFR_IS_FP(y)); if (MPFR_IS_INF(x)) - { - MPFR_CLEAR_FLAGS(z); - if (MPFR_SIGN(y) > 0) - MPFR_SET_INF(z); - else - MPFR_SET_ZERO(z); - MPFR_SET_POS(z); - MPFR_RET(0); - } - else - { - mpfr_t one; - int cmp; - - mpfr_init2(one, BITS_PER_MP_LIMB); - mpfr_set_ui(one, 1, GMP_RNDN); - cmp = mpfr_cmpabs(x, one) * MPFR_SIGN(y); - mpfr_clear(one); - - if (cmp > 0) - { - /* Return +inf. */ - MPFR_CLEAR_NAN(z); - MPFR_SET_INF(z); - MPFR_SET_POS(z); - MPFR_RET(0); - } - else if (cmp < 0) - { - /* Return +0. */ - MPFR_CLEAR_FLAGS(z); - MPFR_SET_ZERO(z); - MPFR_SET_POS(z); - MPFR_RET(0); - } - else - { - /* Return 1. */ - return mpfr_set_ui (z, 1, GMP_RNDN); - } - } - } - - MPFR_ASSERTN(MPFR_IS_FP(y)); - - if (MPFR_IS_INF(x)) - { - int negative; - /* Determine the sign now, in case y and z are the same object */ - negative = MPFR_SIGN(x) < 0 && is_odd(y); - MPFR_CLEAR_FLAGS(z); - if (MPFR_SIGN(y) > 0) - MPFR_SET_INF(z); - else - MPFR_SET_ZERO(z); - if (negative) - MPFR_SET_NEG(z); - else - MPFR_SET_POS(z); - MPFR_RET(0); + { + int negative; + /* Determine the sign now, in case y and z are the same object */ + negative = MPFR_IS_NEG(x) && is_odd(y); + MPFR_CLEAR_FLAGS(z); + if (MPFR_IS_POS(y)) + MPFR_SET_INF(z); + else + MPFR_SET_ZERO(z); + if (negative) + MPFR_SET_NEG(z); + else + MPFR_SET_POS(z); + MPFR_RET(0); + } + MPFR_ASSERTD(MPFR_IS_FP(x)); + if (MPFR_IS_ZERO(x)) + { + int negative; + /* Determine the sign now, in case y and z are the same object */ + negative = MPFR_IS_NEG(x) && is_odd (y); + MPFR_CLEAR_FLAGS(z); + if (MPFR_IS_NEG(y)) + MPFR_SET_INF(z); + else + MPFR_SET_ZERO(z); + if (negative) + MPFR_SET_NEG(z); + else + MPFR_SET_POS(z); + MPFR_RET(0); + } + /* Should never reach this code */ + MPFR_ASSERTN(1); + return 0; } - MPFR_ASSERTN(MPFR_IS_FP(x)); - - if (MPFR_IS_ZERO(x)) - { - int negative; - /* Determine the sign now, in case y and z are the same object */ - negative = MPFR_SIGN(x) < 0 && is_odd (y); - MPFR_CLEAR_FLAGS(z); - if (MPFR_SIGN(y) < 0) - MPFR_SET_INF(z); - else - MPFR_SET_ZERO(z); - if (negative) - MPFR_SET_NEG(z); - else - MPFR_SET_POS(z); - MPFR_RET(0); - } - else if (mpfr_cmp_ui (x, 1) == 0) /* 1^y is always 1 */ + if (mpfr_cmp_ui (x, 1) == 0) /* 1^y is always 1 */ { mpfr_set_ui (z, 1, GMP_RNDN); MPFR_RET(0); @@ -293,7 +289,7 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode) return mpfr_pow_si (z, x, zii, rnd_mode); } - if (MPFR_SIGN(x) < 0) + if (MPFR_IS_NEG(x)) { MPFR_SET_NAN(z); MPFR_RET_NAN; @@ -34,42 +34,39 @@ MA 02111-1307, USA. */ int mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) { - if (n > 0) + if (n >= 0) return mpfr_pow_ui (y, x, n, rnd_mode); else { - if (MPFR_IS_NAN(x)) - { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - if (n == 0) - return mpfr_set_ui (y, 1, GMP_RNDN); - - if (MPFR_IS_INF(x)) - { - MPFR_SET_ZERO(y); - if (MPFR_SIGN(x) > 0 || ((unsigned) n & 1) == 0) - MPFR_SET_POS(y); - else - MPFR_SET_NEG(y); - MPFR_RET(0); - } - - if (MPFR_IS_ZERO(x)) - { - MPFR_SET_INF(y); - if (MPFR_SIGN(x) > 0 || ((unsigned) n & 1) == 0) - MPFR_SET_POS(y); - else - MPFR_SET_NEG(y); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); + MPFR_CLEAR_FLAGS(y); + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) + { + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(x)) + { + MPFR_SET_ZERO(y); + if (MPFR_IS_POS(x) || ((unsigned) n & 1) == 0) + MPFR_SET_POS(y); + else + MPFR_SET_NEG(y); + MPFR_RET(0); + } + else if (MPFR_IS_ZERO(x)) + { + MPFR_SET_INF(y); + if (MPFR_IS_POS(x) || ((unsigned) n & 1) == 0) + MPFR_SET_POS(y); + else + MPFR_SET_NEG(y); + MPFR_RET(0); + } + else + MPFR_ASSERTN(1); + } /* detect exact powers: x^(-n) is exact iff x is a power of 2 */ if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0) @@ -132,3 +129,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode) } } } + + + + @@ -36,46 +36,53 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) int inexact; mp_rnd_t rnd1; - if (MPFR_IS_NAN(y)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(y) )) { - MPFR_SET_NAN(x); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(x); - - if (n == 0) /* y^0 = 1 for any y */ - { - /* The return mpfr_set_ui is important as 1 isn't necessarily - in the exponent range. */ - return mpfr_set_ui (x, 1, rnd); - } - - if (MPFR_IS_INF(y)) - { - /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ - if ((MPFR_SIGN(y) < 0) && (n % 2 == 1)) - MPFR_SET_NEG(x); + if (MPFR_IS_NAN(y)) + { + MPFR_SET_NAN(x); + MPFR_RET_NAN; + } + else if (n == 0) /* y^0 = 1 for any y */ + { + /* The return mpfr_set_ui is important as 1 isn't necessarily + in the exponent range. */ + return mpfr_set_ui (x, 1, rnd); + } + else if (MPFR_IS_INF(y)) + { + /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ + if ((MPFR_IS_NEG(y)) && ((n & 1) == 1)) + MPFR_SET_NEG(x); + else + MPFR_SET_POS(x); + MPFR_SET_INF(x); + MPFR_RET(0); + } + else if (MPFR_IS_ZERO(y)) /* 0^n = 0 for any n */ + { + MPFR_SET_ZERO(x); + MPFR_RET(0); + } else - MPFR_SET_POS(x); - MPFR_SET_INF(x); - MPFR_RET(0); + MPFR_ASSERTN(1); } - - MPFR_CLEAR_INF(x); - - if (MPFR_IS_ZERO(y)) /* 0^n = 0 for any n */ - { - MPFR_SET_ZERO(x); - MPFR_RET(0); + else if (MPFR_UNLIKELY( n <= 1)) + { + if (n == 0) + /* y^0 = 1 for any y */ + return mpfr_set_ui (x, 1, rnd); + MPFR_ASSERTD(n==1); + /* y^1 = y */ + return mpfr_set(x, y, rnd); } - + mpfr_save_emin_emax (); mpfr_init (res); prec = MPFR_PREC(x); - rnd1 = (MPFR_SIGN(y) > 0) ? GMP_RNDU : GMP_RNDD; /* away */ + rnd1 = (MPFR_IS_POS(y)) ? GMP_RNDU : GMP_RNDD; /* away */ do { @@ -96,7 +103,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) } /* check underflow */ - if (MPFR_EXP(res) <= (double) __gmpfr_emin) + if (MPFR_UNLIKELY(MPFR_EXP(res) <= (double) __gmpfr_emin)) { mpfr_clear (res); mpfr_restore_emin_emax (); diff --git a/round_prec.c b/round_prec.c index c86733123..9abe0f9c0 100644 --- a/round_prec.c +++ b/round_prec.c @@ -169,17 +169,20 @@ mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode) ow = MPFR_GET_ALLOC_SIZE(x); if (nw > ow) { + /* Realloc mantissa */ mp_ptr tmp = (mp_ptr) (*__gmp_reallocate_func) (MPFR_GET_REAL_PTR(x), MPFR_ALLOC_SIZE(ow), MPFR_ALLOC_SIZE(nw)); MPFR_SET_MANT_PTR(x, tmp); MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */ } - if (MPFR_IS_NAN(x)) - MPFR_RET_NAN; - - if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) - return 0; /* infinity and zero are exact */ + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) + { + if (MPFR_IS_NAN(x)) + MPFR_RET_NAN; + MPFR_ASSERTD(MPFR_IS_INF(x) || MPFR_IS_ZERO(x)); + return 0; /* infinity and zero are exact */ + } /* x is a non-zero real number */ @@ -190,11 +193,11 @@ mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode) prec, rnd_mode, &inexact); MPFR_PREC(x) = prec; - if (carry) + if (MPFR_UNLIKELY(carry)) { mp_exp_t exp = MPFR_EXP (x); - if (exp == __gmpfr_emax) + if (MPFR_UNLIKELY(exp == __gmpfr_emax)) (void) mpfr_set_overflow(x, rnd_mode, MPFR_SIGN(x)); else { @@ -225,12 +228,11 @@ int mpfr_can_round (mpfr_ptr b, mp_exp_t err, mp_rnd_t rnd1, mp_rnd_t rnd2, mp_prec_t prec) { - return MPFR_IS_FP(b) ? - (MPFR_IS_ZERO(b) ? 0 : /* we cannot round if b=0 */ - mpfr_can_round_raw (MPFR_MANT(b), - (MPFR_PREC(b) - 1)/BITS_PER_MP_LIMB + 1, - MPFR_SIGN(b), err, rnd1, rnd2, prec)) - : 0; /* cannnot round NaN or Inf */ + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b))) + return 0; /* We cannot round if Zero, Nan or Inf */ + else + return mpfr_can_round_raw(MPFR_MANT(b), MPFR_LIMB_SIZE(b), + MPFR_SIGN(b), err, rnd1, rnd2, prec); } int @@ -244,6 +246,7 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, mp_limb_t *tmp; TMP_DECL(marker); + /* FIXME: likely or not ? To check... How?*/ if (err0 < 0 || (mp_exp_unsigned_t) err0 <= prec) return 0; /* can't round */ @@ -253,7 +256,7 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, /* can round only in rounding to the nearest and err0 >= prec + 2 */ } - neg = neg <= 0; + neg = MPFR_IS_NEG_SIGN(neg); /* if the error is smaller than ulp(b), then anyway it will propagate up to ulp(b) */ @@ -281,15 +284,6 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not change bp[bn-1] >> s1, then we can round */ - if (rnd1 == GMP_RNDU) - if (neg) - rnd1 = GMP_RNDZ; - - if (rnd1 == GMP_RNDD) - rnd1 = neg ? GMP_RNDU : GMP_RNDZ; - - /* in the sequel, RNDU = towards infinity, RNDZ = towards zero */ - TMP_MARK(marker); tn = bn; k++; /* since we work with k+1 everywhere */ @@ -297,27 +291,40 @@ mpfr_can_round_raw (mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0, if (bn > k) MPN_COPY (tmp, bp, bn - k); - if (rnd1 != GMP_RNDN) - { /* GMP_RNDZ or GMP_RNDU */ + MPFR_ASSERTD (k > 0); + + switch (rnd1) + { + case GMP_RNDD: + if (neg) + goto round_rndu; + else + goto round_rndz; + break; + case GMP_RNDU: + if (neg) + goto round_rndz; + round_rndu: cc = (bp[bn - 1] >> s1) & 1; cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec); - /* now round b +/- 2^(MPFR_EXP(b)-err) */ - MPFR_ASSERTN (k > 0); - cc2 = rnd1 == GMP_RNDZ ? - mpn_add_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s) : - mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s); - } - else - { /* GMP_RNDN */ - /* first round b+2^(MPFR_EXP(b)-err) */ - MPFR_ASSERTN (k > 0); + cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s); + break; + case GMP_RNDZ: + round_rndz: + cc = (bp[bn - 1] >> s1) & 1; + cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec); + /* now round b +/- 2^(MPFR_EXP(b)-err) */ + cc2 = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s); + break; + case GMP_RNDN: + /* first round b+2^(MPFR_EXP(b)-err) */ cc = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s); cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */ cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec); - /* now round b-2^(MPFR_EXP(b)-err) */ cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MP_LIMB_T_ONE << s); + break; } if (cc2 && cc) @@ -29,53 +29,56 @@ int mpfr_set4 (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode, int signb) { int inex; - - if (MPFR_IS_NAN(b)) + + MPFR_CLEAR_FLAGS(a); + + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(b) )) { - MPFR_CLEAR_FLAGS(a); - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - - if (MPFR_IS_INF(b)) - { - MPFR_CLEAR_FLAGS(a); - MPFR_SET_INF(a); - inex = 0; - } - else - { - MPFR_CLEAR_FLAGS(a); - if (MPFR_IS_ZERO(b)) - { + if (MPFR_IS_NAN(b)) + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(b)) + { + MPFR_CLEAR_FLAGS(a); + MPFR_SET_INF(a); + inex = 0; + } + else if (MPFR_IS_ZERO(b)) + { MPFR_SET_ZERO(a); inex = 0; } else - { - mp_limb_t *ap; - mp_prec_t aq; - int carry; - - ap = MPFR_MANT(a); - aq = MPFR_PREC(a); - - carry = mpfr_round_raw(ap, MPFR_MANT(b), MPFR_PREC(b), (signb < 0), - aq, rnd_mode, &inex); - - if (carry) - { - mp_exp_t exp = MPFR_GET_EXP (b); - - if (exp == __gmpfr_emax) - return mpfr_set_overflow(a, rnd_mode, signb); - - MPFR_SET_EXP(a, exp + 1); - ap[(MPFR_PREC(a)-1)/BITS_PER_MP_LIMB] = MPFR_LIMB_HIGHBIT; - } - else - MPFR_SET_EXP (a, MPFR_GET_EXP (b)); - } + /* Should never reach this code */ + MPFR_ASSERTN(1); + } + else + { + mp_limb_t *ap; + mp_prec_t aq; + int carry; + + ap = MPFR_MANT(a); + aq = MPFR_PREC(a); + + carry = mpfr_round_raw(ap, MPFR_MANT(b), MPFR_PREC(b), + MPFR_IS_NEG_SIGN(signb), + aq, rnd_mode, &inex); + /* FIXME: Carry is likely or not ? */ + if (MPFR_UNLIKELY(carry)) + { + mp_exp_t exp = MPFR_GET_EXP (b); + + if (MPFR_UNLIKELY(exp == __gmpfr_emax)) + return mpfr_set_overflow(a, rnd_mode, signb); + + MPFR_SET_EXP(a, exp + 1); + ap[(MPFR_PREC(a)-1)/BITS_PER_MP_LIMB] = MPFR_LIMB_HIGHBIT; + } + else + MPFR_SET_EXP (a, MPFR_GET_EXP (b)); } MPFR_SET_SIGN(a, signb); @@ -51,9 +51,11 @@ mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode) /* don't forget to put zero in lower limbs */ MPN_ZERO(xp, xn); /* set sign */ - if ((i < 0) ^ (MPFR_SIGN(x) < 0)) - MPFR_CHANGE_SIGN(x); - + if (i < 0) + MPFR_SET_NEG(x); + else + MPFR_SET_POS(x); + nbits = BITS_PER_MP_LIMB - cnt; MPFR_EXP (x) = nbits; /* may be out-of-range, check range below */ inex = mpfr_check_range(x, 0, rnd_mode); @@ -61,18 +63,17 @@ mpfr_set_si (mpfr_ptr x, long i, mp_rnd_t rnd_mode) return inex; /* underflow or overflow */ /* round if MPFR_PREC(x) smaller than length of i */ - if (MPFR_PREC(x) < nbits) + if (MPFR_UNLIKELY(MPFR_PREC(x) < nbits)) { int carry; carry = mpfr_round_raw(xp+xn, xp+xn, nbits, (i < 0), MPFR_PREC(x), rnd_mode, &inex); - if (carry) + if (MPFR_UNLIKELY(carry)) { /* nbits is the current exponent */ - if (nbits == __gmpfr_emax) + if (MPFR_UNLIKELY(nbits == __gmpfr_emax)) return mpfr_set_overflow(x, rnd_mode, (i < 0 ? -1 : 1)); - MPFR_SET_EXP (x, nbits + 1); xp[xn] = MPFR_LIMB_HIGHBIT; } @@ -40,7 +40,7 @@ mpfr_set_ui (mpfr_ptr x, unsigned long i, mp_rnd_t rnd_mode) mp_limb_t *xp; xn = (MPFR_PREC(x)-1)/BITS_PER_MP_LIMB; - MPFR_ASSERTN(i == (mp_limb_t) i); + MPFR_ASSERTD(i == (mp_limb_t) i); count_leading_zeros(cnt, (mp_limb_t) i); xp = MPFR_MANT(x); @@ -55,15 +55,15 @@ mpfr_set_ui (mpfr_ptr x, unsigned long i, mp_rnd_t rnd_mode) return inex; /* underflow or overflow */ /* round if MPFR_PREC(x) smaller than length of i */ - if (MPFR_PREC(x) < nbits) + if (MPFR_UNLIKELY( MPFR_PREC(x) < nbits) ) { int carry; carry = mpfr_round_raw(xp+xn, xp+xn, nbits, 0, MPFR_PREC(x), rnd_mode, &inex); - if (carry) + if (MPFR_UNLIKELY(carry)) { /* nbits is the current exponent */ - if (nbits == __gmpfr_emax) + if (MPFR_UNLIKELY(nbits == __gmpfr_emax)) return mpfr_set_overflow(x, rnd_mode, 1); MPFR_SET_EXP (x, nbits + 1); @@ -90,18 +90,21 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int precy, m, ok, e, inexact, sign; mpfr_t c; - if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - if (MPFR_IS_ZERO(x)) - { - MPFR_CLEAR_FLAGS(y); - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); + if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(x)) + { + MPFR_CLEAR_FLAGS(y); + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + MPFR_ASSERTN(1); } precy = MPFR_PREC(y); @@ -119,8 +122,8 @@ mpfr_sin (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) mpfr_ui_sub (c, 1, c, GMP_RNDN); e = 2 + (- MPFR_GET_EXP (c)) / 2; mpfr_sqrt (c, c, GMP_RNDN); - if (sign < 0) - mpfr_neg (c, c, GMP_RNDN); + if (MPFR_IS_NEG_SIGN(sign)) + MPFR_CHANGE_SIGN(c); /* the absolute error on c is at most 2^(e-m) = 2^(EXP(c)-err) */ e = MPFR_GET_EXP (c) + m - e; @@ -38,33 +38,33 @@ mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode) mp_prec_t Nxt = MPFR_PREC(xt); int flag_neg=0, inexact =0; - if (MPFR_IS_NAN(xt)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(xt)) - { - MPFR_SET_INF(y); - MPFR_SET_SAME_SIGN(y, xt); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(y); - - if (MPFR_IS_ZERO(xt)) - { - MPFR_SET_ZERO(y); /* sinh(0) = 0 */ - MPFR_SET_SAME_SIGN(y, xt); - MPFR_RET(0); + if (MPFR_IS_NAN(xt)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(xt)) + { + MPFR_SET_INF(y); + MPFR_SET_SAME_SIGN(y, xt); + MPFR_RET(0); + } + else if (MPFR_IS_ZERO(xt)) + { + MPFR_SET_ZERO(y); /* sinh(0) = 0 */ + MPFR_SET_SAME_SIGN(y, xt); + MPFR_RET(0); + } + /* Should never reach this */ + MPFR_ASSERTN(1); } mpfr_init2 (x, Nxt); mpfr_set (x, xt, GMP_RNDN); - if(MPFR_SIGN(x)<0) + if (MPFR_IS_NEG(x)) { MPFR_CHANGE_SIGN(x); flag_neg=1; @@ -42,48 +42,44 @@ mpfr_sqrt (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) TMP_DECL(marker); - if (MPFR_IS_NAN(u)) - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - - if (MPFR_SIGN(u) < 0) - { - if (MPFR_IS_INF(u) || MPFR_NOTZERO(u)) - { - MPFR_SET_NAN(r); - MPFR_RET_NAN; - } - else - { /* sqrt(-0) = -0 */ - MPFR_CLEAR_FLAGS(r); - MPFR_SET_ZERO(r); - MPFR_SET_NEG(r); - MPFR_RET(0); - } - } - - MPFR_CLEAR_NAN(r); - MPFR_SET_POS(r); - - if (MPFR_IS_INF(u)) - { - MPFR_SET_INF(r); - MPFR_RET(0); - } - - MPFR_CLEAR_INF(r); - - if (MPFR_IS_ZERO(u)) - { - MPFR_SET_ZERO(r); - MPFR_RET(0); /* zero is exact */ - } - - up = MPFR_MANT(u); - usize = (MPFR_PREC(u) - 1)/BITS_PER_MP_LIMB + 1; - + MPFR_CLEAR_FLAGS(u); + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(u))) + { + if (MPFR_IS_NAN(u)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(u)) + { + /* 0+ or 0- */ + MPFR_SET_SAME_SIGN(r, u); + MPFR_SET_ZERO(r); + MPFR_RET(0); /* zero is exact */ + } + else if (MPFR_IS_INF(u)) + { + /* sqrt(-Inf) = NAN */ + if (MPFR_IS_NEG(u)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + MPFR_SET_POS(r); + MPFR_SET_INF(r); + MPFR_RET(0); + } + } + if (MPFR_UNLIKELY(MPFR_IS_NEG(u))) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + MPFR_SET_POS(r); + + up = MPFR_MANT(u); + usize = MPFR_LIMB_SIZE(u); + #ifdef DEBUG printf("Entering square root : "); for(k = usize - 1; k >= 0; k--) { printf("%lu ", up[k]); } @@ -28,61 +28,60 @@ MA 02111-1307, USA. */ int mpfr_sub (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t rnd_mode) { - if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + if (MPFR_ARE_SINGULAR(b,c)) { - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(a); - - if (MPFR_IS_INF(b)) - { - if (!MPFR_IS_INF(c) || MPFR_SIGN(b) != MPFR_SIGN(c)) - { - MPFR_SET_INF(a); - MPFR_SET_SAME_SIGN(a, b); - MPFR_RET(0); /* exact */ - } + if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + MPFR_CLEAR_NAN(a); + if (MPFR_IS_INF(b)) + { + if (!MPFR_IS_INF(c) || MPFR_SIGN(b) != MPFR_SIGN(c)) + { + MPFR_SET_INF(a); + MPFR_SET_SAME_SIGN(a, b); + MPFR_RET(0); /* exact */ + } + else + { + MPFR_SET_NAN(a); + MPFR_RET_NAN; + } + } else - { - MPFR_SET_NAN(a); - MPFR_RET_NAN; - } - } - else - if (MPFR_IS_INF(c)) - { - MPFR_SET_INF(a); - if (MPFR_SIGN(c) == MPFR_SIGN(a)) - MPFR_CHANGE_SIGN(a); - MPFR_RET(0); /* exact */ - } - - MPFR_ASSERTN(MPFR_IS_FP(b) && MPFR_IS_FP(c)); - - if (MPFR_IS_ZERO(b)) - { + if (MPFR_IS_INF(c)) + { + MPFR_SET_INF(a); + if (MPFR_SIGN(c) == MPFR_SIGN(a)) + MPFR_CHANGE_SIGN(a); + MPFR_RET(0); /* exact */ + } + if (MPFR_IS_ZERO(b)) + { + if (MPFR_IS_ZERO(c)) + { + if (MPFR_SIGN(a) != + (rnd_mode != GMP_RNDD ? + ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) > 0) ? -1 : 1) : + ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) < 0) ? 1 : -1))) + MPFR_CHANGE_SIGN(a); + MPFR_CLEAR_INF(a); + MPFR_SET_ZERO(a); + MPFR_RET(0); /* 0 - 0 is exact */ + } + return mpfr_neg (a, c, rnd_mode); + } if (MPFR_IS_ZERO(c)) - { - if (MPFR_SIGN(a) != - (rnd_mode != GMP_RNDD ? - ((MPFR_SIGN(b) < 0 && MPFR_SIGN(c) > 0) ? -1 : 1) : - ((MPFR_SIGN(b) > 0 && MPFR_SIGN(c) < 0) ? 1 : -1))) - MPFR_CHANGE_SIGN(a); - MPFR_CLEAR_INF(a); - MPFR_SET_ZERO(a); - MPFR_RET(0); /* 0 - 0 is exact */ - } - return mpfr_neg (a, c, rnd_mode); + { + return mpfr_set (a, b, rnd_mode); + } + /* Should never reach here */ + MPFR_ASSERTN(1); } - - if (MPFR_IS_ZERO(c)) - { - return mpfr_set (a, b, rnd_mode); - } - - MPFR_CLEAR_INF(a); + MPFR_CLEAR_FLAGS(a); + MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c)); if (MPFR_SIGN(b) == MPFR_SIGN(c)) { /* signs are equal, it's a real subtraction */ diff --git a/sub_one_ulp.c b/sub_one_ulp.c index 361e5853e..55d80ed5f 100644 --- a/sub_one_ulp.c +++ b/sub_one_ulp.c @@ -31,12 +31,14 @@ mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode) mp_size_t xn; int sh; mp_limb_t *xp; - - if (MPFR_IS_NAN(x)) - MPFR_RET_NAN; - - if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) - return 0; + + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) )) + { + if (MPFR_IS_NAN(x)) + MPFR_RET_NAN; + if (MPFR_IS_INF(x) || MPFR_IS_ZERO(x)) + return 0; + } MPFR_ASSERTN(MPFR_PREC_MIN > 1); @@ -49,7 +51,7 @@ mpfr_sub_one_ulp(mpfr_ptr x, mp_rnd_t rnd_mode) mp_exp_t exp = MPFR_EXP (x); /* Note: In case of underflow and rounding to the nearest mode, x won't be changed. Beware of infinite loops! */ - if (exp == __gmpfr_emin) + if (MPFR_UNLIKELY( exp == __gmpfr_emin )) return mpfr_set_underflow(x, rnd_mode, MPFR_SIGN(x)); else { @@ -31,18 +31,22 @@ mpfr_tan (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode) int precy, m, ok, inexact; mpfr_t s, c; - if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - if (MPFR_IS_ZERO(x)) - { - MPFR_CLEAR_FLAGS(y); - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN(y, x); - MPFR_RET(0); + if (MPFR_IS_NAN(x) || MPFR_IS_INF(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_ZERO(x)) + { + MPFR_CLEAR_FLAGS(y); + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y, x); + MPFR_RET(0); + } + /* Should never reach this point */ + MPFR_ASSERTN(1); } precy = MPFR_PREC(y); @@ -30,14 +30,7 @@ MA 02111-1307, USA. */ */ int -#if __STDC__ mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode) -#else -mpfr_tanh (y, xt, rnd_mode) - mpfr_ptr y; - mpfr_srcptr xt; - mp_rnd_t rnd_mode; -#endif { /****** Declaration ******/ @@ -46,35 +39,35 @@ mpfr_tanh (y, xt, rnd_mode) int flag_neg=0, inexact=0; /* Special value checking */ - - if (MPFR_IS_NAN(xt)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(xt)) - { - if (MPFR_SIGN(xt) > 0) - return mpfr_set_si(y,1,rnd_mode); /* tanh(inf) = 1 */ - else - return mpfr_set_si(y,-1,rnd_mode); /* tanh(-inf) = -1 */ - } - MPFR_CLEAR_INF(y); - - /* tanh(0) = 0 */ - if (MPFR_IS_ZERO(xt)) - { - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN(y,xt); - MPFR_RET(0); + if (MPFR_IS_NAN(xt)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(xt)) + { + if (MPFR_IS_POS(xt)) + return mpfr_set_si(y,1,rnd_mode); /* tanh(inf) = 1 */ + else + return mpfr_set_si(y,-1,rnd_mode); /* tanh(-inf) = -1 */ + } + /* tanh(0) = 0 */ + else if (MPFR_IS_ZERO(xt)) + { + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y,xt); + MPFR_RET(0); + } + /* Should never reach this point */ + MPFR_ASSERTN(1); } mpfr_init2(x,Nxt); mpfr_set(x,xt,GMP_RNDN); - if (MPFR_SIGN(x) < 0) + if (MPFR_IS_NEG(x)) { MPFR_CHANGE_SIGN(x); flag_neg=1; diff --git a/tests/tpow.c b/tests/tpow.c index fa183e77e..2161d80e7 100644 --- a/tests/tpow.c +++ b/tests/tpow.c @@ -206,10 +206,14 @@ special () mpfr_clear (t); } +void f(void); + static void particular_cases (void) { mpfr_t t[11], r; + static const char *name[11] = { + "NaN", "+inf", "-inf", "+0", "-0", "+1", "-1", "+2", "-2", "+0.5", "-0.5"}; int i, j; int error = 0; @@ -244,6 +248,8 @@ particular_cases (void) /* +0.5 */ { 0, 2, 1, 128, 128, 64, 256, 32, 512, 90, 180 }, /* -0.5 */ { 0, 2, 1, 128, 128, -64,-256, 32, 512, 0, 0 } }; + if (i == 5 && j == 1) + f(); mpfr_pow (r, t[i], t[j], GMP_RNDN); p = mpfr_nan_p (r) ? 0 : mpfr_inf_p (r) ? 1 : @@ -253,8 +259,8 @@ particular_cases (void) p = -p; if (p != q[i][j]) { - printf ("Error in mpfr_pow for particular case (%d,%d):\n" - "got %d instead of %d\n", i, j, p, q[i][j]); + printf ("Error in mpfr_pow for particular case (%s)^(%s) (%d,%d):\n" + "got %d instead of %d\n", name[i], name[j], i,j,p, q[i][j]); error = 1; } } @@ -267,6 +273,10 @@ particular_cases (void) exit (1); } +void f(void) +{ +} + static void underflows (void) { @@ -33,25 +33,39 @@ mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) mp_limb_t up[1]; unsigned long cnt; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) /* u/Inf = 0 */ - { - MPFR_CLEAR_INF(y); - MPFR_SET_ZERO(y); - MPFR_SET_SAME_SIGN(y,x); - MPFR_RET(0); + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF(x)) /* u/Inf = 0 */ + { + MPFR_SET_ZERO(y); + MPFR_SET_SAME_SIGN(y,x); + MPFR_RET(0); + } + else if (MPFR_IS_ZERO(x)) /* u / 0 */ + { + if (u) + { + /* u > 0, so y = sign(x) * Inf */ + MPFR_SET_SAME_SIGN(y, x); + MPFR_SET_INF(y); + MPFR_RET(0); + } + else + { + /* 0 / 0 */ + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + } + MPFR_ASSERTN(1); + return 0; /* To avoid warning */ } - - MPFR_CLEAR_INF(y); - - if (u) + else if (u) { MPFR_TMP_INIT1(up, uu, BITS_PER_MP_LIMB); MPFR_ASSERTN(u == (mp_limb_t) u); @@ -60,17 +74,9 @@ mpfr_ui_div (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) MPFR_SET_EXP (uu, BITS_PER_MP_LIMB - cnt); return mpfr_div (y, uu, x, rnd_mode); } - else /* u = 0 */ + else /* u = 0, and x != 0 */ { - if (MPFR_IS_ZERO(x)) /* 0/0 */ - { - MPFR_SET_NAN(y); - MPFR_RET_NAN; - } - else - { - MPFR_SET_ZERO(y); /* if u=0, then set y to 0 */ - MPFR_RET(0); - } + MPFR_SET_ZERO(y); /* if u=0, then set y to 0 */ + MPFR_RET(0); } } @@ -32,23 +32,28 @@ mpfr_ui_sub (mpfr_ptr y, unsigned long int u, mpfr_srcptr x, mp_rnd_t rnd_mode) mp_limb_t up[1]; unsigned long cnt; - if (MPFR_IS_NAN(x)) + if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(x))) { - MPFR_SET_NAN(y); - MPFR_RET_NAN; + if (MPFR_IS_NAN(x)) + { + MPFR_SET_NAN(y); + MPFR_RET_NAN; + } + if (MPFR_IS_INF(x)) + { + /* u - Inf = -Inf and u - -Inf = +Inf */ + MPFR_SET_INF(y); + MPFR_SET_OPPOSITE_SIGN(y,x); + MPFR_RET(0); /* +/-infinity is exact */ + } + if (MPFR_IS_ZERO(x)) + /* u - 0 = u */ + return mpfr_set_ui(y, u, rnd_mode); + /* Should never reach this code */ + MPFR_ASSERTN(1); + return 0; /* To avoid a warning. It isn't reached */ } - - MPFR_CLEAR_NAN(y); - - if (MPFR_IS_INF(x)) - { - MPFR_SET_INF(y); - if (MPFR_SIGN(x) == MPFR_SIGN(y)) - MPFR_CHANGE_SIGN(y); - MPFR_RET(0); /* +/-infinity is exact */ - } - - if (u) + else if (u) { MPFR_TMP_INIT1 (up, uu, BITS_PER_MP_LIMB); MPFR_ASSERTN(u == (mp_limb_t) u); @@ -332,32 +332,35 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode) mp_prec_t precz, prec1, precs, precs1; int inex; - if (mpfr_nan_p (s)) + /* Zero, Nan or Inf ? */ + if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(s) )) { - MPFR_SET_NAN (z); - MPFR_RET_NAN; - } - - if (mpfr_inf_p (s)) - { - if (MPFR_SIGN(s) > 0) - return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */ - MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ - MPFR_RET_NAN; - } - - /* now s is neither NaN nor Infinite */ - - if (mpfr_cmp_ui(s,0) == 0) /* Case s = 0 */ - { - mpfr_set_ui (z, 1, rnd_mode); - mpfr_div_2ui (z, z, 1, rnd_mode); - return mpfr_neg (z, z, rnd_mode); + if (MPFR_IS_NAN(s)) + { + MPFR_SET_NAN (z); + MPFR_RET_NAN; + } + if (MPFR_IS_INF(s)) + { + if (MPFR_SIGN(s) > 0) + return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */ + MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ + MPFR_RET_NAN; + } + if (MPFR_IS_ZERO(s)) + { + mpfr_set_ui (z, 1, rnd_mode); + mpfr_div_2ui (z, z, 1, rnd_mode); + MPFR_CHANGE_SIGN(z); + MPFR_RET(0); + } + MPFR_ASSERTN(1); } + /* s is neither Nan, nor Inf, nor Zero */ mpfr_init2(s2, mpfr_get_prec(s)); mpfr_div_2ui(s2, s, 1, rnd_mode); - if ((MPFR_SIGN(s) == -1) && (mpfr_floor(s2, s2) == 0)) /* Case s = -2n */ + if (MPFR_IS_NEG(s) && mpfr_floor(s2, s2) == 0) /* Case s = -2n */ { mpfr_clear (s2); return mpfr_set_ui (z, 0, rnd_mode); |