[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

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>