[BACK]Return to cmp.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/cmp.c, Revision 1.1

1.1     ! maekawa     1: /* mpfr_cmp -- compare two floating-point numbers
        !             2:
        !             3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
        !             4:
        !             5: This file is part of the MPFR Library.
        !             6:
        !             7: The MPFR Library is free software; you can redistribute it and/or modify
        !             8: it under the terms of the GNU Library General Public License as published by
        !             9: the Free Software Foundation; either version 2 of the License, or (at your
        !            10: option) any later version.
        !            11:
        !            12: The MPFR Library is distributed in the hope that it will be useful, but
        !            13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
        !            14: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
        !            15: License for more details.
        !            16:
        !            17: You should have received a copy of the GNU Library General Public License
        !            18: along with the MPFR Library; see the file COPYING.LIB.  If not, write to
        !            19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
        !            20: MA 02111-1307, USA. */
        !            21:
        !            22: #include <stdio.h>
        !            23: #include "gmp.h"
        !            24: #include "gmp-impl.h"
        !            25: #include "longlong.h"
        !            26: #include "mpfr.h"
        !            27:
        !            28: /* returns 0 iff b = c
        !            29:             a positive value iff b > c
        !            30:             a negative value iff b < c
        !            31:
        !            32: More precisely, in case b and c are of same sign, the absolute value
        !            33: of the result is one plus the absolute difference between the exponents
        !            34: of b and c, i.e. one plus the number of bits shifts to align b and c
        !            35: (this value is useful in mpfr_sub).
        !            36:
        !            37: */
        !            38:
        !            39: /* #define DEBUG */
        !            40:
        !            41: /* compares b and sign(s)*c */
        !            42: int
        !            43: #if __STDC__
        !            44: mpfr_cmp3 ( mpfr_srcptr b, mpfr_srcptr c, long int s)
        !            45: #else
        !            46: mpfr_cmp3(b, c, s)
        !            47:      mpfr_srcptr b;
        !            48:      mpfr_srcptr c;
        !            49:      long int s;
        !            50: #endif
        !            51: {
        !            52:    long int diff_exp;
        !            53:    unsigned long bn, cn;
        !            54:    mp_limb_t *bp, *cp;
        !            55:
        !            56:    if (!NOTZERO(b) && !NOTZERO(c)) { return 0; }
        !            57:    s = s*SIGN(b)*SIGN(c);
        !            58:    if (s<0) return(SIGN(b));
        !            59:
        !            60:    /* now signs are equal */
        !            61:
        !            62:    diff_exp = EXP(b)-EXP(c);
        !            63:    s = (SIGN(b)>0) ? 1 : -1;
        !            64:
        !            65:    if (diff_exp>0) return(s*(1+diff_exp));
        !            66:    else if (diff_exp<0) return(s*(-1+diff_exp));
        !            67:    /* both signs and exponents are equal */
        !            68:
        !            69:    bn = (PREC(b)-1)/mp_bits_per_limb+1;
        !            70:    cn = (PREC(c)-1)/mp_bits_per_limb+1;
        !            71:    bp = MANT(b); cp = MANT(c);
        !            72:
        !            73:    while (bn && cn) {
        !            74:      if (bp[--bn] != cp[--cn])
        !            75:        return((bp[bn]>cp[cn]) ? s : -s);
        !            76:    }
        !            77:
        !            78:    if (bn) { while (bn) if (bp[--bn]) return(s); }
        !            79:    else if (cn) while (cn) if (cp[--cn]) return(-s);
        !            80:
        !            81:    return 0;
        !            82: }
        !            83:
        !            84: /* returns the number of cancelled bits when one subtracts abs(c) from abs(b).
        !            85:    Assumes b>=c, which implies EXP(b)>=EXP(c).
        !            86:    if b=c, returns prec(b).
        !            87: */
        !            88: int
        !            89: #if __STDC__
        !            90: mpfr_cmp2 ( mpfr_srcptr b, mpfr_srcptr c )
        !            91: #else
        !            92: mpfr_cmp2(b, c)
        !            93:      mpfr_srcptr b;
        !            94:      mpfr_srcptr c;
        !            95: #endif
        !            96: {
        !            97:   long int d, bn, cn, k, z;
        !            98:   mp_limb_t *bp, *cp, t, u=0, cc=0;
        !            99:
        !           100: #ifdef DEBUG
        !           101:   printf("b="); mpfr_print_raw(b); putchar('\n');
        !           102:   printf("c="); mpfr_print_raw(c); putchar('\n');
        !           103: #endif
        !           104:   if (NOTZERO(c)==0) return (NOTZERO(b)) ? 0 : PREC(b);
        !           105:   d = EXP(b)-EXP(c);
        !           106:   k = 0; /* result can be d or d+1 if d>1, or >= d otherwise */
        !           107:   /* k is the number of identical bits in the high part,
        !           108:      then z is the number of possibly cancelled bits */
        !           109: #ifdef DEBUG
        !           110:   if (d<0) { printf("assumption EXP(b)<EXP(c) violated\n"); exit(1); }
        !           111: #endif
        !           112:   bn = (PREC(b)-1)/mp_bits_per_limb;
        !           113:   cn = (PREC(c)-1)/mp_bits_per_limb;
        !           114:   bp = MANT(b); cp = MANT(c);
        !           115:   /* subtract c from b from most significant to less significant limbs,
        !           116:      and first determines first non zero limb difference */
        !           117:   if (d)
        !           118:     {
        !           119:       cc = bp[bn--];
        !           120:       if (d<mp_bits_per_limb)
        !           121:        cc -= cp[cn]>>d;
        !           122:     }
        !           123:   else { /* d=0 */
        !           124:     while (bn>=0 && cn>=0 && (cc=(bp[bn--]-cp[cn--]))==0) {
        !           125:       k+=mp_bits_per_limb;
        !           126:     }
        !           127:
        !           128:     if (cc==0) { /* either bn<0 or cn<0 */
        !           129:       while (bn>=0 && (cc=bp[bn--])==0) k+=mp_bits_per_limb;
        !           130:     }
        !           131:     /* now bn<0 or cc<>0 */
        !           132:     if (cc==0 && bn<0) return(PREC(b));
        !           133:   }
        !           134:
        !           135:   /* the first non-zero limb difference is cc, and the number
        !           136:      of cancelled bits in the upper limbs is k */
        !           137:
        !           138:   count_leading_zeros(u, cc);
        !           139:   k += u;
        !           140:
        !           141:   if (cc != (1<<(mp_bits_per_limb-u-1))) return k;
        !           142:
        !           143:   /* now cc is an exact power of two */
        !           144:   if (cc != 1)
        !           145:     /* We just need to compare the following limbs */
        !           146:     /* until two of them differ. The result is either k or k+1. */
        !           147:     {
        !           148:       /* First flush all the unmatched limbs of b ; they all have to
        !           149:         be 0 in order for the process to go on */
        !           150:       while (bn >= 0)
        !           151:        {
        !           152:          if (cn < 0) { return k; }
        !           153:          t = bp[bn--];
        !           154:          if (d < mp_bits_per_limb)
        !           155:            {
        !           156:              if (d)
        !           157:                {
        !           158:                  u = cp[cn--] << (mp_bits_per_limb - d);
        !           159:                  if (cn >= 0) u+=(cp[cn]>>d);
        !           160:                }
        !           161:              else u = cp[cn--];
        !           162:
        !           163:              if (t > u || (t == u && cn < 0)) return k;
        !           164:              if (t < u) return k+1;
        !           165:            }
        !           166:          else
        !           167:            if (t) return k; else d -= mp_bits_per_limb;
        !           168:        }
        !           169:
        !           170:       /* bn < 0; if some limb of c is nonzero, return k+1, otherwise return k*/
        !           171:
        !           172:       if (cn>=0 && (cp[cn--] << (mp_bits_per_limb - d))) { return k+1; }
        !           173:
        !           174:       while (cn >= 0)
        !           175:        if (cp[cn--]) return k+1;
        !           176:       return k;
        !           177:     }
        !           178:
        !           179:   /* cc = 1. Too bad. */
        !           180:   z = 0; /* number of possibly cancelled bits - 1 */
        !           181:   /* thus result is either k if low(b) >= low(c)
        !           182:                         or k+z+1 if low(b) < low(c) */
        !           183:   if (d > mp_bits_per_limb) { return k; }
        !           184:
        !           185:   while (bn >= 0)
        !           186:     {
        !           187:       if (cn < 0) { return k; }
        !           188:
        !           189:       if (d)
        !           190:         {
        !           191:           u = cp[cn--] << (mp_bits_per_limb - d);
        !           192:           if (cn >= 0) u+=(cp[cn]>>d);
        !           193:         }
        !           194:        else u = cp[cn--];
        !           195:
        !           196:       /* bp[bn--] > cp[cn--] : no borrow possible, k unchanged
        !           197:         bp[bn--] = cp[cn--] : need to consider next limbs
        !           198:          bp[bn--] < cp[cn--] : borrow
        !           199:       */
        !           200:        if ((cc = bp[bn--]) != u) {
        !           201:         if (cc>u) return k;
        !           202:         else { count_leading_zeros(u, cc-u); return k + 1 + z + u; }
        !           203:        }
        !           204:        else z += mp_bits_per_limb;
        !           205:     }
        !           206:
        !           207:   if (cn >= 0)
        !           208:     count_leading_zeros(cc, ~(cp[cn--] << (mp_bits_per_limb - d)));
        !           209:   else { cc = 0; }
        !           210:
        !           211:   k += cc;
        !           212:   if (cc < d) return k;
        !           213:
        !           214:   while (cn >= 0 && !~cp[cn--]) { z += mp_bits_per_limb; }
        !           215:   if (cn >= 0) { count_leading_zeros(cc, ~cp[cn--]); return (k + z + cc); }
        !           216:
        !           217:   return k; /* We **need** that the nonsignificant limbs of c are set
        !           218:               to zero there */
        !           219: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>