1 /* mpfr_add1 -- internal function to perform a "real" addition
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #include "mpfr-impl.h"
25 /* compute sign(b) * (|b| + |c|), assuming b and c have same sign,
26 and are not NaN, Inf, nor zero. */
28 mpfr_add1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
30 mp_limb_t *ap, *bp, *cp;
31 mpfr_prec_t aq, bq, cq, aq2;
36 MPFR_TMP_DECL(marker);
38 MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
39 MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
41 MPFR_TMP_MARK(marker);
47 an = (aq-1)/GMP_NUMB_BITS+1; /* number of limbs of a */
48 aq2 = (mpfr_prec_t) an * GMP_NUMB_BITS;
49 sh = aq2 - aq; /* non-significant bits in low limb */
51 bn = (bq-1)/GMP_NUMB_BITS+1; /* number of limbs of b */
52 cn = (cq-1)/GMP_NUMB_BITS+1; /* number of limbs of c */
58 if (MPFR_UNLIKELY(ap == bp))
60 bp = MPFR_TMP_LIMBS_ALLOC (bn);
61 MPN_COPY (bp, ap, bn);
65 else if (MPFR_UNLIKELY(ap == cp))
67 cp = MPFR_TMP_LIMBS_ALLOC (cn);
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);
80 * 1. Compute the significant part A', the non-significant bits of A
81 * are taken into account.
83 * 2. Perform the rounding. At each iteration, we remember:
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.
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.
97 rb = fb = -1; /* means: not initialized */
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);
110 /* copy the most significant limbs of b to a */
111 MPN_COPY(ap, bp + (bn - an), an);
114 else /* aq2 > diff_exp */
115 { /* c overlaps with a' */
122 /* copy c (shifted) into a */
124 dif = aq2 - diff_exp;
125 /* dif is the number of bits of c which overlap with a' */
127 difn = (dif-1)/GMP_NUMB_BITS + 1;
128 /* only the highest difn limbs from c have to be considered */
129 if (MPFR_UNLIKELY(difn > cn))
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);
137 k = diff_exp / GMP_NUMB_BITS;
139 /* zero the most significant k limbs of a */
143 shift = diff_exp % GMP_NUMB_BITS;
145 if (MPFR_LIKELY(shift))
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;
153 MPN_COPY(a2p - difn, cp + (cn - difn), difn);
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);
160 if (MPFR_UNLIKELY(cc)) /* carry */
162 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
164 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
168 rb = (ap[0] >> sh) & 1; /* LSB(a) --> rounding bit after the shift */
173 mask = MPFR_LIMB_MASK (sh);
175 ap[0] &= (~mask) << 1;
181 mpn_rshift(ap, ap, an, 1);
182 ap[an-1] += MPFR_LIMB_HIGHBIT;
186 } /* aq2 > diff_exp */
188 /* non-significant bits of a */
189 if (MPFR_LIKELY(rb < 0 && sh))
193 mask = MPFR_LIMB_MASK (sh);
197 if (MPFR_LIKELY(sh > 1))
210 /* determine rounding and sticky bits (and possible carry) */
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. */
217 if (MPFR_UNLIKELY(bn > an))
218 { /* there are still limbs from b that haven't been taken into account */
221 if (fb == 0 && difw <= 0)
223 fb = 1; /* c hasn't been taken into account ==> sticky bit != 0 */
227 bk = bn - an; /* index of lowest considered limb from b, > 0 */
229 { /* ulp(next limb from b) > msb(c) */
234 MPFR_ASSERTD(fb != 0);
237 if (bb != MP_LIMB_T_MAX)
239 fb = 1; /* c hasn't been taken into account
240 ==> sticky bit != 0 */
244 else /* fb not initialized yet */
246 if (rb < 0) /* rb not initialized yet */
248 rb = bb >> (GMP_NUMB_BITS - 1);
249 bb |= MPFR_LIMB_HIGHBIT;
252 if (bb != MP_LIMB_T_MAX)
257 { /* b has entirely been read */
258 fb = 1; /* c hasn't been taken into account
259 ==> sticky bit != 0 */
265 MPFR_ASSERTD(bk > 0 && difw >= 0);
274 difs = diff_exp % GMP_NUMB_BITS;
276 if (difs == 0 && ck == 0)
279 cprev = ck == cn ? 0 : cp[ck];
287 cc = cprev << (GMP_NUMB_BITS - difs);
299 if (bb < cc /* carry */
300 && (rb < 0 || (rb ^= 1) == 0)
301 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
303 if (exp == __gmpfr_emax)
305 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
309 ap[an-1] = MPFR_LIMB_HIGHBIT;
313 if (rb < 0) /* rb not initialized yet */
315 rb = bb >> (GMP_NUMB_BITS - 1);
317 bb |= bb >> (GMP_NUMB_BITS - 1);
321 if (fb && bb != MP_LIMB_T_MAX)
333 cc = cprev << (GMP_NUMB_BITS - difs);
348 if (bb < cc) /* carry */
354 if (rb == 0 && mpn_add_1(ap, ap, an, MPFR_LIMB_ONE << sh))
356 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
358 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
362 ap[an-1] = MPFR_LIMB_HIGHBIT;
371 if (fb && bb != MP_LIMB_T_MAX)
375 /* b has entirely been read */
379 if (difs && cprev << (GMP_NUMB_BITS - difs))
394 { /* c has entirely been read */
396 if (fb < 0) /* fb not initialized yet */
400 MPFR_ASSERTD(bk > 0);
402 if (rb < 0) /* rb not initialized yet */
404 rb = bb >> (GMP_NUMB_BITS - 1);
405 bb &= ~MPFR_LIMB_HIGHBIT;
421 else if (fb != 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
422 { /* b has entirely been read */
424 { /* c has entirely been read */
429 else if (diff_exp > MPFR_UEXP (aq2))
430 { /* b is followed by at least a zero bit, then by c */
440 MPFR_ASSERTD(difw >= 0 && cn >= difw);
442 difs = diff_exp % GMP_NUMB_BITS;
444 if (difs == 0 && ck == 0)
445 { /* c has entirely been read */
454 cc = difs ? (MPFR_ASSERTD(ck < cn),
455 cp[ck] << (GMP_NUMB_BITS - difs)) : cp[--ck];
458 rb = cc >> (GMP_NUMB_BITS - 1);
459 cc &= ~MPFR_LIMB_HIGHBIT;
476 /* rnd_mode should be one of MPFR_RNDN, MPFR_RNDZ or MPFR_RNDA */
477 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
487 if (ap[0] & (MPFR_LIMB_ONE << sh))
495 inex = MPFR_IS_NEG(a) ? 1 : -1;
501 inex = MPFR_IS_POS(a) ? 1 : -1;
505 else if (rnd_mode == MPFR_RNDZ)
507 inex = rb || fb ? (MPFR_IS_NEG(a) ? 1 : -1) : 0;
512 MPFR_ASSERTN (rnd_mode == MPFR_RNDA);
513 inex = rb || fb ? (MPFR_IS_POS(a) ? 1 : -1) : 0;
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)))
523 if (MPFR_UNLIKELY(exp == __gmpfr_emax))
525 inex = mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
529 ap[an-1] = MPFR_LIMB_HIGHBIT;
533 MPFR_SET_EXP (a, exp);
536 MPFR_TMP_FREE(marker);