1 /* mpfr_add1sp -- internal function to perform a "real" addition
2 All the op must have the same precision
4 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
27 /* Check if we have to check the result of mpfr_add1sp with mpfr_add1 */
31 int mpfr_add1sp2 (mpfr_ptr
, mpfr_srcptr
, mpfr_srcptr
, mpfr_rnd_t
);
32 int mpfr_add1sp (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mpfr_rnd_t rnd_mode
)
34 mpfr_t tmpa
, tmpb
, tmpc
;
35 int inexb
, inexc
, inexact
, inexact2
;
37 mpfr_init2 (tmpa
, MPFR_PREC (a
));
38 mpfr_init2 (tmpb
, MPFR_PREC (b
));
39 mpfr_init2 (tmpc
, MPFR_PREC (c
));
41 inexb
= mpfr_set (tmpb
, b
, MPFR_RNDN
);
42 MPFR_ASSERTN (inexb
== 0);
44 inexc
= mpfr_set (tmpc
, c
, MPFR_RNDN
);
45 MPFR_ASSERTN (inexc
== 0);
47 inexact2
= mpfr_add1 (tmpa
, tmpb
, tmpc
, rnd_mode
);
48 inexact
= mpfr_add1sp2 (a
, b
, c
, rnd_mode
);
50 if (mpfr_cmp (tmpa
, a
) || inexact
!= inexact2
)
52 fprintf (stderr
, "add1 & add1sp return different values for %s\n"
53 "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
54 mpfr_print_rnd_mode (rnd_mode
),
55 (unsigned long) MPFR_PREC (a
),
56 (unsigned long) MPFR_PREC (b
),
57 (unsigned long) MPFR_PREC (c
));
58 mpfr_fprint_binary (stderr
, tmpb
);
59 fprintf (stderr
, "\nC = ");
60 mpfr_fprint_binary (stderr
, tmpc
);
61 fprintf (stderr
, "\n\nadd1 : ");
62 mpfr_fprint_binary (stderr
, tmpa
);
63 fprintf (stderr
, "\nadd1sp: ");
64 mpfr_fprint_binary (stderr
, a
);
65 fprintf (stderr
, "\nInexact sp = %d | Inexact = %d\n",
69 mpfr_clears (tmpa
, tmpb
, tmpc
, (mpfr_ptr
) 0);
72 # define mpfr_add1sp mpfr_add1sp2
76 /* Debugging support */
81 # define DEBUG(x) /**/
84 /* compute sign(b) * (|b| + |c|)
85 Returns 0 iff result is exact,
86 a negative value when the result is less than the exact value,
87 a positive value otherwise. */
89 mpfr_add1sp (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mpfr_rnd_t rnd_mode
)
99 MPFR_TMP_DECL(marker
);
101 MPFR_TMP_MARK(marker
);
103 MPFR_ASSERTD(MPFR_PREC(a
) == MPFR_PREC(b
) && MPFR_PREC(b
) == MPFR_PREC(c
));
104 MPFR_ASSERTD(MPFR_IS_PURE_FP(b
));
105 MPFR_ASSERTD(MPFR_IS_PURE_FP(c
));
106 MPFR_ASSERTD(MPFR_GET_EXP(b
) >= MPFR_GET_EXP(c
));
108 /* Read prec and num of limbs */
110 n
= MPFR_PREC2LIMBS (p
);
111 MPFR_UNSIGNED_MINUS_MODULO(sh
, p
);
112 bx
= MPFR_GET_EXP(b
);
113 d
= (mpfr_uexp_t
) (bx
- MPFR_GET_EXP(c
));
115 DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d
));
117 if (MPFR_UNLIKELY(d
== 0))
120 DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c
), p
) );
121 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b
), p
) );
124 limb
= mpn_add_n(ap
, MPFR_MANT(b
), MPFR_MANT(c
), n
);
125 DEBUG( mpfr_print_mant_binary("A= ", ap
, p
) );
126 MPFR_ASSERTD(limb
!= 0); /* There must be a carry */
127 limb
= ap
[0]; /* Get LSB (In fact, LSW) */
128 mpn_rshift(ap
, ap
, n
, 1); /* Shift mantissa A */
129 ap
[n
-1] |= MPFR_LIMB_HIGHBIT
; /* Set MSB */
130 ap
[0] &= ~MPFR_LIMB_MASK(sh
); /* Clear LSB bit */
131 if (MPFR_LIKELY((limb
&(MPFR_LIMB_ONE
<<sh
)) == 0)) /* Check exact case */
132 { inexact
= 0; goto set_exponent
; }
134 Nearest: Even Rule => truncate or add 1
136 if (MPFR_LIKELY(rnd_mode
==MPFR_RNDN
))
138 if (MPFR_LIKELY((ap
[0]&(MPFR_LIMB_ONE
<<sh
))==0))
139 { inexact
= -1; goto set_exponent
; }
143 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(b
));
144 if (rnd_mode
==MPFR_RNDZ
)
145 { inexact
= -1; goto set_exponent
; }
149 else if (MPFR_UNLIKELY (d
>= p
))
151 if (MPFR_LIKELY (d
> p
))
153 /* d > p : Copy B in A */
157 if (MPFR_LIKELY (rnd_mode
==MPFR_RNDN
158 || MPFR_IS_LIKE_RNDZ (rnd_mode
, MPFR_IS_NEG (b
))))
162 MPN_COPY (ap
, MPFR_MANT(b
), n
);
170 MPN_COPY (ap
, MPFR_MANT(b
), n
);
176 /* d==p : Copy B in A */
178 Nearest: Even Rule if C is a power of 2, else Add 1
180 if (MPFR_LIKELY(rnd_mode
==MPFR_RNDN
))
182 /* Check if C was a power of 2 */
184 if (MPFR_UNLIKELY(cp
[n
-1] == MPFR_LIMB_HIGHBIT
))
189 } while (k
>=0 && cp
[k
]==0);
190 if (MPFR_UNLIKELY(k
<0))
191 /* Power of 2: Even rule */
192 if ((MPFR_MANT (b
)[0]&(MPFR_LIMB_ONE
<<sh
))==0)
193 goto copy_set_exponent
;
195 /* Not a Power of 2 */
196 goto copy_add_one_ulp
;
198 else if (MPFR_IS_LIKE_RNDZ (rnd_mode
, MPFR_IS_NEG (b
)))
199 goto copy_set_exponent
;
201 goto copy_add_one_ulp
;
207 mp_limb_t bcp
, bcp1
; /* Cp and C'p+1 */
209 /* General case: 1 <= d < p */
210 cp
= MPFR_TMP_LIMBS_ALLOC (n
);
212 /* Shift c in temporary allocated place */
217 dm
= d
% GMP_NUMB_BITS
;
218 m
= d
/ GMP_NUMB_BITS
;
219 if (MPFR_UNLIKELY(dm
== 0))
221 /* dm = 0 and m > 0: Just copy */
223 MPN_COPY(cp
, MPFR_MANT(c
)+m
, n
-m
);
226 else if (MPFR_LIKELY(m
== 0))
228 /* dm >=1 and m == 0: just shift */
229 MPFR_ASSERTD(dm
>= 1);
230 mpn_rshift(cp
, MPFR_MANT(c
), n
, dm
);
234 /* dm > 0 and m > 0: shift and zero */
235 mpn_rshift(cp
, MPFR_MANT(c
)+m
, n
-m
, dm
);
240 DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c
), p
) );
241 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b
), p
) );
242 DEBUG( mpfr_print_mant_binary("After ", cp
, p
) );
244 /* Compute bcp=Cp and bcp1=C'p+1 */
245 if (MPFR_LIKELY (sh
> 0))
247 /* Try to compute them from C' rather than C */
248 bcp
= (cp
[0] & (MPFR_LIMB_ONE
<<(sh
-1))) ;
249 if (MPFR_LIKELY(cp
[0]&MPFR_LIMB_MASK(sh
-1)))
253 /* We can't compute C'p+1 from C'. Compute it from C */
254 /* Start from bit x=p-d+sh in mantissa C
255 (+sh since we have already looked sh bits in C'!) */
256 mpfr_prec_t x
= p
-d
+sh
-1;
257 if (MPFR_LIKELY(x
>p
))
258 /* We are already looked at all the bits of c, so C'p+1 = 0*/
262 mp_limb_t
*tp
= MPFR_MANT(c
);
263 mp_size_t kx
= n
-1 - (x
/ GMP_NUMB_BITS
);
264 mpfr_prec_t sx
= GMP_NUMB_BITS
-1-(x
%GMP_NUMB_BITS
);
265 DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
266 (unsigned long) x
, (long) kx
,
267 (unsigned long) sx
));
268 /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
269 if (tp
[kx
] & MPFR_LIMB_MASK(sx
))
274 /*If sx==0, tp[kx] hasn't been checked*/
277 } while (kx
>=0 && tp
[kx
]==0);
285 /* Compute Cp and C'p+1 from C with sh=0 */
286 mp_limb_t
*tp
= MPFR_MANT(c
);
287 /* Start from bit x=p-d in mantissa C */
289 mp_size_t kx
= n
-1 - (x
/ GMP_NUMB_BITS
);
290 mpfr_prec_t sx
= GMP_NUMB_BITS
-1-(x
%GMP_NUMB_BITS
);
291 MPFR_ASSERTD(p
>= d
);
292 bcp
= tp
[kx
] & (MPFR_LIMB_ONE
<<sx
);
293 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
294 if (tp
[kx
]&MPFR_LIMB_MASK(sx
))
300 } while (kx
>=0 && tp
[kx
]==0);
304 DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh
,
305 (unsigned long) bcp
, (unsigned long) bcp1
));
307 /* Clean shifted C' */
308 mask
= ~MPFR_LIMB_MASK(sh
);
311 /* Add the mantissa c from b in a */
313 limb
= mpn_add_n (ap
, MPFR_MANT(b
), cp
, n
);
314 DEBUG( mpfr_print_mant_binary("Add= ", ap
, p
) );
316 /* Check for overflow */
317 if (MPFR_UNLIKELY (limb
))
319 limb
= ap
[0] & (MPFR_LIMB_ONE
<<sh
); /* Get LSB */
320 mpn_rshift (ap
, ap
, n
, 1); /* Shift mantissa*/
321 bx
++; /* Fix exponent */
322 ap
[n
-1] |= MPFR_LIMB_HIGHBIT
; /* Set MSB */
323 ap
[0] &= mask
; /* Clear LSB bit */
324 bcp1
|= bcp
; /* Recompute C'p+1 */
325 bcp
= limb
; /* Recompute Cp */
326 DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
327 (unsigned long) bcp
, (unsigned long) bcp1
));
328 DEBUG (mpfr_print_mant_binary ("Add= ", ap
, p
));
332 Zero: Truncate but could be exact.
333 Away: Add 1 if Cp or C'p+1 !=0
334 Nearest: Truncate but could be exact if Cp==0
337 if (MPFR_LIKELY(rnd_mode
== MPFR_RNDN
))
339 if (MPFR_LIKELY(bcp
== 0))
340 { inexact
= MPFR_LIKELY(bcp1
) ? -1 : 0; goto set_exponent
; }
341 else if (MPFR_UNLIKELY(bcp1
==0) && (ap
[0]&(MPFR_LIMB_ONE
<<sh
))==0)
342 { inexact
= -1; goto set_exponent
; }
346 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(b
));
347 if (rnd_mode
== MPFR_RNDZ
)
349 inexact
= MPFR_LIKELY(bcp
|| bcp1
) ? -1 : 0;
354 if (MPFR_UNLIKELY(bcp
==0 && bcp1
==0))
355 { inexact
= 0; goto set_exponent
; }
363 /* add one unit in last place to a */
364 DEBUG( printf("AddOneUlp\n") );
365 if (MPFR_UNLIKELY( mpn_add_1(ap
, ap
, n
, MPFR_LIMB_ONE
<<sh
) ))
367 /* Case 100000x0 = 0x1111x1 + 1*/
368 DEBUG( printf("Pow of 2\n") );
370 ap
[n
-1] = MPFR_LIMB_HIGHBIT
;
375 if (MPFR_UNLIKELY(bx
> __gmpfr_emax
)) /* Check for overflow */
377 DEBUG( printf("Overflow\n") );
378 MPFR_TMP_FREE(marker
);
379 MPFR_SET_SAME_SIGN(a
,b
);
380 return mpfr_overflow(a
, rnd_mode
, MPFR_SIGN(a
));
382 MPFR_SET_EXP (a
, bx
);
383 MPFR_SET_SAME_SIGN(a
,b
);
385 MPFR_TMP_FREE(marker
);
386 MPFR_RET (inexact
* MPFR_INT_SIGN (a
));