1 /* mpfr_add1 -- internal function to perform a "real" addition
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #include "mpfr-impl.h"
25 /* compute sign(b) * (|b| + |c|) */
27 mpfr_add1 (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mp_rnd_t rnd_mode
)
29 mp_limb_t
*ap
, *bp
, *cp
;
30 mp_prec_t aq
, bq
, cq
, aq2
;
34 mp_exp_unsigned_t diff_exp
;
35 MPFR_TMP_DECL(marker
);
37 MPFR_ASSERTD(MPFR_IS_PURE_FP(b
) && MPFR_IS_PURE_FP(c
));
39 MPFR_TMP_MARK(marker
);
45 an
= (aq
-1)/BITS_PER_MP_LIMB
+1; /* number of limbs of a */
46 aq2
= (mp_prec_t
) an
* BITS_PER_MP_LIMB
;
47 sh
= aq2
- aq
; /* non-significant bits in low limb */
49 bn
= (bq
-1)/BITS_PER_MP_LIMB
+1; /* number of limbs of b */
50 cn
= (cq
-1)/BITS_PER_MP_LIMB
+1; /* number of limbs of c */
56 if (MPFR_UNLIKELY(ap
== bp
))
58 bp
= (mp_ptr
) MPFR_TMP_ALLOC (bn
* BYTES_PER_MP_LIMB
);
59 MPN_COPY (bp
, ap
, bn
);
63 else if (MPFR_UNLIKELY(ap
== cp
))
65 cp
= (mp_ptr
) MPFR_TMP_ALLOC (cn
* BYTES_PER_MP_LIMB
);
69 exp
= MPFR_GET_EXP (b
);
70 MPFR_SET_SAME_SIGN(a
, b
);
71 diff_exp
= (mp_exp_unsigned_t
) exp
- MPFR_GET_EXP(c
);
74 * 1. Compute the significant part A', the non-significant bits of A
75 * are taken into account.
77 * 2. Perform the rounding. At each iteration, we remember:
79 * _ f = following bits (same value)
80 * where the result has the form: [number A]rfff...fff + a remaining
81 * value in the interval [0,2) ulp. We consider the most significant
82 * bits of the remaining value to update the result; a possible carry
83 * is immediately taken into account and A is updated accordingly. As
84 * soon as the bits f don't have the same value, A can be rounded.
86 * _ rb = rounding bit (0 or 1).
87 * _ fb = following bits (0 or 1), then sticky bit.
88 * If fb == 0, the only thing that can change is the sticky bit.
91 rb
= fb
= -1; /* means: not initialized */
93 if (MPFR_UNLIKELY(aq2
<= diff_exp
))
94 { /* c does not overlap with a' */
95 if (MPFR_UNLIKELY(an
> bn
))
96 { /* a has more limbs than b */
97 /* copy b to the most significant limbs of a */
98 MPN_COPY(ap
+ (an
- bn
), bp
, bn
);
99 /* zero the least significant limbs of a */
100 MPN_ZERO(ap
, an
- bn
);
104 /* copy the most significant limbs of b to a */
105 MPN_COPY(ap
, bp
+ (bn
- an
), an
);
108 else /* aq2 > diff_exp */
109 { /* c overlaps with a' */
116 /* copy c (shifted) into a */
118 dif
= aq2
- diff_exp
;
119 /* dif is the number of bits of c which overlap with a' */
121 difn
= (dif
-1)/BITS_PER_MP_LIMB
+ 1;
122 /* only the highest difn limbs from c have to be considered */
123 if (MPFR_UNLIKELY(difn
> cn
))
125 /* c doesn't have enough limbs; take into account the virtual
126 zero limbs now by zeroing the least significant limbs of a' */
127 MPFR_ASSERTD(difn
- cn
<= an
);
128 MPN_ZERO(ap
, difn
- cn
);
131 k
= diff_exp
/ BITS_PER_MP_LIMB
;
133 /* zero the most significant k limbs of a */
137 shift
= diff_exp
% BITS_PER_MP_LIMB
;
139 if (MPFR_LIKELY(shift
))
141 MPFR_ASSERTD(a2p
- difn
>= ap
);
142 cc
= mpn_rshift(a2p
- difn
, cp
+ (cn
- difn
), difn
, shift
);
143 if (MPFR_UNLIKELY(a2p
- difn
> ap
))
144 *(a2p
- difn
- 1) = cc
;
147 MPN_COPY(a2p
- difn
, cp
+ (cn
- difn
), difn
);
150 cc
= MPFR_UNLIKELY(an
> bn
)
151 ? mpn_add_n(ap
+ (an
- bn
), ap
+ (an
- bn
), bp
, bn
)
152 : mpn_add_n(ap
, ap
, bp
+ (bn
- an
), an
);
154 if (MPFR_UNLIKELY(cc
)) /* carry */
156 if (MPFR_UNLIKELY(exp
== __gmpfr_emax
))
158 inex
= mpfr_overflow(a
, rnd_mode
, MPFR_SIGN(a
));
162 rb
= (ap
[0] >> sh
) & 1; /* LSB(a) --> rounding bit after the shift */
167 mask
= MPFR_LIMB_MASK (sh
);
169 ap
[0] &= (~mask
) << 1;
175 mpn_rshift(ap
, ap
, an
, 1);
176 ap
[an
-1] += MPFR_LIMB_HIGHBIT
;
180 } /* aq2 > diff_exp */
182 /* non-significant bits of a */
183 if (MPFR_LIKELY(rb
< 0 && sh
))
187 mask
= MPFR_LIMB_MASK (sh
);
191 if (MPFR_LIKELY(sh
> 1))
204 /* determine rounding and sticky bits (and possible carry) */
206 difw
= (mp_exp_t
) an
- (mp_exp_t
) (diff_exp
/ BITS_PER_MP_LIMB
);
207 /* difw is the number of limbs from b (regarded as having an infinite
208 precision) that have already been combined with c; -n if the next
209 n limbs from b won't be combined with c. */
211 if (MPFR_UNLIKELY(bn
> an
))
212 { /* there are still limbs from b that haven't been taken into account */
215 if (fb
== 0 && difw
<= 0)
217 fb
= 1; /* c hasn't been taken into account ==> sticky bit != 0 */
221 bk
= bn
- an
; /* index of lowest considered limb from b, > 0 */
223 { /* ulp(next limb from b) > msb(c) */
228 MPFR_ASSERTD(fb
!= 0);
231 if (bb
!= MP_LIMB_T_MAX
)
233 fb
= 1; /* c hasn't been taken into account
234 ==> sticky bit != 0 */
238 else /* fb not initialized yet */
240 if (rb
< 0) /* rb not initialized yet */
242 rb
= bb
>> (BITS_PER_MP_LIMB
- 1);
243 bb
|= MPFR_LIMB_HIGHBIT
;
246 if (bb
!= MP_LIMB_T_MAX
)
251 { /* b has entirely been read */
252 fb
= 1; /* c hasn't been taken into account
253 ==> sticky bit != 0 */
259 MPFR_ASSERTD(bk
> 0 && difw
>= 0);
268 difs
= diff_exp
% BITS_PER_MP_LIMB
;
270 if (difs
== 0 && ck
== 0)
273 cprev
= ck
== cn
? 0 : cp
[ck
];
281 cc
= cprev
<< (BITS_PER_MP_LIMB
- difs
);
293 if (bb
< cc
/* carry */
294 && (rb
< 0 || (rb
^= 1) == 0)
295 && mpn_add_1(ap
, ap
, an
, MPFR_LIMB_ONE
<< sh
))
297 if (exp
== __gmpfr_emax
)
299 inex
= mpfr_overflow(a
, rnd_mode
, MPFR_SIGN(a
));
303 ap
[an
-1] = MPFR_LIMB_HIGHBIT
;
307 if (rb
< 0) /* rb not initialized yet */
309 rb
= bb
>> (BITS_PER_MP_LIMB
- 1);
311 bb
|= bb
>> (BITS_PER_MP_LIMB
- 1);
315 if (fb
&& bb
!= MP_LIMB_T_MAX
)
327 cc
= cprev
<< (BITS_PER_MP_LIMB
- difs
);
342 if (bb
< cc
) /* carry */
348 if (rb
== 0 && mpn_add_1(ap
, ap
, an
, MPFR_LIMB_ONE
<< sh
))
350 if (MPFR_UNLIKELY(exp
== __gmpfr_emax
))
352 inex
= mpfr_overflow(a
, rnd_mode
, MPFR_SIGN(a
));
356 ap
[an
-1] = MPFR_LIMB_HIGHBIT
;
365 if (fb
&& bb
!= MP_LIMB_T_MAX
)
369 /* b has entirely been read */
373 if (difs
&& cprev
<< (BITS_PER_MP_LIMB
- difs
))
388 { /* c has entirely been read */
390 if (fb
< 0) /* fb not initialized yet */
394 MPFR_ASSERTD(bk
> 0);
396 if (rb
< 0) /* rb not initialized yet */
398 rb
= bb
>> (BITS_PER_MP_LIMB
- 1);
399 bb
&= ~MPFR_LIMB_HIGHBIT
;
415 else if (fb
!= 1) /* if fb == 1, the sticky bit is 1 (no possible carry) */
416 { /* b has entirely been read */
418 { /* c has entirely been read */
423 else if (diff_exp
> aq2
)
424 { /* b is followed by at least a zero bit, then by c */
434 MPFR_ASSERTD(difw
>= 0 && cn
>= difw
);
436 difs
= diff_exp
% BITS_PER_MP_LIMB
;
438 if (difs
== 0 && ck
== 0)
439 { /* c has entirely been read */
448 cc
= difs
? (MPFR_ASSERTD(ck
< cn
),
449 cp
[ck
] << (BITS_PER_MP_LIMB
- difs
)) : cp
[--ck
];
452 rb
= cc
>> (BITS_PER_MP_LIMB
- 1);
453 cc
&= ~MPFR_LIMB_HIGHBIT
;
481 if (ap
[0] & (MPFR_LIMB_ONE
<< sh
))
489 inex
= MPFR_IS_NEG(a
) ? 1 : -1;
495 inex
= MPFR_IS_POS(a
) ? 1 : -1;
500 inex
= rb
|| fb
? (MPFR_IS_NEG(a
) ? 1 : -1) : 0;
505 if (inex
&& MPFR_IS_POS(a
))
512 if (inex
&& MPFR_IS_NEG(a
))
523 add_one_ulp
: /* add one unit in last place to a */
524 if (MPFR_UNLIKELY(mpn_add_1(ap
, ap
, an
, MPFR_LIMB_ONE
<< sh
)))
526 if (MPFR_UNLIKELY(exp
== __gmpfr_emax
))
528 inex
= mpfr_overflow(a
, rnd_mode
, MPFR_SIGN(a
));
532 ap
[an
-1] = MPFR_LIMB_HIGHBIT
;
536 MPFR_SET_EXP (a
, exp
);
539 MPFR_TMP_FREE(marker
);