diff options
Diffstat (limited to 'mpfr/cmp.c')
| -rw-r--r-- | mpfr/cmp.c | 255 |
1 files changed, 50 insertions, 205 deletions
diff --git a/mpfr/cmp.c b/mpfr/cmp.c index 5b7b73cf5..ec6fdd166 100644 --- a/mpfr/cmp.c +++ b/mpfr/cmp.c @@ -1,249 +1,94 @@ /* mpfr_cmp -- compare two floating-point numbers -Copyright (C) 1999 Free Software Foundation. +Copyright (C) 1999, 2001 Free Software Foundation. This file is part of the MPFR Library. The MPFR Library is free software; you can redistribute it and/or modify -it under the terms of the GNU Library General Public License as published by -the Free Software Foundation; either version 2 of the License, or (at your +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. The 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 Library General Public +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 Library General Public License +You should have received a copy of the GNU Lesser General Public License 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" #include "mpfr.h" #include "mpfr-impl.h" -/* returns 0 iff b = c - a positive value iff b > c - a negative value iff b < c - -More precisely, in case b and c are of same sign, the absolute value -of the result is one plus the absolute difference between the exponents -of b and c, i.e. one plus the number of bits shifts to align b and c -(this value is useful in mpfr_sub). - +/* returns 0 iff b = sign(s) * c + a positive value iff b > sign(s) * c + a negative value iff b < sign(s) * c */ -/* #define DEBUG */ - -/* compares b and sign(s)*c */ -int -#if __STDC__ -mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, long int s) -#else -mpfr_cmp3 (b, c, s) - mpfr_srcptr b; - mpfr_srcptr c; - long int s; -#endif +int +mpfr_cmp3 (mpfr_srcptr b, mpfr_srcptr c, int s) { - long int diff_exp; - unsigned long bn, cn; - mp_limb_t *bp, *cp; - - if (MPFR_IS_NAN(b) || MPFR_IS_NAN(c)) return 1; + mp_exp_t be, ce; + mp_size_t bn, cn; + mp_limb_t *bp, *cp; - if (MPFR_IS_INF(b)) { - if (MPFR_IS_INF(c) && (MPFR_SIGN(b) * s * MPFR_SIGN(c) > 0)) - return 0; - else - return MPFR_SIGN(b); - } + MPFR_ASSERTN(!MPFR_IS_NAN(b)); + MPFR_ASSERTN(!MPFR_IS_NAN(c)); + s *= MPFR_SIGN(c); - if (!MPFR_NOTZERO(b)) { - if (!MPFR_NOTZERO(c)) return 0; else return -(s*MPFR_SIGN(c)); - } - else if (!MPFR_NOTZERO(c)) return MPFR_SIGN(b); + if (MPFR_IS_INF(b)) + return (MPFR_IS_INF(c) && s * MPFR_SIGN(b) > 0) ? 0 : MPFR_SIGN(b); - s = s * MPFR_SIGN(b) * MPFR_SIGN(c); - if (s<0) return(MPFR_SIGN(b)); + if (MPFR_IS_INF(c)) + return -s; - /* now signs are equal */ - diff_exp = MPFR_EXP(b)-MPFR_EXP(c); - s = (MPFR_SIGN(b) > 0) ? 1 : -1; + /* b and c are real numbers */ - if (diff_exp>0) return(s*(1+diff_exp)); - else if (diff_exp<0) return(s*(-1+diff_exp)); - /* both signs and exponents are equal */ + if (MPFR_IS_ZERO(b)) + return MPFR_IS_ZERO(c) ? 0 : -s; - bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB+1; - cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB+1; - bp = MPFR_MANT(b); cp = MPFR_MANT(c); + if (MPFR_IS_ZERO(c)) + return MPFR_SIGN(b); - while (bn && cn) { - if (bp[--bn] != cp[--cn]) - return((bp[bn]>cp[cn]) ? s : -s); - } + if (s * MPFR_SIGN(b) < 0) + return MPFR_SIGN(b); - if (bn) { while (bn) if (bp[--bn]) return(s); } - else if (cn) while (cn) if (cp[--cn]) return(-s); + /* now signs are equal */ - return 0; -} + be = MPFR_EXP(b); + ce = MPFR_EXP(c); + if (be > ce) + return s; + if (be < ce) + return -s; -/* 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). + /* both signs and exponents are equal */ - In other terms mpfr_cmp2 (b, c) returns EXP(b) - EXP(b-c). -*/ -int -#if __STDC__ -mpfr_cmp2 ( mpfr_srcptr b, mpfr_srcptr c ) -#else -mpfr_cmp2(b, c) - mpfr_srcptr b; - mpfr_srcptr c; -#endif -{ - long int d, bn, cn, k, z; - mp_limb_t *bp, *cp, t, u=0, cc=0; - -#ifdef DEBUG - printf("b="); mpfr_print_raw(b); putchar('\n'); - printf("c="); mpfr_print_raw(c); putchar('\n'); -#endif - if (MPFR_NOTZERO(c)==0) return (MPFR_NOTZERO(b)) ? 0 : MPFR_PREC(b); - d = MPFR_EXP(b)-MPFR_EXP(c); - k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */ - /* k is the number of identical bits in the high part, - then z is the number of possibly cancelled bits */ -#ifdef DEBUG - printf("d=%u\n", d); - if (d<0) { printf("assumption MPFR_EXP(b)<MPFR_EXP(c) violated\n"); exit(1); } -#endif bn = (MPFR_PREC(b)-1)/BITS_PER_MP_LIMB; cn = (MPFR_PREC(c)-1)/BITS_PER_MP_LIMB; - bp = MPFR_MANT(b); cp = MPFR_MANT(c); - /* subtract c from b from most significant to less significant limbs, - and first determines first non zero limb difference */ - if (d) - { - cc = bp[bn--]; - if (d < BITS_PER_MP_LIMB) - cc -= cp[cn] >> d; - } - else { /* d=0 */ - while (bn>=0 && cn>=0 && (cc=(bp[bn--]-cp[cn--]))==0) { - k+=BITS_PER_MP_LIMB; - } - - if (cc==0) { /* either bn<0 or cn<0 */ - while (bn>=0 && (cc=bp[bn--])==0) k+=BITS_PER_MP_LIMB; - } - /* now bn<0 or cc<>0 */ - if (cc==0 && bn<0) return(MPFR_PREC(b)); - } - - /* the first non-zero limb difference is cc, and the number - of cancelled bits in the upper limbs is k */ - count_leading_zeros(u, cc); - k += u; + bp = MPFR_MANT(b); + cp = MPFR_MANT(c); - if (cc != ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - u - 1))) return k; - - /* now cc is an exact power of two */ - if (cc != 1) - /* We just need to compare the following limbs */ - /* until two of them differ. The result is either k or k+1. */ + for ( ; bn >= 0 && cn >= 0; bn--, cn--) { - /* First flush all the unmatched limbs of b ; they all have to - be 0 in order for the process to go on */ - while (bn >= 0) - { - if (cn < 0) { return k; } - t = bp[bn--]; - if (d < BITS_PER_MP_LIMB) - { - if (d) - { - u = cp[cn--] << (BITS_PER_MP_LIMB - d); - if (cn >= 0) u+=(cp[cn]>>d); - } - else u = cp[cn--]; - - if (t > u || (t == u && cn < 0)) return k; - if (t < u) return k+1; - } - else - if (t) return k; else d -= BITS_PER_MP_LIMB; - } - - /* bn < 0; if some limb of c is nonzero, return k+1, otherwise return k*/ - - /* if d < BITS_PER_MP_LIMB, only the last d bits of cp[cn] have to be - considered, otherwise all bits */ - if (d < BITS_PER_MP_LIMB) - if (cn>=0 && (cp[cn--] << (BITS_PER_MP_LIMB - d))) { return k+1; } - - while (cn >= 0) - if (cp[cn--]) return k+1; - return k; + if (bp[bn] > cp[cn]) + return s; + if (bp[bn] < cp[cn]) + return -s; } - /* cc = 1. Too bad. */ - z = 0; /* number of possibly cancelled bits - 1 */ - /* thus result is either k if low(b) >= low(c) - or k+z+1 if low(b) < low(c) */ - if (d > BITS_PER_MP_LIMB) return k; - - while (bn >= 0) /* the next limb of b to be considered is b[bn] */ - { - /* for c we have to consider the low d bits of c[cn] - and the high (BITS_PER_MP_LIMB-d) bits of c[cn-1] */ - if (cn < 0) return k; - - if (d) - { - u = cp[cn--] << (BITS_PER_MP_LIMB - d); - if (cn >= 0) u += cp[cn] >> d; - } - else u = cp[cn--]; - - /* bp[bn--] > cp[cn--] : no borrow possible, k unchanged - bp[bn--] = cp[cn--] : need to consider next limbs - bp[bn--] < cp[cn--] : borrow - */ - if ((cc = bp[bn--]) != u) { - if (cc > u) return k; - else { - count_leading_zeros(u, cc-u); - z += u + 1; - if (u + 1 < BITS_PER_MP_LIMB) return k + z; - } - } - else z += BITS_PER_MP_LIMB; - } + for ( ; bn >= 0; bn--) + if (bp[bn]) + return s; - /* warning: count_leading_zeros doesn't work with zero */ - if ((cn >= 0) && d && (u = ~(cp[cn--] << (BITS_PER_MP_LIMB - d)))) - count_leading_zeros(cc, u); - else - cc = 0; - - /* here d=0 or d=1: if d=1, we have one more cancelled bit if we don't - shift cp[cn] */ - k += cc; - if (cc < d) return k; - - while (cn >= 0 && !~cp[cn]) { z += BITS_PER_MP_LIMB; cn--; } - /* now either cn<0 or cp[cn] is not 111...111 */ - if (cn >= 0) { count_leading_zeros(cc, ~cp[cn]); return (k + 1 + z + cc); } - - return k; /* We **need** that the nonsignificant limbs of c are set - to zero there */ + for ( ; cn >= 0; cn--) + if (cp[cn]) + return -s; + + return 0; } |