1 /* mpfr_add1 -- internal function to perform a "real" addition
2
3 Copyright 1999-2020 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|), assuming that b and c
26 are not NaN, Inf, nor zero. Assumes EXP(b) >= EXP(c).
27 */
28 MPFR_HOT_FUNCTION_ATTR int
mpfr_add1(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)29 mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
30 {
31 mp_limb_t *ap, *bp, *cp;
32 mpfr_prec_t aq, bq, cq, aq2;
33 mp_size_t an, bn, cn;
34 mpfr_exp_t difw, exp, diff_exp;
35 int sh, rb, fb, inex;
36 MPFR_TMP_DECL(marker);
37
38 MPFR_ASSERTD (MPFR_IS_PURE_UBF (b));
39 MPFR_ASSERTD (MPFR_IS_PURE_UBF (c));
40 MPFR_ASSERTD (! MPFR_UBF_EXP_LESS_P (b, c));
41
42 if (MPFR_UNLIKELY (MPFR_IS_UBF (b)))
43 {
44 exp = MPFR_UBF_GET_EXP (b);
45 if (exp > __gmpfr_emax)
46 return mpfr_overflow (a, rnd_mode, MPFR_SIGN (b));;
47 }
48 else
49 exp = MPFR_GET_EXP (b);
50
51 MPFR_ASSERTD (exp <= __gmpfr_emax);
52
53 MPFR_TMP_MARK(marker);
54
55 aq = MPFR_GET_PREC (a);
56 bq = MPFR_GET_PREC (b);
57 cq = MPFR_GET_PREC (c);
58
59 an = MPFR_PREC2LIMBS (aq); /* number of limbs of a */
60 aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS;
61 sh = aq2 - aq; /* non-significant bits in low limb */
62
63 bn = MPFR_PREC2LIMBS (bq); /* number of limbs of b */
64 cn = MPFR_PREC2LIMBS (cq); /* number of limbs of c */
65
66 ap = MPFR_MANT(a);
67 bp = MPFR_MANT(b);
68 cp = MPFR_MANT(c);
69
70 if (MPFR_UNLIKELY(ap == bp))
71 {
72 bp = MPFR_TMP_LIMBS_ALLOC (bn);
73 MPN_COPY (bp, ap, bn);
74 if (ap == cp)
75 { cp = bp; }
76 }
77 else if (ap == cp)
78 {
79 cp = MPFR_TMP_LIMBS_ALLOC (cn);
80 MPN_COPY(cp, ap, cn);
81 }
82
83 MPFR_SET_SAME_SIGN(a, b);
84 MPFR_UPDATE2_RND_MODE (rnd_mode, MPFR_SIGN (b));
85 /* now rnd_mode is either MPFR_RNDN, MPFR_RNDZ, MPFR_RNDA or MPFR_RNDF. */
86 if (MPFR_UNLIKELY (MPFR_IS_UBF (c)))
87 {
88 MPFR_STAT_STATIC_ASSERT (MPFR_EXP_MAX > MPFR_PREC_MAX);
89 diff_exp = mpfr_ubf_diff_exp (b, c);
90 }
91 else
92 diff_exp = exp - MPFR_GET_EXP (c);
93
94 MPFR_ASSERTD (diff_exp >= 0);
95
96 /*
97 * 1. Compute the significant part A', the non-significant bits of A
98 * are taken into account.
99 *
100 * 2. Perform the rounding. At each iteration, we remember:
101 * _ r = rounding bit
102 * _ f = following bits (same value)
103 * where the result has the form: [number A]rfff...fff + a remaining
104 * value in the interval [0,2) ulp. We consider the most significant
105 * bits of the remaining value to update the result; a possible carry
106 * is immediately taken into account and A is updated accordingly. As
107 * soon as the bits f don't have the same value, A can be rounded.
108 * Variables:
109 * _ rb = rounding bit (0 or 1).
110 * _ fb = following bits (0 or 1), then sticky bit.
111 * If fb == 0, the only thing that can change is the sticky bit.
112 */
113
114 rb = fb = -1; /* means: not initialized */
115
116 if (MPFR_UNLIKELY (MPFR_UEXP (aq2) <= diff_exp))
117 { /* c does not overlap with a' */
118 if (MPFR_UNLIKELY(an > bn))
119 { /* a has more limbs than b */
120 /* copy b to the most significant limbs of a */
121 MPN_COPY(ap + (an - bn), bp, bn);
122 /* zero the least significant limbs of a */
123 MPN_ZERO(ap, an - bn);
124 }
125 else /* an <= bn */
126 {
127 /* copy the most significant limbs of b to a */
128 MPN_COPY(ap, bp + (bn - an), an);
129 }
130 }
131 else /* aq2 > diff_exp */
132 { /* c overlaps with a' */
133 mp_limb_t *a2p;
134 mp_limb_t cc;
135 mpfr_prec_t dif;
136 mp_size_t difn, k;
137 int shift;
138
139 /* copy c (shifted) into a */
140
141 dif = aq2 - diff_exp;
142 /* dif is the number of bits of c which overlap with a' */
143
144 difn = MPFR_PREC2LIMBS (dif);
145 /* only the highest difn limbs from c have to be considered */
146 if (MPFR_UNLIKELY(difn > cn))
147 {
148 /* c doesn't have enough limbs; take into account the virtual
149 zero limbs now by zeroing the least significant limbs of a' */
150 MPFR_ASSERTD(difn - cn <= an);
151 MPN_ZERO(ap, difn - cn);
152 difn = cn;
153 }
154 k = diff_exp / GMP_NUMB_BITS;
155
156 /* zero the most significant k limbs of a */
157 a2p = ap + (an - k);
158 MPN_ZERO(a2p, k);
159
160 shift = diff_exp % GMP_NUMB_BITS;
161
162 if (MPFR_LIKELY(shift))
163 {
164 MPFR_ASSERTD(a2p - difn >= ap);
165 cc = mpn_rshift(a2p - difn, cp + (cn - difn), difn, shift);
166 if (MPFR_UNLIKELY(a2p - difn > ap))
167 *(a2p - difn - 1) = cc;
168 }
169 else
170 MPN_COPY(a2p - difn, cp + (cn - difn), difn);
171
172 /* add b to a */
173 cc = an > bn
174 ? mpn_add_n(ap + (an - bn), ap + (an - bn), bp, bn)
175 : mpn_add_n(ap, ap, bp + (bn - an), an);
176
177 if (MPFR_UNLIKELY(cc)) /* carry */
178 {
179 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
180 {
181 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
182 goto end_of_add;
183 }
184 exp++;
185 rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
186 if (MPFR_LIKELY(sh))
187 {
188 mp_limb_t mask, bb;
189
190 mask = MPFR_LIMB_MASK (sh);
191 bb = ap[0] & mask;
192 ap[0] &= MPFR_LIMB_LSHIFT (~mask, 1);
193 if (bb == 0)
194 fb = 0;
195 else if (bb == mask)
196 fb = 1;
197 }
198 mpn_rshift(ap, ap, an, 1);
199 ap[an-1] += MPFR_LIMB_HIGHBIT;
200 if (sh && fb < 0)
201 goto rounding;
202 } /* cc */
203 } /* aq2 > diff_exp */
204
205 /* zero the non-significant bits of a */
206 if (MPFR_LIKELY(rb < 0 && sh))
207 {
208 mp_limb_t mask, bb;
209
210 mask = MPFR_LIMB_MASK (sh);
211 bb = ap[0] & mask;
212 ap[0] &= ~mask;
213 rb = bb >> (sh - 1);
214 if (MPFR_LIKELY(sh > 1))
215 {
216 mask >>= 1;
217 bb &= mask;
218 if (bb == 0)
219 fb = 0;
220 else if (bb == mask)
221 fb = 1;
222 else
223 goto rounding;
224 }
225 }
226
227 /* Determine rounding and sticky bits (and possible carry).
228 In faithful rounding, we may stop two bits after ulp(a):
229 the approximation is regarded as the number formed by a,
230 the rounding bit rb and an additional bit fb; and the
231 corresponding error is < 1/2 ulp of the unrounded result. */
232
233 difw = (mpfr_exp_t) an - (mpfr_exp_t) (diff_exp / GMP_NUMB_BITS);
234 /* difw is the number of limbs from b (regarded as having an infinite
235 precision) that have already been combined with c; -n if the next
236 n limbs from b won't be combined with c. */
237
238 if (MPFR_UNLIKELY(bn > an))
239 { /* there are still limbs from b that haven't been taken into account */
240 mp_size_t bk;
241
242 if (fb == 0 && difw <= 0)
243 {
244 fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
245 goto rounding;
246 }
247
248 bk = bn - an; /* index of lowest considered limb from b, > 0 */
249 while (difw < 0)
250 { /* ulp(next limb from b) > msb(c) */
251 mp_limb_t bb;
252
253 bb = bp[--bk];
254
255 MPFR_ASSERTD(fb != 0);
256 if (fb > 0)
257 {
258 /* Note: Here, we can round to nearest, but the loop may still
259 be necessary to determine whether there is a carry from c,
260 which will have an effect on the ternary value. However, in
261 faithful rounding, we do not have to determine the ternary
262 value, so that we can end the loop here. */
263 if (bb != MPFR_LIMB_MAX || rnd_mode == MPFR_RNDF)
264 goto rounding;
265 }
266 else /* fb not initialized yet */
267 {
268 if (rb < 0) /* rb not initialized yet */
269 {
270 rb = bb >> (GMP_NUMB_BITS - 1);
271 bb |= MPFR_LIMB_HIGHBIT;
272 }
273 fb = 1;
274 if (bb != MPFR_LIMB_MAX)
275 goto rounding;
276 }
277
278 if (bk == 0)
279 { /* b has entirely been read */
280 fb = 1; /* c hasn't been taken into account
281 ==> sticky bit != 0 */
282 goto rounding;
283 }
284
285 difw++;
286 } /* while */
287 MPFR_ASSERTD(bk > 0 && difw >= 0);
288
289 if (difw <= cn)
290 {
291 mp_size_t ck;
292 mp_limb_t cprev;
293 int difs;
294
295 ck = cn - difw;
296 difs = diff_exp % GMP_NUMB_BITS;
297
298 if (difs == 0 && ck == 0)
299 goto c_read;
300
301 cprev = ck == cn ? 0 : cp[ck];
302
303 if (fb < 0)
304 {
305 mp_limb_t bb, cc;
306
307 if (difs)
308 {
309 cc = cprev << (GMP_NUMB_BITS - difs);
310 if (--ck >= 0)
311 {
312 cprev = cp[ck];
313 cc += cprev >> difs;
314 }
315 }
316 else
317 cc = cp[--ck];
318
319 bb = bp[--bk] + cc;
320
321 if (bb < cc /* carry */
322 && (rb < 0 || (rb ^= 1) == 0)
323 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
324 {
325 if (exp == __gmpfr_emax)
326 {
327 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
328 goto end_of_add;
329 }
330 exp++;
331 ap[an-1] = MPFR_LIMB_HIGHBIT;
332 rb = 0;
333 }
334
335 if (rb < 0) /* rb not initialized yet */
336 {
337 rb = bb >> (GMP_NUMB_BITS - 1);
338 bb <<= 1;
339 bb |= bb >> (GMP_NUMB_BITS - 1);
340 }
341
342 fb = bb != 0;
343 if (fb && bb != MPFR_LIMB_MAX)
344 goto rounding;
345 } /* fb < 0 */
346
347 /* At least two bits after ulp(a) have been read, which is
348 sufficient for faithful rounding, as we do not need to
349 determine on which side of a breakpoint the result is. */
350 if (rnd_mode == MPFR_RNDF)
351 goto rounding;
352
353 while (bk > 0)
354 {
355 mp_limb_t bb, cc;
356
357 if (difs)
358 {
359 if (ck < 0)
360 goto c_read;
361 cc = cprev << (GMP_NUMB_BITS - difs);
362 if (--ck >= 0)
363 {
364 cprev = cp[ck];
365 cc += cprev >> difs;
366 }
367 }
368 else
369 {
370 if (ck == 0)
371 goto c_read;
372 cc = cp[--ck];
373 }
374
375 bb = bp[--bk] + cc;
376 if (bb < cc) /* carry */
377 {
378 fb ^= 1;
379 if (fb)
380 goto rounding;
381 rb ^= 1;
382 if (rb == 0 && mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh))
383 {
384 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
385 {
386 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
387 goto end_of_add;
388 }
389 exp++;
390 ap[an-1] = MPFR_LIMB_HIGHBIT;
391 }
392 } /* bb < cc */
393
394 if (!fb && bb != 0)
395 {
396 fb = 1;
397 goto rounding;
398 }
399 if (fb && bb != MPFR_LIMB_MAX)
400 goto rounding;
401 } /* while */
402
403 /* b has entirely been read */
404
405 if (fb || ck < 0)
406 goto rounding;
407 if (difs && MPFR_LIMB_LSHIFT(cprev, GMP_NUMB_BITS - difs) != 0)
408 {
409 fb = 1;
410 goto rounding;
411 }
412 while (ck)
413 {
414 if (cp[--ck])
415 {
416 fb = 1;
417 goto rounding;
418 }
419 } /* while */
420 } /* difw <= cn */
421 else
422 { /* c has entirely been read */
423 c_read:
424 if (fb < 0) /* fb not initialized yet */
425 {
426 mp_limb_t bb;
427
428 MPFR_ASSERTD(bk > 0);
429 bb = bp[--bk];
430 if (rb < 0) /* rb not initialized yet */
431 {
432 rb = bb >> (GMP_NUMB_BITS - 1);
433 bb &= ~MPFR_LIMB_HIGHBIT;
434 }
435 fb = bb != 0;
436 } /* fb < 0 */
437 if (fb || rnd_mode == MPFR_RNDF)
438 goto rounding;
439 while (bk)
440 {
441 if (bp[--bk])
442 {
443 fb = 1;
444 goto rounding;
445 }
446 } /* while */
447 } /* difw > cn */
448 } /* bn > an */
449 else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
450 { /* b has entirely been read */
451 if (difw > cn)
452 { /* c has entirely been read */
453 if (rb < 0)
454 rb = 0;
455 fb = 0;
456 }
457 else if (diff_exp > MPFR_UEXP (aq2))
458 { /* b is followed by at least a zero bit, then by c */
459 if (rb < 0)
460 rb = 0;
461 fb = 1;
462 }
463 else
464 {
465 mp_size_t ck;
466 int difs;
467
468 MPFR_ASSERTD(difw >= 0 && cn >= difw);
469 ck = cn - difw;
470 difs = diff_exp % GMP_NUMB_BITS;
471
472 if (difs == 0 && ck == 0)
473 { /* c has entirely been read */
474 if (rb < 0)
475 rb = 0;
476 fb = 0;
477 }
478 else
479 {
480 mp_limb_t cc;
481
482 cc = difs ? (MPFR_ASSERTD(ck < cn),
483 cp[ck] << (GMP_NUMB_BITS - difs)) : cp[--ck];
484 if (rb < 0)
485 {
486 rb = cc >> (GMP_NUMB_BITS - 1);
487 cc &= ~MPFR_LIMB_HIGHBIT;
488 }
489 if (cc == 0 && rnd_mode == MPFR_RNDF)
490 {
491 fb = 0;
492 goto rounding;
493 }
494 while (cc == 0)
495 {
496 if (ck == 0)
497 {
498 fb = 0;
499 goto rounding;
500 }
501 cc = cp[--ck];
502 } /* while */
503 fb = 1;
504 }
505 }
506 } /* fb != 1 */
507
508 rounding:
509 /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDF, MPFR_RNDZ or MPFR_RNDA */
510 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDF))
511 {
512 if (fb == 0)
513 {
514 if (rb == 0)
515 {
516 inex = 0;
517 goto set_exponent;
518 }
519 /* round to even */
520 if (ap[0] & (MPFR_LIMB_ONE << sh))
521 goto rndn_away;
522 else
523 goto rndn_zero;
524 }
525 if (rb == 0)
526 {
527 rndn_zero:
528 inex = MPFR_IS_NEG(a) ? 1 : -1;
529 goto set_exponent;
530 }
531 else
532 {
533 rndn_away:
534 inex = MPFR_IS_POS(a) ? 1 : -1;
535 goto add_one_ulp;
536 }
537 }
538 else if (rnd_mode == MPFR_RNDZ)
539 {
540 inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
541 goto set_exponent;
542 }
543 else
544 {
545 MPFR_ASSERTN (rnd_mode == MPFR_RNDA);
546 inex = rb || fb ? (MPFR_IS_POS(a) ? 1 : -1) : 0;
547 if (inex)
548 goto add_one_ulp;
549 else
550 goto set_exponent;
551 }
552
553 add_one_ulp: /* add one unit in last place to a */
554 if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
555 {
556 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
557 {
558 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
559 goto end_of_add;
560 }
561 exp++;
562 ap[an-1] = MPFR_LIMB_HIGHBIT;
563 }
564
565 set_exponent:
566 if (MPFR_UNLIKELY (exp < __gmpfr_emin)) /* possible if b and c are UBF's */
567 {
568 if (rnd_mode == MPFR_RNDN &&
569 (exp < __gmpfr_emin - 1 ||
570 (inex >= 0 && mpfr_powerof2_raw (a))))
571 rnd_mode = MPFR_RNDZ;
572 inex = mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
573 goto end_of_add;
574 }
575 MPFR_SET_EXP (a, exp);
576
577 end_of_add:
578 MPFR_TMP_FREE(marker);
579 MPFR_RET (inex);
580 }
581