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;
  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
 135   if (MPFR_LIKELY (diff_exp < BITS_PER_MP_LIMB))
 136     {
 137       cc = cp[cn] >> diff_exp;
 138       /* warning: a shift by BITS_PER_MP_LIMB may give wrong results */
 139       if (diff_exp)
 140         lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
 141       cn--;
 142     }
 143   else
 144     diff_exp -= BITS_PER_MP_LIMB;
 145
 146   dif = bp[bn--] - cc; /* necessarily dif >= 1 */
 147
 148   while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
 149                         && (high_dif == 0) && (dif == 1)))
 150     { /* dif=1 implies diff_exp = 0 or 1 */
 151       bb = (bn >= 0) ? bp[bn] : 0;
 152       cc = lastc;
 153       if (cn >= 0)
 154         {
 155           if (diff_exp == 0)
 156             {
 157               cc += cp[cn];
 158             }
 159           else /* diff_exp = 1 */
 160             {
 161               cc += cp[cn] >> 1;
 162               lastc = cp[cn] << (BITS_PER_MP_LIMB - 1);
 163             }
 164         }
 165       else
 166         lastc = 0;
 167       high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
 168       bn--;
 169       cn--;
 170       res += BITS_PER_MP_LIMB;
 171     }
 172
 173   /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
 174
 175   if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */
 176     {
 177       res--;
 178       if (dif != 0)
 179         {
 180           *cancel = res;
 181           return sign;
 182         }
 183     }
 184   else /* high_dif == 0 */
 185     {
 186       unsigned int z;
 187
 188       count_leading_zeros(z, dif); /* dif > 1 here */
 189       res += z;
 190       if (MPFR_LIKELY(dif != (MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - z - 1))))
 191         { /* dif is not a power of two */
 192           *cancel = res;
 193           return sign;
 194         }
 195     }
 196
 197   /* now result is res + (low(b) < low(c)) */
 198   while (MPFR_UNLIKELY (bn >= 0 && (cn >= 0 || lastc != 0)))
 199     {
 200       if (diff_exp >= BITS_PER_MP_LIMB)
 201         diff_exp -= BITS_PER_MP_LIMB;
 202       else
 203         {
 204           cc = lastc;
 205           if (cn >= 0)
 206             {
 207               cc += cp[cn] >> diff_exp;
 208               if (diff_exp != 0)
 209                 lastc = cp[cn] << (BITS_PER_MP_LIMB - diff_exp);
 210             }
 211           else
 212             lastc = 0;
 213           cn--;
 214         }
 215       if (bp[bn] != cc)
 216         {
 217           *cancel = res + (bp[bn] < cc);
 218           return sign;
 219         }
 220       bn--;
 221     }
 222
 223   if (bn < 0)
 224     {
 225       if (lastc != 0)
 226         res++;
 227       else
 228         {
 229           while (cn >= 0 && cp[cn] == 0)
 230             cn--;
 231           if (cn >= 0)
 232             res++;
 233         }
 234     }
 235
 236   *cancel = res;
 237   return sign;
 238 }