1 /* mpfr_cmp2 -- exponent shift when subtracting two numbers.
2
3 Copyright 1999-2004, 2006-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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) stores
33 EXP(max(|b|,|c|)) - EXP(|b| - |c|) in *cancel.
34
35 One necessarily has 0 <= cancel <= max(PREC(b),PREC(c)), so that this
36 value is representable in a mpfr_prec_t. Note that in the code, the
37 maximum intermediate value is cancel + 1, but since MPFR_PREC_MAX is
38 not the maximum value of mpfr_prec_t, there is no integer overflow.
39 */
40
41 int
42 mpfr_cmp2 (mpfr_srcptr b, mpfr_srcptr c, mpfr_prec_t *cancel)
43 {
44 mp_limb_t *bp, *cp, bb, cc, lastc, dif;
45 int high_dif; /* manipulated like a boolean */
46 mp_size_t bn, cn;
47 mpfr_exp_t sdiff_exp;
48 mpfr_uexp_t diff_exp;
49 mpfr_prec_t res = 0; /* will be the number of canceled bits (*cancel) */
50 int sign;
51
52 /* b=c should not happen, since cmp2 is called only from agm (with
53 different variables) and from sub1 (if b=c, then sub1sp would be
54 called instead). So, no need for a particular optimization here. */
55
56 /* the cases b=0 or c=0 are also treated apart in agm and sub
57 (which calls sub1) */
58 MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
59 MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
60
61 sdiff_exp = MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)) ?
62 mpfr_ubf_diff_exp (b, c) : MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
63
64 /* The returned result is saturated to [MPFR_EXP_MIN,MPFR_EXP_MAX],
65 which is the range of the mpfr_exp_t type. But under the condition
66 below, since |MPFR_EXP_MIN| >= MPFR_EXP_MAX, the value of cancel
67 will not be affected: by symmetry (as done in the code), assume
68 |b| >= |c|; if EXP(b) - EXP(c) >= MPFR_EXP_MAX, then |c| < ulp(b),
69 so that the value of cancel is 0, or 1 if |b| is a power of 2,
70 whatever the exact value of EXP(b) - EXP(c). */
71 MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX);
72
73 if (sdiff_exp >= 0)
74 {
75 sign = 1; /* assumes |b| > |c|; will be changed if not. */
76 diff_exp = sdiff_exp;
77
78 bp = MPFR_MANT(b);
79 cp = MPFR_MANT(c);
80
81 /* index of the most significant limb of b and c */
82 bn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
83 cn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS;
84
85 /* If diff_exp != 0, i.e. diff_exp > 0, then |b| > |c|. Otherwise... */
86 if (diff_exp == 0)
87 {
88 /* Skip the identical most significant limbs, adding GMP_NUMB_BITS
89 to the number of canceled bits at each iteration. */
90 while (bn >= 0 && cn >= 0 && bp[bn] == cp[cn])
91 {
92 bn--;
93 cn--;
94 res += GMP_NUMB_BITS;
95 }
96
97 if (MPFR_UNLIKELY (bn < 0))
98 {
99 if (MPFR_LIKELY (cn < 0)) /* |b| = |c| */
100 return 0;
101
102 /* b has been read entirely, but not c. Thus |b| <= |c|.
103 Swap (bp,bn) and (cp,cn), and take the opposite sign
104 for the symmetric case below (simulating a swap).
105 Note: cp will not be used, thus is not assigned; and
106 "cn = -1;" is necessary to enter the following "if"
107 (probably less confusing than a "goto"). */
108 bp = cp;
109 bn = cn;
110 cn = -1;
111 sign = -1;
112 }
113
114 if (MPFR_UNLIKELY (cn < 0))
115 /* c discards exactly the upper part of b */
116 {
117 int z;
118
119 MPFR_ASSERTD (bn >= 0);
120
121 /* Skip null limbs of b (= non-represented null limbs of c),
122 adding GMP_NUMB_BITS to the number of canceled bits at
123 each iteration. */
124 while (bp[bn] == 0)
125 {
126 if (--bn < 0) /* |b| = |c| */
127 return 0;
128 res += GMP_NUMB_BITS;
129 }
130
131 count_leading_zeros (z, bp[bn]); /* bp[bn] != 0 */
132 *cancel = res + z;
133 return sign;
134 }
135
136 MPFR_ASSERTD (bn >= 0);
137 MPFR_ASSERTD (cn >= 0);
138 MPFR_ASSERTD (bp[bn] != cp[cn]);
139
140 /* |b| != |c|. If |b| < |c|: swap (bp,bn) and (cp,cn),
141 and take the opposite sign. */
142 if (bp[bn] < cp[cn])
143 {
144 mp_limb_t *tp;
145 mp_size_t tn;
146
147 tp = bp; bp = cp; cp = tp;
148 tn = bn; bn = cn; cn = tn;
149 sign = -1;
150 }
151 }
152 } /* MPFR_EXP(b) >= MPFR_EXP(c) */
153 else /* MPFR_EXP(b) < MPFR_EXP(c) */
154 {
155 /* We necessarily have |b| < |c|. Simulate a swap by reading the
156 parameters so that |(bp,bn)| > |(cp,cn)|. */
157
158 sign = -1;
159 diff_exp = - (mpfr_uexp_t) sdiff_exp;
160
161 bp = MPFR_MANT(c);
162 cp = MPFR_MANT(b);
163
164 bn = (MPFR_PREC(c) - 1) / GMP_NUMB_BITS;
165 cn = (MPFR_PREC(b) - 1) / GMP_NUMB_BITS;
166 }
167
168 /* Now we have removed the identical upper limbs of b and c
169 (when diff_exp = 0), and after the possible swap, we have |b| > |c|,
170 where b is represented by (bp,bn) and c is represented by (cp,cn).
171 The value diff_exp = EXP(b) - EXP(c) can be regarded as the number
172 of leading zeros of c, when aligned with b. */
173
174 /* When a limb of c is read from memory, the part that is not taken
175 into account for the operation with a limb bp[bn] of b will be put
176 in lastc, shifted to the leftmost part (for alignment with b):
177 [-------- bp[bn] --------][------- bp[bn-1] -------]
178 [-- old_lastc --][-------- cp[cn] --------]
179 [-- new_lastc --]
180 Note: if diff_exp == 0, then lastc will always remain 0. */
181 lastc = 0;
182
183 /* Compute the next limb difference, which cannot be 0 (dif >= 1). */
184
185 if (MPFR_LIKELY (diff_exp < GMP_NUMB_BITS))
186 {
187 cc = cp[cn] >> diff_exp;
188 /* warning: a shift by GMP_NUMB_BITS is not allowed by ISO C */
189 if (diff_exp != 0)
190 lastc = cp[cn] << (GMP_NUMB_BITS - diff_exp);
191 cn--;
192 }
193 else
194 {
195 cc = 0;
196 diff_exp -= GMP_NUMB_BITS; /* remove GMP_NUMB_BITS leading zeros */
197 }
198
199 MPFR_ASSERTD (bp[bn] >= cc); /* no borrow out in subtraction below */
200 dif = bp[bn--] - cc;
201 MPFR_ASSERTD (dif >= 1);
202 high_dif = 0;
203
204 /* The current difference, here and later, is expressed under the form
205 [high_dif][dif], where high_dif is 0 or 1, and dif is a limb.
206 Here, since we have computed a difference of limbs (with b >= c),
207 high_dif = 0. */
208
209 /* One needs to accumulate canceled bits for the remaining case where
210 b and c are close to each other due to a long borrow propagation:
211 b = [common part]1000...000[low(b)]
212 c = [common part]0111...111[low(c)]
213 After eliminating the common part above, we have computed a difference
214 of the most significant parts, which has been stored in [high_dif][dif]
215 with high_dif = 0. We will loop as long as the currently computed
216 difference [high_dif][dif] = 1 (it is >= 1 by construction). The
217 computation of the difference will be:
218 1bbb...bbb
219 - ccc...ccc
220 where the leading 1 before bbb...bbb corresponds to [high_dif][dif]
221 at the beginning of the loop. We will exit the loop also when c has
222 entirely been taken into account as cancellation is no longer possible
223 in this case (it is no longer possible to cancel the leading 1).
224 Note: We can enter the loop only with diff_exp = 0 (with a non-empty
225 common part, partly or entirely removed) or with diff_exp = 1 (with
226 an empty common part). Indeed, if diff_exp > 1, then no limbs have
227 been skipped, so that bp[bn] had its MSB equal to 1 and the most two
228 significant bits of cc are 0, which implies that dif > 1. */
229
230 while (MPFR_UNLIKELY ((cn >= 0 || lastc != 0)
231 && high_dif == 0 && dif == 1))
232 {
233 /* Since we consider the next limb, we assume a cancellation of
234 GMP_NUMB_BITS (the new exponent of the difference now being the
235 one of the MSB of the next limb). But if the leading 1 remains
236 1 in the difference (i.e. high_dif = 1 at the end of the loop),
237 then we will need to decrease res. */
238 res += GMP_NUMB_BITS;
239 MPFR_ASSERTD (diff_exp <= 1); /* see comment before the loop */
240 bb = bn >= 0 ? bp[bn--] : 0; /* next limb of b or non-represented 0 */
241 if (MPFR_UNLIKELY (cn < 0))
242 {
243 cc = lastc;
244 lastc = 0;
245 }
246 else if (diff_exp == 0)
247 {
248 cc = cp[cn--];
249 }
250 else
251 {
252 MPFR_ASSERTD (diff_exp == 1);
253 MPFR_ASSERTD (lastc == 0 || lastc == MPFR_LIMB_HIGHBIT);
254 cc = lastc + (cp[cn] >> 1);
255 lastc = cp[cn--] << (GMP_NUMB_BITS - 1);
256 }
257 dif = bb - cc;
258 high_dif = bb >= cc;
259 }
260
261 /* Now, c has entirely been taken into account or [high_dif][dif] > 1.
262 In any case, [high_dif][dif] >= 1 by construction.
263 First, we determine the currently number of canceled bits,
264 corresponding to the exponent of the current difference.
265 The trailing bits of c, if any, can still decrease the exponent of
266 the difference when [high_dif][dif] is a power of two, but since
267 [high_dif][dif] > 1 in this case, by not more than 1. */
268
269 if (high_dif != 0) /* high_dif == 1 */
270 {
271 res--; /* see comment at the beginning of the above loop */
272 MPFR_ASSERTD (res >= 0);
273 /* Terminate if [high_dif][dif] is not a power of two. */
274 if (MPFR_LIKELY (dif != 0))
275 goto end;
276 }
277 else /* high_dif == 0 */
278 {
279 int z;
280
281 MPFR_ASSERTD (dif >= 1); /* [high_dif][dif] >= 1 */
282 count_leading_zeros (z, dif);
283 res += z;
284 /* Terminate if [high_dif][dif] is not a power of two. */
285 if (MPFR_LIKELY (NOT_POW2 (dif)))
286 goto end;
287 }
288
289 /* Now, the result will be res + (low(b) < low(c)). */
290
291 /* If c has entirely been taken into account, it can no longer modify
292 the current result. */
293 if (cn < 0 && lastc == 0)
294 goto end;
295
296 for (; bn >= 0 ; bn--)
297 {
298 if (diff_exp >= GMP_NUMB_BITS)
299 {
300 diff_exp -= GMP_NUMB_BITS;
301 MPFR_ASSERTD (cc == 0);
302 }
303 else if (MPFR_UNLIKELY (cn < 0))
304 {
305 cc = lastc;
306 lastc = 0;
307 }
308 else if (diff_exp == 0)
309 {
310 cc = cp[cn--];
311 }
312 else
313 {
314 MPFR_ASSERTD (diff_exp >= 1 && diff_exp < GMP_NUMB_BITS);
315 cc = lastc + (cp[cn] >> diff_exp);
316 lastc = cp[cn--] << (GMP_NUMB_BITS - diff_exp);
317 }
318
319 if (bp[bn] != cc)
320 {
321 res += bp[bn] < cc;
322 goto end;
323 }
324 }
325
326 /* b has entirely been read. Determine whether the trailing part of c
327 is non-zero. */
328
329 if (lastc != 0)
330 res++;
331 else
332 {
333 while (cn >= 0 && cp[cn] == 0)
334 cn--;
335 if (cn >= 0)
336 res++;
337 }
338
339 end:
340 *cancel = res;
341 return sign;
342 }