contrib/mpfr/src/cmp2.c

   1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel 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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, 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 (-1, 0, 1).
  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, mpfr_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   mpfr_uexp_t diff_exp;
  42   mpfr_prec_t res = 0;
  43   int sign;
  44
  45   /* b=c should not happen, since cmp2 is called only from agm (with
  46      different variables) and from sub1 (if b=c, then sub1sp would be
  47      called instead). So, no need for a particular optimization here. */
  48
  49   /* the cases b=0 or c=0 are also treated apart in agm and sub
  50      (which calls sub1) */
  51   MPFR_ASSERTD (MPFR_IS_PURE_FP(b));
  52   MPFR_ASSERTD (MPFR_IS_PURE_FP(c));
  53
  54   if (MPFR_GET_EXP (b) >= MPFR_GET_EXP (c))
  55     {
  56       sign = 1;
  57       diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
  58
  59       bp = MPFR_MANT(b);
  60       cp = MPFR_MANT(c);
  61
  62       bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
  63       cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS; /* # of limbs of c minus 1 */
  64
  65       if (MPFR_UNLIKELY( diff_exp == 0 ))
  66         {
  67           while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
  68             {
  69               bn--;
  70               cn--;
  71               res += GMP_NUMB_BITS;
  72             }
  73
  74           if (MPFR_UNLIKELY (bn < 0))
  75             {
  76               if (MPFR_LIKELY (cn < 0)) /* b = c */
  77                 return 0;
  78
  79               bp = cp;
  80               bn = cn;
  81               cn = -1;
  82               sign = -1;
  83             }
  84
  85           if (MPFR_UNLIKELY (cn < 0))
  86             /* c discards exactly the upper part of b */
  87             {
  88               unsigned int z;
  89
  90               MPFR_ASSERTD (bn >= 0);
  91
  92               while (bp[bn] == 0)
  93                 {
  94                   if (--bn < 0) /* b = c */
  95                     return 0;
  96                   res += GMP_NUMB_BITS;
  97                 }
  98
  99               count_leading_zeros(z, bp[bn]); /* bp[bn] <> 0 */
 100               *cancel = res + z;
 101               return sign;
 102             }
 103
 104           MPFR_ASSERTD (bn >= 0);
 105           MPFR_ASSERTD (cn >= 0);
 106           MPFR_ASSERTD (bp[bn] != cp[cn]);
 107           if (bp[bn] < cp[cn])
 108             {
 109               mp_limb_t *tp;
 110               mp_size_t tn;
 111
 112               tp = bp; bp = cp; cp = tp;
 113               tn = bn; bn = cn; cn = tn;
 114               sign = -1;
 115             }
 116         }
 117     } /* MPFR_EXP(b) >= MPFR_EXP(c) */
 118   else /* MPFR_EXP(b) < MPFR_EXP(c) */
 119     {
 120       sign = -1;
 121       diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (c) - MPFR_GET_EXP (b);
 122
 123       bp = MPFR_MANT(c);
 124       cp = MPFR_MANT(b);
 125
 126       bn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS;
 127       cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
 128     }
 129
 130   /* now we have removed the identical upper limbs of b and c
 131      (can happen only when diff_exp = 0), and after the possible
 132      swap, we have |b| > |c|: bp[bn] > cc, bn >= 0, cn >= 0,
 133      diff_exp = EXP(b) - EXP(c).
 134   */
 135
 136   if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS))
 137     {
 138       cc = cp[cn] >> diff_exp;
 139       /* warning: a shift by GMP_NUMB_BITS may give wrong results */
 140       if (diff_exp)
 141         lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp);
 142       cn--;
 143     }
 144   else
 145     diff_exp -= GMP_NUMB_BITS; /* cc = 0 */
 146
 147   dif = bp[bn--] - cc; /* necessarily dif >= 1 */
 148   MPFR_ASSERTD(dif >= 1);
 149
 150   /* now high_dif = 0, dif >= 1, lastc is the neglected part of cp[cn+1] */
 151
 152   while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
 153                         && (high_dif == 0) && (dif == 1)))
 154     { /* dif=1 implies diff_exp = 0 or 1 */
 155       bb = (bn >= 0) ? bp[bn] : 0;
 156       cc = lastc;
 157       if (cn >= 0)
 158         {
 159           if (diff_exp == 0)
 160             {
 161               cc += cp[cn];
 162             }
 163           else /* diff_exp = 1 */
 164             {
 165               cc += cp[cn] >> 1;
 166               lastc = cp[cn] << (GMP_NUMB_BITS - 1);
 167             }
 168         }
 169       else
 170         lastc = 0;
 171       high_dif = 1 - mpn_sub_n (&dif, &bb, &cc, 1);
 172       bn--;
 173       cn--;
 174       res += GMP_NUMB_BITS;
 175     }
 176
 177   /* (cn<0 and lastc=0) or (high_dif,dif)<>(0,1) */
 178
 179   if (MPFR_UNLIKELY (high_dif != 0)) /* high_dif == 1 */
 180     {
 181       res--;
 182       MPFR_ASSERTD (res >= 0);
 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 << (GMP_NUMB_BITS - 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 >= GMP_NUMB_BITS)
 206         diff_exp -= GMP_NUMB_BITS;
 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] << (GMP_NUMB_BITS - 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 }