1/* mpfr_cmp2 -- exponent shift when subtracting two numbers.
2
3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 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
36int
37mpfr_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}
244