contrib/mpfr/cmp2.c

   1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23
  24 #define MPFR_NEED_LONGLONG_H
  25 #include "mpfr-impl.h"
  26
  27 /* If |b| != |c|, puts the number of canceled bits when one subtracts |c|
  28    from |b| in *cancel. Returns the sign of the difference.
  29
  30    Assumes neither of b or c is NaN, +/- infinity, or +/- 0.
  31
  32    In other terms, if |b| != |c|, mpfr_cmp2 (b, c) returns
  33    EXP(max(|b|,|c|)) - EXP(|b| - |c|).
  34 */
  35
  36 int
  37 mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mp_prec_t *cancel)
  38 {
  39   mp_limb_t *bp, *cp, bb, cc = 0, lastc = 0, dif, high_dif = 0;
  40   mp_size_t bn, cn;
  41   mp_exp_unsigned_t diff_exp;
  42   mp_prec_t res = 0;
  43   int sign;
  44
  45   /* b=c should not happen, since cmp2 is called only from agm
  46      (with different variables), and from sub1 (if same b=c, then
  47      sub1sp would be called instead */
  48   MPFR_ASSERTD (b != c);
  49
  50   /* the cases b=0 or c=0 are also treated apart in agm and sub
  51      (which calls sub1) */
  52   MPFR_ASSERTD (MPFR_IS_PURE_FP(b));
  53   MPFR_ASSERTD (MPFR_IS_PURE_FP(c));
  54
  55   if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c))
  56     {
  57       sign = 1;
  58       diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
  59
  60       bp = MPFR_MANT(b);
  61       cp = MPFR_MANT(c);
  62
  63       bn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
  64       cn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB; /* # of limbs of c minus 1 */
  65
  66       if (MPFR_UNLIKELY( diff_exp == 0 ))
  67         {
  68           while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
  69             {
  70               bn--;
  71               cn--;
  72               res += BITS_PER_MP_LIMB;
  73             }
  74
  75           if (MPFR_UNLIKELY (bn < 0))
  76             {
  77               if (MPFR_LIKELY (cn < 0)) /* b = c */
  78                 return 0;
  79
  80               bp = cp;
  81               bn = cn;
  82               cn = -1;
  83               sign = -1;
  84             }
  85
  86           if (MPFR_UNLIKELY (cn < 0))
  87             /* c discards exactly the upper part of b */
  88             {
  89               unsigned int z;
  90
  91               MPFR_ASSERTD (bn >= 0);
  92
  93               while (bp[bn] == 0)
  94                 {
  95                   if (--bn < 0) /* b = c */
  96                     return 0;
  97                   res += BITS_PER_MP_LIMB;
  98                 }
  99
 100               count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
 101               *cancel = res + z;
 102               return sign;
 103             }
 104
 105           MPFR_ASSERTD (bn >= 0);
 106           MPFR_ASSERTD (cn >= 0);
 107           MPFR_ASSERTD (bp[bn] != cp[cn]);
 108           if (bp[bn] < cp[cn])
 109             {
 110               mp_limb_t *tp;
 111               mp_size_t tn;
 112
 113               tp = bp; bp = cp; cp = tp;
 114               tn = bn; bn = cn; cn = tn;
 115               sign = -1;
 116             }
 117         }
 118     } /* MPFR_EXP(b) >= MPFR_EXP(c) */
 119   else /* MPFR_EXP(b) < MPFR_EXP(c) */
 120     {
 121       sign = -1;
 122       diff_exp = (mp_exp_unsigned_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b);
 123
 124       bp = MPFR_MANT(c);
 125       cp = MPFR_MANT(b);
 126
 127       bn = (MPFR_PREC(c) - 1) / BITS_PER_MP_LIMB;
 128       cn = (MPFR_PREC(b) - 1) / BITS_PER_MP_LIMB;
 129     }
 130
 131   /* now we have removed the identical upper limbs of b and c
 132      (can happen only when diff_exp = 0), and after the possible
 133      swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0,
 134      diff_exp = EXP(b) - EXP(c).
 135   */
 136
 137   if (MPFR_LIKELY (diff_exp < BITS_PER_MP_LIMB))
 138     {
 139       cc = cp[cn] >> diff_exp;
 140       /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
 141       if (diff_exp)
 142         lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
 143       cn--;
 144     }
 145   else
 146     diff_exp -= BITS_PER_MP_LIMB; /* cc = 0 */
 147
 148   dif = bp[bn--] - cc; /* necessarily dif >= 1 */
 149   MPFR_ASSERTD(dif >= 1);
 150
 151   /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */
 152
 153   while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
 154                         && (high_dif == 0) && (dif == 1)))
 155     { /* dif=1 implies diff_exp = 0 or 1 */
 156       bb = (bn >= 0) ? bp[bn] : 0;
 157       cc = lastc;
 158       if (cn >= 0)
 159         {
 160           if (diff_exp == 0)
 161             {
 162               cc += cp[cn];
 163             }
 164           else /* diff_exp = 1 */
 165             {
 166               cc += cp[cn] >> 1;
 167               lastc = cp[cn] << (BITS_PER_MP_LIMB - 1);
 168             }
 169         }
 170       else
 171         lastc = 0;
 172       high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
 173       bn--;
 174       cn--;
 175       res += BITS_PER_MP_LIMB;
 176     }
 177
 178   /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
 179
 180   if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */
 181     {
 182       res--;
 183       if (dif != 0)
 184         {
 185           *cancel = res;
 186           return sign;
 187         }
 188     }
 189   else /* high_dif == 0 */
 190     {
 191       unsigned int z;
 192
 193       count_leading_zeros(z, dif); /* dif > 1 here */
 194       res += z;
 195       if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - z - 1))))
 196         { /* dif is not a power of two */
 197           *cancel = res;
 198           return sign;
 199         }
 200     }
 201
 202   /* now result is res + (low(b) < low(c)) */
 203   while (MPFR_UNLIKELY (bn >= 0 && (cn >= 0 || lastc != 0)))
 204     {
 205       if (diff_exp >= BITS_PER_MP_LIMB)
 206         diff_exp -= BITS_PER_MP_LIMB;
 207       else
 208         {
 209           cc = lastc;
 210           if (cn >= 0)
 211             {
 212               cc += cp[cn] >> diff_exp;
 213               if (diff_exp != 0)
 214                 lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
 215             }
 216           else
 217             lastc = 0;
 218           cn--;
 219         }
 220       if (bp[bn] != cc)
 221         {
 222           *cancel = res + (bp[bn] < cc);
 223           return sign;
 224         }
 225       bn--;
 226     }
 227
 228   if (bn < 0)
 229     {
 230       if (lastc != 0)
 231         res++;
 232       else
 233         {
 234           while (cn >= 0 && cp[cn] == 0)
 235             cn--;
 236           if (cn >= 0)
 237             res++;
 238         }
 239     }
 240
 241   *cancel = res;
 242   return sign;
 243 }