1 /* mpfr_sub1 -- internal function to perform a "real" subtraction
2
3 Copyright 2001-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 #include "mpfr-impl.h"
24
25 /* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
26 Returns 0 iff result is exact,
27 a negative value when the result is less than the exact value,
28 a positive value otherwise.
29 */
30
31 int
32 mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
33 {
34 int sign;
35 mpfr_exp_t diff_exp, exp_a, exp_b;
36 mpfr_prec_t cancel, cancel1;
37 mp_size_t cancel2, an, bn, cn, cn0;
38 mp_limb_t *ap, *bp, *cp;
39 mp_limb_t carry, bb, cc;
40 mpfr_prec_t aq, bq;
41 int inexact, shift_b, shift_c, add_exp = 0;
42 int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c),
43 negative if low(b) < low(c), positive if low(b) > low(c) */
44 int sh, k;
45 MPFR_TMP_DECL(marker);
46
47 MPFR_TMP_MARK(marker);
48 ap = MPFR_MANT(a);
49 an = MPFR_LIMB_SIZE(a);
50
51 (void) MPFR_GET_PREC (a);
52 (void) MPFR_GET_PREC (b);
53 (void) MPFR_GET_PREC (c);
54
55 sign = mpfr_cmp2 (b, c, &cancel);
56
57 if (MPFR_UNLIKELY(sign == 0))
58 {
59 MPFR_LOG_MSG (("sign=0\n", 0));
60 if (rnd_mode == MPFR_RNDD)
61 MPFR_SET_NEG (a);
62 else
63 MPFR_SET_POS (a);
64 MPFR_SET_ZERO (a);
65 MPFR_RET (0);
66 }
67
68 /* sign != 0, so that cancel has a valid value. */
69 MPFR_LOG_MSG (("sign=%d cancel=%Pd\n", sign, cancel));
70 MPFR_ASSERTD (cancel >= 0 && cancel <= MPFR_PREC_MAX);
71
72 /*
73 * If subtraction: sign(a) = sign * sign(b)
74 * If addition: sign(a) = sign of the larger argument in absolute value.
75 *
76 * Both cases can be simplified in:
77 * if (sign>0)
78 * if addition: sign(a) = sign * sign(b) = sign(b)
79 * if subtraction, b is greater, so sign(a) = sign(b)
80 * else
81 * if subtraction, sign(a) = - sign(b)
82 * if addition, sign(a) = sign(c) (since c is greater)
83 * But if it is an addition, sign(b) and sign(c) are opposed!
84 * So sign(a) = - sign(b)
85 */
86
87 if (sign < 0) /* swap b and c so that |b| > |c| */
88 {
89 mpfr_srcptr t;
90 MPFR_SET_OPPOSITE_SIGN (a,b);
91 t = b; b = c; c = t;
92 }
93 else
94 MPFR_SET_SAME_SIGN (a,b);
95
96 if (MPFR_UNLIKELY (MPFR_IS_UBF (b) || MPFR_IS_UBF (c)))
97 {
98 exp_b = MPFR_UBF_GET_EXP (b);
99 /* Early underflow detection. Rare, but a test is needed anyway
100 since in the "MAX (aq, bq) + 2 <= diff_exp" branch, the exponent
101 may decrease and MPFR_EXP_MIN would yield an integer overflow. */
102 if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1))
103 {
104 if (rnd_mode == MPFR_RNDN)
105 rnd_mode = MPFR_RNDZ;
106 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
107 }
108 diff_exp = mpfr_ubf_diff_exp (b, c);
109 MPFR_LOG_MSG (("UBF: exp_b=%" MPFR_EXP_FSPEC "d%s "
110 "diff_exp=%" MPFR_EXP_FSPEC "d%s\n",
111 (mpfr_eexp_t) exp_b,
112 exp_b == MPFR_EXP_MAX ? "=MPFR_EXP_MAX" : "",
113 (mpfr_eexp_t) diff_exp,
114 diff_exp == MPFR_EXP_MAX ? "=MPFR_EXP_MAX" : ""));
115 /* If diff_exp == MPFR_EXP_MAX, the actual value can be larger,
116 but anyway, since mpfr_exp_t >= mp_size_t, this will be the
117 case c small below, and the exact value does not matter. */
118 /* mpfr_set4 below used with MPFR_RNDF does not support UBF. */
119 if (rnd_mode == MPFR_RNDF)
120 rnd_mode = MPFR_RNDN;
121 }
122 else
123 {
124 exp_b = MPFR_GET_EXP (b);
125 diff_exp = exp_b - MPFR_GET_EXP (c);
126 }
127 MPFR_ASSERTD (diff_exp >= 0);
128
129 aq = MPFR_GET_PREC (a);
130 bq = MPFR_GET_PREC (b);
131
132 /* Check if c is too small.
133 A more precise test is to replace 2 by
134 (rnd == MPFR_RNDN) + mpfr_power2_raw (b)
135 but it is more expensive and not very useful */
136 if (MPFR_UNLIKELY (MAX (aq, bq) + 2 <= diff_exp))
137 {
138 MPFR_LOG_MSG (("case c small\n", 0));
139
140 /* Remember, we can't have an exact result! */
141 /* A.AAAAAAAAAAAAAAAAA
142 = B.BBBBBBBBBBBBBBB
143 - C.CCCCCCCCCCCCC */
144 /* A = S*ABS(B) +/- ulp(a) */
145
146 /* since we can't have an exact result, for RNDF we can truncate b */
147 if (rnd_mode == MPFR_RNDF)
148 return mpfr_set4 (a, b, MPFR_RNDZ, MPFR_SIGN (a));
149
150 exp_a = exp_b; /* may be any out-of-range value due to UBF */
151 MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), bq,
152 rnd_mode, MPFR_SIGN (a),
153 if (exp_a != MPFR_EXP_MAX)
154 exp_a ++);
155 MPFR_LOG_MSG (("inexact=%d\n", inexact));
156 if (inexact == 0 &&
157 /* a = b, but the exact value of b - c is a bit below. Then,
158 except for directed rounding similar to toward zero and
159 before overflow checking: a is the correctly rounded value
160 and since |b| - |c| < |a|, the ternary value value is given
161 by the sign of a. */
162 ! MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
163 {
164 MPFR_LOG_MSG (("c small, case 1\n", 0));
165 inexact = MPFR_INT_SIGN (a);
166 }
167 else if (inexact != 0 &&
168 /* A.AAAAAAAAAAAAAA
169 = B.BBBBBBBBBBBBBBB
170 - C.CCCCCCCCCCCCC */
171 /* It isn't exact, so PREC(b) > PREC(a) and the last
172 PREC(b)-PREC(a) bits of b are not all zeros.
173 Subtracting c from b will not have an effect on the rounding
174 except in case of a midpoint in the round-to-nearest mode,
175 when the even rounding was done away from zero instead of
176 toward zero.
177 In case of even rounding:
178 1.BBBBBBBBBBBBBx10
179 - 1.CCCCCCCCCCCC
180 = 1.BBBBBBBBBBBBBx01 Rounded to PREC(b)
181 = 1.BBBBBBBBBBBBBx Nearest / Rounded to PREC(a)
182 Set gives:
183 1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0)
184 1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
185 which means we get a wrong rounded result if x == 1,
186 i.e. inexact == MPFR_EVEN_INEX (for positive numbers). */
187 MPFR_LIKELY (inexact != MPFR_EVEN_INEX * MPFR_INT_SIGN (a)))
188 {
189 MPFR_LOG_MSG (("c small, case 2\n", 0));
190 /* nothing to do */
191 }
192 else
193 {
194 /* We need to take the value preceding |a|. We can't use
195 mpfr_nexttozero due to a possible out-of-range exponent.
196 But this will allow us to have more specific code. */
197 MPFR_LOG_MSG (("c small, case 3: correcting the value of a\n", 0));
198 sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq;
199 mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
200 if (MPFR_UNLIKELY (MPFR_LIMB_MSB (ap[an-1]) == 0))
201 {
202 exp_a --;
203 /* The following is valid whether an = 1 or an > 1. */
204 ap[an-1] |= MPFR_LIMB_HIGHBIT;
205 }
206 inexact = - MPFR_INT_SIGN (a);
207 }
208 /* The underflow case is possible only with UBF. The overflow case
209 is also possible with normal FP due to rounding. */
210 if (MPFR_UNLIKELY (exp_a > __gmpfr_emax))
211 return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
212 if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
213 {
214 if (rnd_mode == MPFR_RNDN &&
215 (exp_a < __gmpfr_emin - 1 ||
216 (inexact * MPFR_INT_SIGN (a) >= 0 && mpfr_powerof2_raw (a))))
217 rnd_mode = MPFR_RNDZ;
218 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
219 }
220 MPFR_SET_EXP (a, exp_a);
221 MPFR_RET (inexact);
222 }
223
224 /* reserve a space to store b aligned with the result, i.e. shifted by
225 (-cancel) % GMP_NUMB_BITS to the right */
226 bn = MPFR_LIMB_SIZE (b);
227 MPFR_UNSIGNED_MINUS_MODULO (shift_b, cancel);
228 cancel1 = (cancel + shift_b) / GMP_NUMB_BITS;
229
230 /* the high cancel1 limbs from b should not be taken into account */
231 if (MPFR_UNLIKELY (shift_b == 0))
232 {
233 bp = MPFR_MANT(b); /* no need of an extra space */
234 /* Ensure ap != bp */
235 if (MPFR_UNLIKELY (ap == bp))
236 {
237 bp = MPFR_TMP_LIMBS_ALLOC (bn);
238 MPN_COPY (bp, ap, bn);
239 }
240 }
241 else
242 {
243 bp = MPFR_TMP_LIMBS_ALLOC (bn + 1);
244 bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
245 }
246
247 /* reserve a space to store c aligned with the result, i.e. shifted by
248 (diff_exp-cancel) % GMP_NUMB_BITS to the right */
249 cn = MPFR_LIMB_SIZE (c);
250 if (IS_POW2 (GMP_NUMB_BITS))
251 shift_c = ((mpfr_uexp_t) diff_exp - cancel) % GMP_NUMB_BITS;
252 else
253 {
254 /* The above operation does not work if diff_exp - cancel < 0. */
255 shift_c = diff_exp - (cancel % GMP_NUMB_BITS);
256 shift_c = (shift_c + GMP_NUMB_BITS) % GMP_NUMB_BITS;
257 }
258 MPFR_ASSERTD (shift_c >= 0 && shift_c < GMP_NUMB_BITS);
259
260 if (MPFR_UNLIKELY(shift_c == 0))
261 {
262 cp = MPFR_MANT(c);
263 /* Ensure ap != cp */
264 if (ap == cp)
265 {
266 cp = MPFR_TMP_LIMBS_ALLOC (cn);
267 MPN_COPY(cp, ap, cn);
268 }
269 }
270 else
271 {
272 cp = MPFR_TMP_LIMBS_ALLOC (cn + 1);
273 cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
274 }
275
276 #if 0
277 MPFR_LOG_MSG (("rnd=%s shift_b=%d shift_c=%d diffexp=%" MPFR_EXP_FSPEC
278 "d\n", mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c,
279 (mpfr_eexp_t) diff_exp));
280 #endif
281
282 MPFR_ASSERTD (ap != cp);
283 MPFR_ASSERTD (bp != cp);
284
285 /* here we have shift_c = (diff_exp - cancel) % GMP_NUMB_BITS,
286 0 <= shift_c < GMP_NUMB_BITS
287 thus we want cancel2 = ceil((cancel - diff_exp) / GMP_NUMB_BITS) */
288
289 /* Possible optimization with a C99 compiler (i.e. well-defined
290 integer division): if MPFR_PREC_MAX is reduced to
291 ((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1) - GMP_NUMB_BITS + 1)
292 and diff_exp is of type mpfr_exp_t (no need for mpfr_uexp_t, since
293 the sum or difference of 2 exponents must be representable, as used
294 by the multiplication code), then the computation of cancel2 could
295 be simplified to
296 cancel2 = (cancel - (diff_exp - shift_c)) / GMP_NUMB_BITS;
297 because cancel, diff_exp and shift_c are all non-negative and
298 these variables are signed. */
299
300 MPFR_ASSERTD (cancel >= 0);
301 if (cancel >= diff_exp)
302 /* Note that cancel is signed and will be converted to mpfr_uexp_t
303 (type of diff_exp) in the expression below, so that this will
304 work even if cancel is very large and diff_exp = 0. */
305 cancel2 = (cancel - diff_exp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS;
306 else
307 cancel2 = - (mp_size_t) ((diff_exp - cancel) / GMP_NUMB_BITS);
308 /* the high cancel2 limbs from b should not be taken into account */
309 #if 0
310 MPFR_LOG_MSG (("cancel=%Pd cancel1=%Pd cancel2=%Pd\n",
311 cancel, cancel1, cancel2));
312 #endif
313
314 /* ap[an-1] ap[0]
315 <----------------+-----------|---->
316 <----------PREC(a)----------><-sh->
317 cancel1
318 limbs bp[bn-cancel1-1]
319 <--...-----><----------------+-----------+----------->
320 cancel2
321 limbs cp[cn-cancel2-1] cancel2 >= 0
322 <--...--><----------------+----------------+---------------->
323 (-cancel2) cancel2 < 0
324 limbs <----------------+---------------->
325 */
326
327 /* first part: put in ap[0..an-1] the value of high(b) - high(c),
328 where high(b) consists of the high an+cancel1 limbs of b,
329 and high(c) consists of the high an+cancel2 limbs of c.
330 */
331
332 /* copy high(b) into a */
333 if (MPFR_LIKELY(an + (mp_size_t) cancel1 <= bn))
334 /* a: <----------------+-----------|---->
335 b: <-----------------------------------------> */
336 MPN_COPY (ap, bp + bn - (an + cancel1), an);
337 else
338 /* a: <----------------+-----------|---->
339 b: <-------------------------> */
340 if ((mp_size_t) cancel1 < bn) /* otherwise b does not overlap with a */
341 {
342 MPN_ZERO (ap, an + cancel1 - bn);
343 MPN_COPY (ap + (an + cancel1 - bn), bp, bn - cancel1);
344 }
345 else
346 MPN_ZERO (ap, an);
347
348 /* subtract high(c) */
349 if (MPFR_LIKELY(an + cancel2 > 0)) /* otherwise c does not overlap with a */
350 {
351 mp_limb_t *ap2;
352
353 if (cancel2 >= 0)
354 {
355 if (an + cancel2 <= cn)
356 /* a: <----------------------------->
357 c: <-----------------------------------------> */
358 mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
359 else
360 /* a: <---------------------------->
361 c: <-------------------------> */
362 {
363 ap2 = ap + an + (cancel2 - cn);
364 if (cn > cancel2)
365 mpn_sub_n (ap2, ap2, cp, cn - cancel2);
366 }
367 }
368 else /* cancel2 < 0 */
369 {
370 mp_limb_t borrow;
371
372 if (an + cancel2 <= cn)
373 /* a: <----------------------------->
374 c: <-----------------------------> */
375 borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2),
376 an + cancel2);
377 else
378 /* a: <---------------------------->
379 c: <----------------> */
380 {
381 ap2 = ap + an + cancel2 - cn;
382 borrow = mpn_sub_n (ap2, ap2, cp, cn);
383 }
384 ap2 = ap + an + cancel2;
385 mpn_sub_1 (ap2, ap2, -cancel2, borrow);
386 }
387 }
388
389 /* now perform rounding */
390 sh = (mpfr_prec_t) an * GMP_NUMB_BITS - aq;
391 /* last unused bits from a */
392 carry = ap[0] & MPFR_LIMB_MASK (sh);
393 ap[0] -= carry;
394
395 if (rnd_mode == MPFR_RNDF)
396 {
397 inexact = 0;
398 /* truncating is always correct since -1 ulp < low(b) - low(c) < 1 ulp */
399 goto truncate;
400 }
401 else if (rnd_mode == MPFR_RNDN)
402 {
403 if (MPFR_LIKELY(sh))
404 {
405 /* can decide except when carry = 2^(sh-1) [middle]
406 or carry = 0 [truncate, but cannot decide inexact flag] */
407 if (carry > (MPFR_LIMB_ONE << (sh - 1)))
408 goto add_one_ulp;
409 else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1))))
410 {
411 inexact = -1; /* result if smaller than exact value */
412 goto truncate;
413 }
414 /* now carry = 2^(sh-1), in which case cmp_low=2,
415 or carry = 0, in which case cmp_low=0 */
416 cmp_low = (carry == 0) ? 0 : 2;
417 }
418 }
419 else /* directed rounding: set rnd_mode to RNDZ iff toward zero */
420 {
421 if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a)))
422 rnd_mode = MPFR_RNDZ;
423
424 if (carry)
425 {
426 if (rnd_mode == MPFR_RNDZ)
427 {
428 inexact = -1;
429 goto truncate;
430 }
431 else /* round away */
432 goto add_one_ulp;
433 }
434 }
435
436 /* we have to consider the low (bn - (an+cancel1)) limbs from b,
437 and the (cn - (an+cancel2)) limbs from c. */
438 bn -= an + cancel1;
439 cn0 = cn;
440 cn -= an + cancel2;
441
442 #if 0
443 MPFR_LOG_MSG (("last sh=%d bits from a are %Mu, bn=%Pd, cn=%Pd\n",
444 sh, carry, (mpfr_prec_t) bn, (mpfr_prec_t) cn));
445 #endif
446
447 /* for rounding to nearest, we couldn't conclude up to here in the following
448 cases:
449 1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp
450 or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp
451 2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1):
452 -0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp
453 we can't decide the rounding, in that case cmp_low=2:
454 either we truncate and flag=-1, or we add one ulp and flag=1
455 3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to
456 truncate but we can't decide the ternary value, here cmp_low=0:
457 -0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp
458 we always truncate and inexact can be any of -1,0,1
459 */
460
461 /* note: here cn might exceed cn0, in which case we consider a zero limb */
462 for (k = 0; (bn > 0) || (cn > 0); k = 1)
463 {
464 /* if cmp_low < 0, we know low(b) - low(c) < 0
465 if cmp_low > 0, we know low(b) - low(c) > 0
466 (more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far)
467 if cmp_low = 0, so far low(b) - low(c) = 0 */
468
469 /* get next limbs */
470 bb = (bn > 0) ? bp[--bn] : 0;
471 if ((cn > 0) && (cn-- <= cn0))
472 cc = cp[cn];
473 else
474 cc = 0;
475
476 /* cmp_low compares low(b) and low(c) */
477 if (cmp_low == 0) /* case 1 or 3 */
478 cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0;
479
480 /* Case 1 for k=0 splits into 7 subcases:
481 1a: bb > cc + half
482 1b: bb = cc + half
483 1c: 0 < bb - cc < half
484 1d: bb = cc
485 1e: -half < bb - cc < 0
486 1f: bb - cc = -half
487 1g: bb - cc < -half
488
489 Case 2 splits into 3 subcases:
490 2a: bb > cc
491 2b: bb = cc
492 2c: bb < cc
493
494 Case 3 splits into 3 subcases:
495 3a: bb > cc
496 3b: bb = cc
497 3c: bb < cc
498 */
499
500 /* the case rounding to nearest with sh=0 is special since one couldn't
501 subtract above 1/2 ulp in the trailing limb of the result */
502 if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */
503 {
504 mp_limb_t half = MPFR_LIMB_HIGHBIT;
505
506 /* add one ulp if bb > cc + half
507 truncate if cc - half < bb < cc + half
508 sub one ulp if bb < cc - half
509 */
510
511 if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0,
512 cases 1e, 1f and 1g */
513 {
514 if (cc >= half)
515 cc -= half;
516 else /* since bb < cc < half, bb+half < 2*half */
517 bb += half;
518 /* now we have bb < cc + half:
519 we have to subtract one ulp if bb < cc,
520 and truncate if bb > cc */
521 }
522 else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */
523 {
524 if (cc < half)
525 cc += half;
526 else /* since bb >= cc >= half, bb - half >= 0 */
527 bb -= half;
528 /* now we have bb > cc - half: we have to add one ulp if bb > cc,
529 and truncate if bb < cc */
530 if (cmp_low > 0)
531 cmp_low = 2;
532 }
533 }
534
535 #if 0
536 MPFR_LOG_MSG (("k=%d bb=%Mu cc=%Mu cmp_low=%d\n", k, bb, cc, cmp_low));
537 #endif
538
539 if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract
540 one ulp */
541 {
542 if (rnd_mode == MPFR_RNDZ)
543 goto sub_one_ulp; /* set inexact=-1 */
544 else if (rnd_mode != MPFR_RNDN) /* round away */
545 {
546 inexact = 1;
547 goto truncate;
548 }
549 else /* round to nearest */
550 {
551 /* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0,
552 whatever the value of sh.
553 If sh>0, then cmp_low < 0 implies that the initial neglected
554 sh bits were 0 (otherwise cmp_low=2 initially), thus the
555 weight of the new bits is less than 0.5 ulp too.
556 If k > 0 (and sh=0) this means that either the first neglected
557 limbs bb and cc were equal (thus cmp_low was 0 for k=0),
558 or we had bb - cc = -0.5 ulp or 0.5 ulp.
559 The last case is not possible here since we would have
560 cmp_low > 0 which is sticky.
561 In the first case (where we have cmp_low = -1), we truncate,
562 whereas in the 2nd case we have cmp_low = -2 and we subtract
563 one ulp.
564 */
565 if (bb > cc || sh > 0 || cmp_low == -1)
566 { /* -0.5 ulp < low(b)-low(c) < 0,
567 bb > cc corresponds to cases 1e and 1f1
568 sh > 0 corresponds to cases 3c and 3b3
569 cmp_low = -1 corresponds to case 1d3 (also 3b3) */
570 inexact = 1;
571 goto truncate;
572 }
573 else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp,
574 this corresponds to cases 1g and 1f3 */
575 goto sub_one_ulp;
576 /* the only case where we can't conclude is sh=0 and bb=cc,
577 i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus
578 we don't know if we must truncate or subtract one ulp.
579 Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to
580 now, since low(b) - low(c) > 1/2^sh */
581 }
582 }
583 else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or
584 add one ulp */
585 {
586 if (rnd_mode == MPFR_RNDZ)
587 {
588 inexact = -1;
589 goto truncate;
590 }
591 else if (rnd_mode != MPFR_RNDN) /* round away */
592 goto add_one_ulp;
593 else /* round to nearest */
594 {
595 if (bb > cc)
596 {
597 /* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp,
598 and similarly when cmp_low=2 */
599 if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */
600 goto add_one_ulp;
601 /* sh > 0 and cmp_low > 0: this implies that the sh initial
602 neglected bits were 0, and the remaining low(b)-low(c)>0,
603 but its weight is less than 0.5 ulp */
604 else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to
605 cases 3a, 1d1 and 3b1 */
606 {
607 inexact = -1;
608 goto truncate;
609 }
610 }
611 else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c,
612 1b3, 2b3 and 2c */
613 {
614 inexact = -1;
615 goto truncate;
616 }
617 /* the only case where we can't conclude is bb=cc, i.e.,
618 low(b) - low(c) = 0.5 ulp (up to now), thus we don't know
619 if we must truncate or add one ulp. */
620 }
621 }
622 /* after k=0, we cannot conclude in the following cases, we split them
623 according to the values of bb and cc for k=1:
624 1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp]
625 1b1. bb > cc: add one ulp, inex = 1
626 1b2: bb = cc: cannot conclude
627 1b3: bb < cc: truncate, inex = -1
628 1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0]
629 1d1: bb > cc: truncate, inex = -1
630 1d2: bb = cc: cannot conclude
631 1d3: bb < cc: truncate, inex = +1
632 1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp]
633 1f1: bb > cc: truncate, inex = +1
634 1f2: bb = cc: cannot conclude
635 1f3: bb < cc: sub one ulp, inex = -1
636 2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp]
637 2b1. bb > cc: add one ulp, inex = 1
638 2b2: bb = cc: cannot conclude
639 2b3: bb < cc: truncate, inex = -1
640 3b. sh > 0 and cmp_low = 0 [around 0]
641 3b1. bb > cc: truncate, inex = -1
642 3b2: bb = cc: cannot conclude
643 3b3: bb < cc: truncate, inex = +1
644 */
645 }
646
647 if ((rnd_mode == MPFR_RNDN) && cmp_low != 0)
648 {
649 /* even rounding rule */
650 if ((ap[0] >> sh) & 1)
651 {
652 if (cmp_low < 0)
653 goto sub_one_ulp;
654 else
655 goto add_one_ulp;
656 }
657 else
658 inexact = (cmp_low > 0) ? -1 : 1;
659 }
660 else
661 inexact = 0;
662 goto truncate;
663
664 sub_one_ulp: /* sub one unit in last place to a */
665 mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
666 inexact = -1;
667 goto end_of_sub;
668
669 add_one_ulp: /* add one unit in last place to a */
670 if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
671 /* result is a power of 2: 11111111111111 + 1 = 1000000000000000 */
672 {
673 ap[an-1] = MPFR_LIMB_HIGHBIT;
674 add_exp = 1;
675 }
676 inexact = 1; /* result larger than exact value */
677
678 truncate:
679 if (MPFR_UNLIKELY((ap[an-1] >> (GMP_NUMB_BITS - 1)) == 0))
680 /* case 1 - epsilon */
681 {
682 ap[an-1] = MPFR_LIMB_HIGHBIT;
683 add_exp = 1;
684 }
685
686 end_of_sub:
687 /* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
688 care of underflows/overflows in that computation, and of the allowed
689 exponent range */
690 MPFR_TMP_FREE (marker);
691 if (MPFR_LIKELY(cancel))
692 {
693 cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
694 MPFR_ASSERTD (cancel >= 0);
695 /* Detect an underflow case to avoid a possible integer overflow
696 with UBF in the computation of exp_a. */
697 if (MPFR_UNLIKELY (exp_b < __gmpfr_emin - 1))
698 {
699 if (rnd_mode == MPFR_RNDN)
700 rnd_mode = MPFR_RNDZ;
701 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
702 }
703 exp_a = exp_b - cancel;
704 /* The following assertion corresponds to a limitation of the MPFR
705 implementation. It may fail with a 32-bit ABI and huge precisions,
706 but this is practically impossible with a 64-bit ABI. This kind
707 of issue is not specific to this function. */
708 MPFR_ASSERTN (exp_b != MPFR_EXP_MAX || exp_a > __gmpfr_emax);
709 if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
710 {
711 underflow:
712 if (rnd_mode == MPFR_RNDN &&
713 (exp_a < __gmpfr_emin - 1 ||
714 (inexact >= 0 && mpfr_powerof2_raw (a))))
715 rnd_mode = MPFR_RNDZ;
716 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
717 }
718 /* We cannot have an overflow here, except for UBFs. Indeed:
719 exp_a = exp_b - cancel + add_exp <= emax - 1 + 1 <= emax.
720 For UBFs, we can have exp_b > emax. */
721 if (exp_a > __gmpfr_emax)
722 {
723 MPFR_ASSERTD (exp_b > __gmpfr_emax); /* since exp_b >= exp_a */
724 return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
725 }
726 }
727 else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
728 {
729 /* in case cancel = 0, add_exp can still be 1, in case b is just
730 below a power of two, c is very small, prec(a) < prec(b),
731 and rnd=away or nearest */
732 MPFR_ASSERTD (add_exp == 0 || add_exp == 1);
733 /* Overflow iff exp_b + add_exp > __gmpfr_emax in Z, but we do
734 a subtraction below to avoid a potential integer overflow in
735 the case exp_b == MPFR_EXP_MAX. */
736 if (MPFR_UNLIKELY (exp_b > __gmpfr_emax - add_exp))
737 return mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
738 exp_a = exp_b + add_exp;
739 /* Warning: an underflow can happen for UBFs, for example when
740 mpfr_add is called from mpfr_fmma or mpfr_fmms. */
741 if (MPFR_UNLIKELY (exp_a < __gmpfr_emin))
742 goto underflow;
743 MPFR_ASSERTD (exp_a >= __gmpfr_emin);
744 }
745 MPFR_SET_EXP (a, exp_a);
746 /* check that result is msb-normalized */
747 MPFR_ASSERTD(ap[an-1] > ~ap[an-1]);
748 MPFR_RET (inexact * MPFR_INT_SIGN (a));
749 }