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 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Cacao 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 MPFR_PREC (a), MPFR_PREC (b), MPFR_PREC (c));
56 mpfr_fprint_binary (stderr, tmpb);
57 fprintf (stderr, "\nC = ");
58 mpfr_fprint_binary (stderr, tmpc);
59 fprintf (stderr, "\n\nadd1 : ");
60 mpfr_fprint_binary (stderr, tmpa);
61 fprintf (stderr, "\nadd1sp: ");
62 mpfr_fprint_binary (stderr, a);
63 fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
67 mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
70 # define mpfr_add1sp mpfr_add1sp2
74 /* Debugging support */
79 # define DEBUG(x) /**/
82 /* compute sign(b) * (|b| + |c|)
83 Returns 0 iff result is exact,
84 a negative value when the result is less than the exact value,
85 a positive value otherwise. */
87 mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
97 MPFR_TMP_DECL(marker);
99 MPFR_TMP_MARK(marker);
101 MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
102 MPFR_ASSERTD(MPFR_IS_PURE_FP(b) && MPFR_IS_PURE_FP(c));
103 MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));
105 /* Read prec and num of limbs */
107 n = (p+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
108 MPFR_UNSIGNED_MINUS_MODULO(sh, p);
109 bx = MPFR_GET_EXP(b);
110 d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
112 DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d));
114 if (MPFR_UNLIKELY(d == 0))
117 DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) );
118 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
121 limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n);
122 DEBUG( mpfr_print_mant_binary("A= ", ap, p) );
123 MPFR_ASSERTD(limb != 0); /* There must be a carry */
124 limb = ap[0]; /* Get LSB (In fact, LSW) */
125 mpn_rshift(ap, ap, n, 1); /* Shift mantissa A */
126 ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
127 ap[0] &= ~MPFR_LIMB_MASK(sh); /* Clear LSB bit */
128 if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */
129 { inexact = 0; goto set_exponent; }
131 Nearest: Even Rule => truncate or add 1
133 if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
135 if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0))
136 { inexact = -1; goto set_exponent; }
140 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
141 if (rnd_mode==MPFR_RNDZ)
142 { inexact = -1; goto set_exponent; }
146 else if (MPFR_UNLIKELY (d >= p))
148 if (MPFR_LIKELY (d > p))
150 /* d > p : Copy B in A */
154 if (MPFR_LIKELY (rnd_mode==MPFR_RNDN
155 || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))))
159 MPN_COPY (ap, MPFR_MANT(b), n);
167 MPN_COPY (ap, MPFR_MANT(b), n);
173 /* d==p : Copy B in A */
175 Nearest: Even Rule if C is a power of 2, else Add 1
177 if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
179 /* Check if C was a power of 2 */
181 if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
186 } while (k>=0 && cp[k]==0);
187 if (MPFR_UNLIKELY(k<0))
188 /* Power of 2: Even rule */
189 if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0)
190 goto copy_set_exponent;
192 /* Not a Power of 2 */
193 goto copy_add_one_ulp;
195 else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))
196 goto copy_set_exponent;
198 goto copy_add_one_ulp;
204 mp_limb_t bcp, bcp1; /* Cp and C'p+1 */
206 /* General case: 1 <= d < p */
207 cp = (mp_limb_t*) MPFR_TMP_ALLOC(n * BYTES_PER_MP_LIMB);
209 /* Shift c in temporary allocated place */
214 dm = d % GMP_NUMB_BITS;
215 m = d / GMP_NUMB_BITS;
216 if (MPFR_UNLIKELY(dm == 0))
218 /* dm = 0 and m > 0: Just copy */
220 MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
223 else if (MPFR_LIKELY(m == 0))
225 /* dm >=1 and m == 0: just shift */
226 MPFR_ASSERTD(dm >= 1);
227 mpn_rshift(cp, MPFR_MANT(c), n, dm);
231 /* dm > 0 and m > 0: shift and zero */
232 mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
237 DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
238 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
239 DEBUG( mpfr_print_mant_binary("After ", cp, p) );
241 /* Compute bcp=Cp and bcp1=C'p+1 */
242 if (MPFR_LIKELY (sh > 0))
244 /* Try to compute them from C' rather than C */
245 bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
246 if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1)))
250 /* We can't compute C'p+1 from C'. Compute it from C */
251 /* Start from bit x=p-d+sh in mantissa C
252 (+sh since we have already looked sh bits in C'!) */
253 mpfr_prec_t x = p-d+sh-1;
254 if (MPFR_LIKELY(x>p))
255 /* We are already looked at all the bits of c, so C'p+1 = 0*/
259 mp_limb_t *tp = MPFR_MANT(c);
260 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
261 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
262 DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
263 (unsigned long) x, (long) kx,
264 (unsigned long) sx));
265 /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
266 if (tp[kx] & MPFR_LIMB_MASK(sx))
271 /*If sx==0, tp[kx] hasn't been checked*/
274 } while (kx>=0 && tp[kx]==0);
282 /* Compute Cp and C'p+1 from C with sh=0 */
283 mp_limb_t *tp = MPFR_MANT(c);
284 /* Start from bit x=p-d in mantissa C */
286 mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
287 mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
288 MPFR_ASSERTD(p >= d);
289 bcp = tp[kx] & (MPFR_LIMB_ONE<<sx);
290 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
291 if (tp[kx]&MPFR_LIMB_MASK(sx))
297 } while (kx>=0 && tp[kx]==0);
301 DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh,
302 (unsigned long) bcp, (unsigned long) bcp1));
304 /* Clean shifted C' */
305 mask = ~MPFR_LIMB_MASK(sh);
308 /* Add the mantissa c from b in a */
310 limb = mpn_add_n (ap, MPFR_MANT(b), cp, n);
311 DEBUG( mpfr_print_mant_binary("Add= ", ap, p) );
313 /* Check for overflow */
314 if (MPFR_UNLIKELY (limb))
316 limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */
317 mpn_rshift (ap, ap, n, 1); /* Shift mantissa*/
318 bx++; /* Fix exponent */
319 ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
320 ap[0] &= mask; /* Clear LSB bit */
321 bcp1 |= bcp; /* Recompute C'p+1 */
322 bcp = limb; /* Recompute Cp */
323 DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
324 (unsigned long) bcp, (unsigned long) bcp1));
325 DEBUG (mpfr_print_mant_binary ("Add= ", ap, p));
329 Zero: Truncate but could be exact.
330 Away: Add 1 if Cp or C'p+1 !=0
331 Nearest: Truncate but could be exact if Cp==0
334 if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
336 if (MPFR_LIKELY(bcp == 0))
337 { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; }
338 else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0)
339 { inexact = -1; goto set_exponent; }
343 MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
344 if (rnd_mode == MPFR_RNDZ)
346 inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0;
351 if (MPFR_UNLIKELY(bcp==0 && bcp1==0))
352 { inexact = 0; goto set_exponent; }
360 /* add one unit in last place to a */
361 DEBUG( printf("AddOneUlp\n") );
362 if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) ))
364 /* Case 100000x0 = 0x1111x1 + 1*/
365 DEBUG( printf("Pow of 2\n") );
367 ap[n-1] = MPFR_LIMB_HIGHBIT;
372 if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
374 DEBUG( printf("Overflow\n") );
375 MPFR_TMP_FREE(marker);
376 MPFR_SET_SAME_SIGN(a,b);
377 return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
379 MPFR_SET_EXP (a, bx);
380 MPFR_SET_SAME_SIGN(a,b);
382 MPFR_TMP_FREE(marker);
383 MPFR_RET (inexact * MPFR_INT_SIGN (a));