diff options
Diffstat (limited to 'cmp2.c')
-rw-r--r-- | cmp2.c | 56 |
1 files changed, 39 insertions, 17 deletions
@@ -19,7 +19,6 @@ along with the MPFR Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ -#include <stdio.h> #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" @@ -27,8 +26,8 @@ MA 02111-1307, USA. */ #include "mpfr-impl.h" /* returns the number of cancelled bits when one subtracts abs(c) from abs(b). - Assumes b>=c, which implies MPFR_EXP(b)>=MPFR_EXP(c). - if b=c, returns prec(b). + Assumes |b| >= |c|, which implies MPFR_EXP(b)>=MPFR_EXP(c). + if |b| = |c|, returns prec(b). Assumes neither of b or c is NaN or +/- infinity. @@ -38,10 +37,17 @@ mp_prec_t mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) { mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0; - long bn, cn, z; - unsigned long diff_exp, res = 0; + mp_size_t bn, cn; + mp_exp_unsigned_t diff_exp; + mp_prec_t res = 0; - if (MPFR_IS_ZERO(c)) return 0; + MPFR_ASSERTN(MPFR_IS_FP(b)); + MPFR_ASSERTN(MPFR_IS_FP(c)); + + if (MPFR_IS_ZERO(c)) + return 0; + + MPFR_ASSERTN(MPFR_NOTZERO(b)); bp = MPFR_MANT(b); cp = MPFR_MANT(c); @@ -49,7 +55,8 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB; cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB; - diff_exp = MPFR_EXP(b) - MPFR_EXP(c); + MPFR_ASSERTN(MPFR_EXP(b) >= MPFR_EXP(c)); + diff_exp = (mp_exp_unsigned_t) MPFR_EXP(b) - MPFR_EXP(c); if (diff_exp == 0) /* otherwise the shifted most significant limb of c cannot match bp[bn] */ @@ -66,6 +73,8 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) if (cn < 0) /* c discards exactly the upper part of b */ { + unsigned int z; + while (bn>=0 && bp[bn]==0) { bn--; @@ -82,11 +91,13 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) /* now we have removed the identical upper limbs of b and c (can happen only when diff_exp = 0): bp[bn] > cc, bn>=0, cn>=0 */ + if (diff_exp < BITS_PER_MP_LIMB) { cc = cp[cn] >> diff_exp; /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */ - if (diff_exp) lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp); + if (diff_exp) + lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp); cn--; } else @@ -94,14 +105,21 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) dif = bp[bn--] - cc; /* necessarily dif >= 1 */ - while ((cn>=0 || lastc) && (high_dif==0) && (dif==1)) + while ((cn >= 0 || lastc != 0) && (high_dif == 0) && (dif == 1)) { /* dif=1 implies diff_exp = 0 or 1 */ bb = (bn >= 0) ? bp[bn] : 0; cc = lastc; if (cn >= 0) { - cc += cp[cn] >> diff_exp; - if (diff_exp) lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp); + if (diff_exp == 0) + { + cc += cp[cn]; + } + else /* diff_exp = 1 */ + { + cc += cp[cn] >> 1; + lastc = cp[cn] << (BITS_PER_MP_LIMB - 1); + } } else lastc = 0; @@ -113,14 +131,16 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */ - if (high_dif) /* necessarily high_dif = 1 */ + if (high_dif != 0) /* necessarily high_dif = 1 */ { res--; - if (dif) + if (dif != 0) return res; } else /* high_dif = 0 */ { + unsigned int z; + count_leading_zeros(z, dif); /* dif > 1 here */ res += z; if (dif != ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - z - 1))) @@ -128,7 +148,7 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) } /* now result is res + (low(b) < low(c)) */ - while (bn>=0 && (cn>=0 || lastc)) + while (bn >= 0 && (cn >= 0 || lastc != 0)) { if (diff_exp >= BITS_PER_MP_LIMB) diff_exp -= BITS_PER_MP_LIMB; @@ -138,7 +158,8 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) if (cn >= 0) { cc += cp[cn] >> diff_exp; - if (diff_exp) lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp); + if (diff_exp != 0) + lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp); } else lastc = 0; @@ -151,9 +172,10 @@ mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c) if (bn < 0) { - if (lastc) + if (lastc != 0) return res + 1; - while (cn>=0 && cp[cn]==0) cn--; + while (cn >= 0 && cp[cn] == 0) + cn--; return res + (cn >= 0); } |