1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
2
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC 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
mpfr_cmp2(mpfr_srcptr b,mpfr_srcptr c,mpfr_prec_t * cancel)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 }
244