14a238c70SJohn Marino /* mpfr_add1sp -- internal function to perform a "real" addition
24a238c70SJohn Marino All the op must have the same precision
34a238c70SJohn Marino
4*ab6d115fSJohn Marino Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
64a238c70SJohn Marino
74a238c70SJohn Marino This file is part of the GNU MPFR Library.
84a238c70SJohn Marino
94a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
104a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
114a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
124a238c70SJohn Marino option) any later version.
134a238c70SJohn Marino
144a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
154a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
164a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
174a238c70SJohn Marino License for more details.
184a238c70SJohn Marino
194a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
204a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
214a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
224a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
234a238c70SJohn Marino
244a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H
254a238c70SJohn Marino #include "mpfr-impl.h"
264a238c70SJohn Marino
274a238c70SJohn Marino /* Check if we have to check the result of mpfr_add1sp with mpfr_add1 */
284a238c70SJohn Marino #ifdef WANT_ASSERT
294a238c70SJohn Marino # if WANT_ASSERT >= 2
304a238c70SJohn Marino
314a238c70SJohn Marino int mpfr_add1sp2 (mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t);
mpfr_add1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)324a238c70SJohn Marino int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
334a238c70SJohn Marino {
344a238c70SJohn Marino mpfr_t tmpa, tmpb, tmpc;
354a238c70SJohn Marino int inexb, inexc, inexact, inexact2;
364a238c70SJohn Marino
374a238c70SJohn Marino mpfr_init2 (tmpa, MPFR_PREC (a));
384a238c70SJohn Marino mpfr_init2 (tmpb, MPFR_PREC (b));
394a238c70SJohn Marino mpfr_init2 (tmpc, MPFR_PREC (c));
404a238c70SJohn Marino
414a238c70SJohn Marino inexb = mpfr_set (tmpb, b, MPFR_RNDN);
424a238c70SJohn Marino MPFR_ASSERTN (inexb == 0);
434a238c70SJohn Marino
444a238c70SJohn Marino inexc = mpfr_set (tmpc, c, MPFR_RNDN);
454a238c70SJohn Marino MPFR_ASSERTN (inexc == 0);
464a238c70SJohn Marino
474a238c70SJohn Marino inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode);
484a238c70SJohn Marino inexact = mpfr_add1sp2 (a, b, c, rnd_mode);
494a238c70SJohn Marino
504a238c70SJohn Marino if (mpfr_cmp (tmpa, a) || inexact != inexact2)
514a238c70SJohn Marino {
524a238c70SJohn Marino fprintf (stderr, "add1 & add1sp return different values for %s\n"
534a238c70SJohn Marino "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
544a238c70SJohn Marino mpfr_print_rnd_mode (rnd_mode),
554a238c70SJohn Marino (unsigned long) MPFR_PREC (a),
564a238c70SJohn Marino (unsigned long) MPFR_PREC (b),
574a238c70SJohn Marino (unsigned long) MPFR_PREC (c));
584a238c70SJohn Marino mpfr_fprint_binary (stderr, tmpb);
594a238c70SJohn Marino fprintf (stderr, "\nC = ");
604a238c70SJohn Marino mpfr_fprint_binary (stderr, tmpc);
614a238c70SJohn Marino fprintf (stderr, "\n\nadd1 : ");
624a238c70SJohn Marino mpfr_fprint_binary (stderr, tmpa);
634a238c70SJohn Marino fprintf (stderr, "\nadd1sp: ");
644a238c70SJohn Marino mpfr_fprint_binary (stderr, a);
654a238c70SJohn Marino fprintf (stderr, "\nInexact sp = %d | Inexact = %d\n",
664a238c70SJohn Marino inexact, inexact2);
674a238c70SJohn Marino MPFR_ASSERTN (0);
684a238c70SJohn Marino }
694a238c70SJohn Marino mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
704a238c70SJohn Marino return inexact;
714a238c70SJohn Marino }
724a238c70SJohn Marino # define mpfr_add1sp mpfr_add1sp2
734a238c70SJohn Marino # endif
744a238c70SJohn Marino #endif
754a238c70SJohn Marino
764a238c70SJohn Marino /* Debugging support */
774a238c70SJohn Marino #ifdef DEBUG
784a238c70SJohn Marino # undef DEBUG
794a238c70SJohn Marino # define DEBUG(x) (x)
804a238c70SJohn Marino #else
814a238c70SJohn Marino # define DEBUG(x) /**/
824a238c70SJohn Marino #endif
834a238c70SJohn Marino
844a238c70SJohn Marino /* compute sign(b) * (|b| + |c|)
854a238c70SJohn Marino Returns 0 iff result is exact,
864a238c70SJohn Marino a negative value when the result is less than the exact value,
874a238c70SJohn Marino a positive value otherwise. */
884a238c70SJohn Marino int
mpfr_add1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)894a238c70SJohn Marino mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
904a238c70SJohn Marino {
914a238c70SJohn Marino mpfr_uexp_t d;
924a238c70SJohn Marino mpfr_prec_t p;
934a238c70SJohn Marino unsigned int sh;
944a238c70SJohn Marino mp_size_t n;
954a238c70SJohn Marino mp_limb_t *ap, *cp;
964a238c70SJohn Marino mpfr_exp_t bx;
974a238c70SJohn Marino mp_limb_t limb;
984a238c70SJohn Marino int inexact;
994a238c70SJohn Marino MPFR_TMP_DECL(marker);
1004a238c70SJohn Marino
1014a238c70SJohn Marino MPFR_TMP_MARK(marker);
1024a238c70SJohn Marino
1034a238c70SJohn Marino MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
1044a238c70SJohn Marino MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
1054a238c70SJohn Marino MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
1064a238c70SJohn Marino MPFR_ASSERTD(MPFR_GET_EXP(b) >= MPFR_GET_EXP(c));
1074a238c70SJohn Marino
1084a238c70SJohn Marino /* Read prec and num of limbs */
1094a238c70SJohn Marino p = MPFR_PREC(b);
110*ab6d115fSJohn Marino n = MPFR_PREC2LIMBS (p);
1114a238c70SJohn Marino MPFR_UNSIGNED_MINUS_MODULO(sh, p);
1124a238c70SJohn Marino bx = MPFR_GET_EXP(b);
1134a238c70SJohn Marino d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
1144a238c70SJohn Marino
1154a238c70SJohn Marino DEBUG (printf ("New add1sp with diff=%lu\n", (unsigned long) d));
1164a238c70SJohn Marino
1174a238c70SJohn Marino if (MPFR_UNLIKELY(d == 0))
1184a238c70SJohn Marino {
1194a238c70SJohn Marino /* d==0 */
1204a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("C= ", MPFR_MANT(c), p) );
1214a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
1224a238c70SJohn Marino bx++; /* exp + 1 */
1234a238c70SJohn Marino ap = MPFR_MANT(a);
1244a238c70SJohn Marino limb = mpn_add_n(ap, MPFR_MANT(b), MPFR_MANT(c), n);
1254a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("A= ", ap, p) );
1264a238c70SJohn Marino MPFR_ASSERTD(limb != 0); /* There must be a carry */
1274a238c70SJohn Marino limb = ap[0]; /* Get LSB (In fact, LSW) */
1284a238c70SJohn Marino mpn_rshift(ap, ap, n, 1); /* Shift mantissa A */
1294a238c70SJohn Marino ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
1304a238c70SJohn Marino ap[0] &= ~MPFR_LIMB_MASK(sh); /* Clear LSB bit */
1314a238c70SJohn Marino if (MPFR_LIKELY((limb&(MPFR_LIMB_ONE<<sh)) == 0)) /* Check exact case */
1324a238c70SJohn Marino { inexact = 0; goto set_exponent; }
1334a238c70SJohn Marino /* Zero: Truncate
1344a238c70SJohn Marino Nearest: Even Rule => truncate or add 1
1354a238c70SJohn Marino Away: Add 1 */
1364a238c70SJohn Marino if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
1374a238c70SJohn Marino {
1384a238c70SJohn Marino if (MPFR_LIKELY((ap[0]&(MPFR_LIMB_ONE<<sh))==0))
1394a238c70SJohn Marino { inexact = -1; goto set_exponent; }
1404a238c70SJohn Marino else
1414a238c70SJohn Marino goto add_one_ulp;
1424a238c70SJohn Marino }
1434a238c70SJohn Marino MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
1444a238c70SJohn Marino if (rnd_mode==MPFR_RNDZ)
1454a238c70SJohn Marino { inexact = -1; goto set_exponent; }
1464a238c70SJohn Marino else
1474a238c70SJohn Marino goto add_one_ulp;
1484a238c70SJohn Marino }
1494a238c70SJohn Marino else if (MPFR_UNLIKELY (d >= p))
1504a238c70SJohn Marino {
1514a238c70SJohn Marino if (MPFR_LIKELY (d > p))
1524a238c70SJohn Marino {
1534a238c70SJohn Marino /* d > p : Copy B in A */
1544a238c70SJohn Marino /* Away: Add 1
1554a238c70SJohn Marino Nearest: Trunc
1564a238c70SJohn Marino Zero: Trunc */
1574a238c70SJohn Marino if (MPFR_LIKELY (rnd_mode==MPFR_RNDN
1584a238c70SJohn Marino || MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b))))
1594a238c70SJohn Marino {
1604a238c70SJohn Marino copy_set_exponent:
1614a238c70SJohn Marino ap = MPFR_MANT (a);
1624a238c70SJohn Marino MPN_COPY (ap, MPFR_MANT(b), n);
1634a238c70SJohn Marino inexact = -1;
1644a238c70SJohn Marino goto set_exponent;
1654a238c70SJohn Marino }
1664a238c70SJohn Marino else
1674a238c70SJohn Marino {
1684a238c70SJohn Marino copy_add_one_ulp:
1694a238c70SJohn Marino ap = MPFR_MANT(a);
1704a238c70SJohn Marino MPN_COPY (ap, MPFR_MANT(b), n);
1714a238c70SJohn Marino goto add_one_ulp;
1724a238c70SJohn Marino }
1734a238c70SJohn Marino }
1744a238c70SJohn Marino else
1754a238c70SJohn Marino {
1764a238c70SJohn Marino /* d==p : Copy B in A */
1774a238c70SJohn Marino /* Away: Add 1
1784a238c70SJohn Marino Nearest: Even Rule if C is a power of 2, else Add 1
1794a238c70SJohn Marino Zero: Trunc */
1804a238c70SJohn Marino if (MPFR_LIKELY(rnd_mode==MPFR_RNDN))
1814a238c70SJohn Marino {
1824a238c70SJohn Marino /* Check if C was a power of 2 */
1834a238c70SJohn Marino cp = MPFR_MANT(c);
1844a238c70SJohn Marino if (MPFR_UNLIKELY(cp[n-1] == MPFR_LIMB_HIGHBIT))
1854a238c70SJohn Marino {
1864a238c70SJohn Marino mp_size_t k = n-1;
1874a238c70SJohn Marino do {
1884a238c70SJohn Marino k--;
1894a238c70SJohn Marino } while (k>=0 && cp[k]==0);
1904a238c70SJohn Marino if (MPFR_UNLIKELY(k<0))
1914a238c70SJohn Marino /* Power of 2: Even rule */
1924a238c70SJohn Marino if ((MPFR_MANT (b)[0]&(MPFR_LIMB_ONE<<sh))==0)
1934a238c70SJohn Marino goto copy_set_exponent;
1944a238c70SJohn Marino }
1954a238c70SJohn Marino /* Not a Power of 2 */
1964a238c70SJohn Marino goto copy_add_one_ulp;
1974a238c70SJohn Marino }
1984a238c70SJohn Marino else if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (b)))
1994a238c70SJohn Marino goto copy_set_exponent;
2004a238c70SJohn Marino else
2014a238c70SJohn Marino goto copy_add_one_ulp;
2024a238c70SJohn Marino }
2034a238c70SJohn Marino }
2044a238c70SJohn Marino else
2054a238c70SJohn Marino {
2064a238c70SJohn Marino mp_limb_t mask;
2074a238c70SJohn Marino mp_limb_t bcp, bcp1; /* Cp and C'p+1 */
2084a238c70SJohn Marino
2094a238c70SJohn Marino /* General case: 1 <= d < p */
2104a238c70SJohn Marino cp = MPFR_TMP_LIMBS_ALLOC (n);
2114a238c70SJohn Marino
2124a238c70SJohn Marino /* Shift c in temporary allocated place */
2134a238c70SJohn Marino {
2144a238c70SJohn Marino mpfr_uexp_t dm;
2154a238c70SJohn Marino mp_size_t m;
2164a238c70SJohn Marino
2174a238c70SJohn Marino dm = d % GMP_NUMB_BITS;
2184a238c70SJohn Marino m = d / GMP_NUMB_BITS;
2194a238c70SJohn Marino if (MPFR_UNLIKELY(dm == 0))
2204a238c70SJohn Marino {
2214a238c70SJohn Marino /* dm = 0 and m > 0: Just copy */
2224a238c70SJohn Marino MPFR_ASSERTD(m!=0);
2234a238c70SJohn Marino MPN_COPY(cp, MPFR_MANT(c)+m, n-m);
2244a238c70SJohn Marino MPN_ZERO(cp+n-m, m);
2254a238c70SJohn Marino }
2264a238c70SJohn Marino else if (MPFR_LIKELY(m == 0))
2274a238c70SJohn Marino {
2284a238c70SJohn Marino /* dm >=1 and m == 0: just shift */
2294a238c70SJohn Marino MPFR_ASSERTD(dm >= 1);
2304a238c70SJohn Marino mpn_rshift(cp, MPFR_MANT(c), n, dm);
2314a238c70SJohn Marino }
2324a238c70SJohn Marino else
2334a238c70SJohn Marino {
2344a238c70SJohn Marino /* dm > 0 and m > 0: shift and zero */
2354a238c70SJohn Marino mpn_rshift(cp, MPFR_MANT(c)+m, n-m, dm);
2364a238c70SJohn Marino MPN_ZERO(cp+n-m, m);
2374a238c70SJohn Marino }
2384a238c70SJohn Marino }
2394a238c70SJohn Marino
2404a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c), p) );
2414a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b), p) );
2424a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("After ", cp, p) );
2434a238c70SJohn Marino
2444a238c70SJohn Marino /* Compute bcp=Cp and bcp1=C'p+1 */
2454a238c70SJohn Marino if (MPFR_LIKELY (sh > 0))
2464a238c70SJohn Marino {
2474a238c70SJohn Marino /* Try to compute them from C' rather than C */
2484a238c70SJohn Marino bcp = (cp[0] & (MPFR_LIMB_ONE<<(sh-1))) ;
2494a238c70SJohn Marino if (MPFR_LIKELY(cp[0]&MPFR_LIMB_MASK(sh-1)))
2504a238c70SJohn Marino bcp1 = 1;
2514a238c70SJohn Marino else
2524a238c70SJohn Marino {
2534a238c70SJohn Marino /* We can't compute C'p+1 from C'. Compute it from C */
2544a238c70SJohn Marino /* Start from bit x=p-d+sh in mantissa C
2554a238c70SJohn Marino (+sh since we have already looked sh bits in C'!) */
2564a238c70SJohn Marino mpfr_prec_t x = p-d+sh-1;
2574a238c70SJohn Marino if (MPFR_LIKELY(x>p))
2584a238c70SJohn Marino /* We are already looked at all the bits of c, so C'p+1 = 0*/
2594a238c70SJohn Marino bcp1 = 0;
2604a238c70SJohn Marino else
2614a238c70SJohn Marino {
2624a238c70SJohn Marino mp_limb_t *tp = MPFR_MANT(c);
2634a238c70SJohn Marino mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
2644a238c70SJohn Marino mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
2654a238c70SJohn Marino DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
2664a238c70SJohn Marino (unsigned long) x, (long) kx,
2674a238c70SJohn Marino (unsigned long) sx));
2684a238c70SJohn Marino /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
2694a238c70SJohn Marino if (tp[kx] & MPFR_LIMB_MASK(sx))
2704a238c70SJohn Marino bcp1 = 1;
2714a238c70SJohn Marino else
2724a238c70SJohn Marino {
2734a238c70SJohn Marino /*kx += (sx==0);*/
2744a238c70SJohn Marino /*If sx==0, tp[kx] hasn't been checked*/
2754a238c70SJohn Marino do {
2764a238c70SJohn Marino kx--;
2774a238c70SJohn Marino } while (kx>=0 && tp[kx]==0);
2784a238c70SJohn Marino bcp1 = (kx >= 0);
2794a238c70SJohn Marino }
2804a238c70SJohn Marino }
2814a238c70SJohn Marino }
2824a238c70SJohn Marino }
2834a238c70SJohn Marino else /* sh == 0 */
2844a238c70SJohn Marino {
2854a238c70SJohn Marino /* Compute Cp and C'p+1 from C with sh=0 */
2864a238c70SJohn Marino mp_limb_t *tp = MPFR_MANT(c);
2874a238c70SJohn Marino /* Start from bit x=p-d in mantissa C */
2884a238c70SJohn Marino mpfr_prec_t x = p-d;
2894a238c70SJohn Marino mp_size_t kx = n-1 - (x / GMP_NUMB_BITS);
2904a238c70SJohn Marino mpfr_prec_t sx = GMP_NUMB_BITS-1-(x%GMP_NUMB_BITS);
2914a238c70SJohn Marino MPFR_ASSERTD(p >= d);
2924a238c70SJohn Marino bcp = tp[kx] & (MPFR_LIMB_ONE<<sx);
2934a238c70SJohn Marino /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
2944a238c70SJohn Marino if (tp[kx]&MPFR_LIMB_MASK(sx))
2954a238c70SJohn Marino bcp1 = 1;
2964a238c70SJohn Marino else
2974a238c70SJohn Marino {
2984a238c70SJohn Marino do {
2994a238c70SJohn Marino kx--;
3004a238c70SJohn Marino } while (kx>=0 && tp[kx]==0);
3014a238c70SJohn Marino bcp1 = (kx>=0);
3024a238c70SJohn Marino }
3034a238c70SJohn Marino }
3044a238c70SJohn Marino DEBUG (printf("sh=%u Cp=%lu C'p+1=%lu\n", sh,
3054a238c70SJohn Marino (unsigned long) bcp, (unsigned long) bcp1));
3064a238c70SJohn Marino
3074a238c70SJohn Marino /* Clean shifted C' */
3084a238c70SJohn Marino mask = ~MPFR_LIMB_MASK(sh);
3094a238c70SJohn Marino cp[0] &= mask;
3104a238c70SJohn Marino
3114a238c70SJohn Marino /* Add the mantissa c from b in a */
3124a238c70SJohn Marino ap = MPFR_MANT(a);
3134a238c70SJohn Marino limb = mpn_add_n (ap, MPFR_MANT(b), cp, n);
3144a238c70SJohn Marino DEBUG( mpfr_print_mant_binary("Add= ", ap, p) );
3154a238c70SJohn Marino
3164a238c70SJohn Marino /* Check for overflow */
3174a238c70SJohn Marino if (MPFR_UNLIKELY (limb))
3184a238c70SJohn Marino {
3194a238c70SJohn Marino limb = ap[0] & (MPFR_LIMB_ONE<<sh); /* Get LSB */
3204a238c70SJohn Marino mpn_rshift (ap, ap, n, 1); /* Shift mantissa*/
3214a238c70SJohn Marino bx++; /* Fix exponent */
3224a238c70SJohn Marino ap[n-1] |= MPFR_LIMB_HIGHBIT; /* Set MSB */
3234a238c70SJohn Marino ap[0] &= mask; /* Clear LSB bit */
3244a238c70SJohn Marino bcp1 |= bcp; /* Recompute C'p+1 */
3254a238c70SJohn Marino bcp = limb; /* Recompute Cp */
3264a238c70SJohn Marino DEBUG (printf ("(Overflow) Cp=%lu C'p+1=%lu\n",
3274a238c70SJohn Marino (unsigned long) bcp, (unsigned long) bcp1));
3284a238c70SJohn Marino DEBUG (mpfr_print_mant_binary ("Add= ", ap, p));
3294a238c70SJohn Marino }
3304a238c70SJohn Marino
3314a238c70SJohn Marino /* Round:
3324a238c70SJohn Marino Zero: Truncate but could be exact.
3334a238c70SJohn Marino Away: Add 1 if Cp or C'p+1 !=0
3344a238c70SJohn Marino Nearest: Truncate but could be exact if Cp==0
3354a238c70SJohn Marino Add 1 if C'p+1 !=0,
3364a238c70SJohn Marino Even rule else */
3374a238c70SJohn Marino if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
3384a238c70SJohn Marino {
3394a238c70SJohn Marino if (MPFR_LIKELY(bcp == 0))
3404a238c70SJohn Marino { inexact = MPFR_LIKELY(bcp1) ? -1 : 0; goto set_exponent; }
3414a238c70SJohn Marino else if (MPFR_UNLIKELY(bcp1==0) && (ap[0]&(MPFR_LIMB_ONE<<sh))==0)
3424a238c70SJohn Marino { inexact = -1; goto set_exponent; }
3434a238c70SJohn Marino else
3444a238c70SJohn Marino goto add_one_ulp;
3454a238c70SJohn Marino }
3464a238c70SJohn Marino MPFR_UPDATE_RND_MODE(rnd_mode, MPFR_IS_NEG(b));
3474a238c70SJohn Marino if (rnd_mode == MPFR_RNDZ)
3484a238c70SJohn Marino {
3494a238c70SJohn Marino inexact = MPFR_LIKELY(bcp || bcp1) ? -1 : 0;
3504a238c70SJohn Marino goto set_exponent;
3514a238c70SJohn Marino }
3524a238c70SJohn Marino else
3534a238c70SJohn Marino {
3544a238c70SJohn Marino if (MPFR_UNLIKELY(bcp==0 && bcp1==0))
3554a238c70SJohn Marino { inexact = 0; goto set_exponent; }
3564a238c70SJohn Marino else
3574a238c70SJohn Marino goto add_one_ulp;
3584a238c70SJohn Marino }
3594a238c70SJohn Marino }
3604a238c70SJohn Marino MPFR_ASSERTN(0);
3614a238c70SJohn Marino
3624a238c70SJohn Marino add_one_ulp:
3634a238c70SJohn Marino /* add one unit in last place to a */
3644a238c70SJohn Marino DEBUG( printf("AddOneUlp\n") );
3654a238c70SJohn Marino if (MPFR_UNLIKELY( mpn_add_1(ap, ap, n, MPFR_LIMB_ONE<<sh) ))
3664a238c70SJohn Marino {
3674a238c70SJohn Marino /* Case 100000x0 = 0x1111x1 + 1*/
3684a238c70SJohn Marino DEBUG( printf("Pow of 2\n") );
3694a238c70SJohn Marino bx++;
3704a238c70SJohn Marino ap[n-1] = MPFR_LIMB_HIGHBIT;
3714a238c70SJohn Marino }
3724a238c70SJohn Marino inexact = 1;
3734a238c70SJohn Marino
3744a238c70SJohn Marino set_exponent:
3754a238c70SJohn Marino if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
3764a238c70SJohn Marino {
3774a238c70SJohn Marino DEBUG( printf("Overflow\n") );
3784a238c70SJohn Marino MPFR_TMP_FREE(marker);
3794a238c70SJohn Marino MPFR_SET_SAME_SIGN(a,b);
3804a238c70SJohn Marino return mpfr_overflow(a, rnd_mode, MPFR_SIGN(a));
3814a238c70SJohn Marino }
3824a238c70SJohn Marino MPFR_SET_EXP (a, bx);
3834a238c70SJohn Marino MPFR_SET_SAME_SIGN(a,b);
3844a238c70SJohn Marino
3854a238c70SJohn Marino MPFR_TMP_FREE(marker);
3864a238c70SJohn Marino MPFR_RET (inexact * MPFR_INT_SIGN (a));
3874a238c70SJohn Marino }
388