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, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC 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
= 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 */
51 bn
= MPFR_PREC2LIMBS (bq
); /* number of limbs of b */
52 cn
= MPFR_PREC2LIMBS (cq
); /* 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
= MPFR_PREC2LIMBS (dif
);
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
);