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