/* mpfr_sub1 -- internal function to perform a "real" subtraction
Copyright 2001-2017 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include "mpfr-impl.h"
/* compute sign(b) * (|b| - |c|), with |b| > |c|, diff_exp = EXP(b) - EXP(c)
Returns 0 iff result is exact,
a negative value when the result is less than the exact value,
a positive value otherwise.
*/
int
mpfr_sub1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
{
int sign;
mpfr_uexp_t diff_exp;
mpfr_prec_t cancel, cancel1;
mp_size_t cancel2, an, bn, cn, cn0;
mp_limb_t *ap, *bp, *cp;
mp_limb_t carry, bb, cc;
int inexact, shift_b, shift_c, add_exp = 0;
int cmp_low = 0; /* used for rounding to nearest: 0 if low(b) = low(c),
negative if low(b) < low(c), positive if low(b)>low(c) */
int sh, k;
MPFR_TMP_DECL(marker);
MPFR_TMP_MARK(marker);
ap = MPFR_MANT(a);
an = MPFR_LIMB_SIZE(a);
sign = mpfr_cmp2 (b, c, &cancel);
if (MPFR_UNLIKELY(sign == 0))
{
if (rnd_mode == MPFR_RNDD)
MPFR_SET_NEG (a);
else
MPFR_SET_POS (a);
MPFR_SET_ZERO (a);
MPFR_RET (0);
}
/*
* If subtraction: sign(a) = sign * sign(b)
* If addition: sign(a) = sign of the larger argument in absolute value.
*
* Both cases can be simplidied in:
* if (sign>0)
* if addition: sign(a) = sign * sign(b) = sign(b)
* if subtraction, b is greater, so sign(a) = sign(b)
* else
* if subtraction, sign(a) = - sign(b)
* if addition, sign(a) = sign(c) (since c is greater)
* But if it is an addition, sign(b) and sign(c) are opposed!
* So sign(a) = - sign(b)
*/
if (sign < 0) /* swap b and c so that |b| > |c| */
{
mpfr_srcptr t;
MPFR_SET_OPPOSITE_SIGN (a,b);
t = b; b = c; c = t;
}
else
MPFR_SET_SAME_SIGN (a,b);
/* Check if c is too small.
A more precise test is to replace 2 by
(rnd == MPFR_RNDN) + mpfr_power2_raw (b)
but it is more expensive and not very useful */
if (MPFR_UNLIKELY (MPFR_GET_EXP (c) <= MPFR_GET_EXP (b)
- (mpfr_exp_t) MAX (MPFR_PREC (a), MPFR_PREC (b)) - 2))
{
/* Remember, we can't have an exact result! */
/* A.AAAAAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
/* A = S*ABS(B) +/- ulp(a) */
MPFR_SET_EXP (a, MPFR_GET_EXP (b));
MPFR_RNDRAW_EVEN (inexact, a, MPFR_MANT (b), MPFR_PREC (b),
rnd_mode, MPFR_SIGN (a), ++ MPFR_EXP (a));
if (inexact == 0)
{
/* a = b (Exact)
But we know it isn't (Since we have to remove `c')
So if we round to Zero, we have to remove one ulp.
Otherwise the result is correctly rounded. */
/* An overflow is not possible. */
MPFR_ASSERTD (MPFR_EXP (a) <= __gmpfr_emax);
if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG (a)))
{
mpfr_nexttozero (a);
MPFR_RET (- MPFR_INT_SIGN (a));
}
MPFR_RET (MPFR_INT_SIGN (a));
}
else
{
/* A.AAAAAAAAAAAAAA
= B.BBBBBBBBBBBBBBB
- C.CCCCCCCCCCCCC */
/* It isn't exact so Prec(b) > Prec(a) and the last
Prec(b)-Prec(a) bits of `b' are not zeros.
Which means that removing c from b can't generate a carry
execpt in case of even rounding.
In all other case the result and the inexact flag should be
correct (We can't have an exact result).
In case of EVEN rounding:
1.BBBBBBBBBBBBBx10
- 1.CCCCCCCCCCCC
= 1.BBBBBBBBBBBBBx01 Rounded to Prec(b)
= 1.BBBBBBBBBBBBBx Nearest / Rounded to Prec(a)
Set gives:
1.BBBBBBBBBBBBB0 if inexact == EVEN_INEX (x == 0)
1.BBBBBBBBBBBBB1+1 if inexact == -EVEN_INEX (x == 1)
which means we get a wrong rounded result if x==1,
i.e. inexact= MPFR_EVEN_INEX */
if (MPFR_UNLIKELY (inexact == MPFR_EVEN_INEX*MPFR_INT_SIGN (a)))
{
if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
mpfr_setmax (a, __gmpfr_emax);
else
mpfr_nexttozero (a);
inexact = -MPFR_INT_SIGN (a);
}
else if (MPFR_UNLIKELY (MPFR_EXP (a) > __gmpfr_emax))
inexact = mpfr_overflow (a, rnd_mode, MPFR_SIGN (a));
MPFR_RET (inexact);
}
}
diff_exp = (mpfr_uexp_t) MPFR_GET_EXP (b) - MPFR_GET_EXP (c);
/* reserve a space to store b aligned with the result, i.e. shifted by
(-cancel) % GMP_NUMB_BITS to the right */
bn = MPFR_LIMB_SIZE (b);
MPFR_UNSIGNED_MINUS_MODULO (shift_b, cancel);
cancel1 = (cancel + shift_b) / GMP_NUMB_BITS;
/* the high cancel1 limbs from b should not be taken into account */
if (MPFR_UNLIKELY (shift_b == 0))
{
bp = MPFR_MANT(b); /* no need of an extra space */
/* Ensure ap != bp */
if (MPFR_UNLIKELY (ap == bp))
{
bp = MPFR_TMP_LIMBS_ALLOC (bn);
MPN_COPY (bp, ap, bn);
}
}
else
{
bp = MPFR_TMP_LIMBS_ALLOC (bn + 1);
bp[0] = mpn_rshift (bp + 1, MPFR_MANT(b), bn++, shift_b);
}
/* reserve a space to store c aligned with the result, i.e. shifted by
(diff_exp-cancel) % GMP_NUMB_BITS to the right */
cn = MPFR_LIMB_SIZE(c);
if ((UINT_MAX % GMP_NUMB_BITS) == (GMP_NUMB_BITS-1)
&& ((-(unsigned) 1)%GMP_NUMB_BITS > 0))
shift_c = ((mpfr_uexp_t) diff_exp - cancel) % GMP_NUMB_BITS;
else
{
shift_c = diff_exp - (cancel % GMP_NUMB_BITS);
shift_c = (shift_c + GMP_NUMB_BITS) % GMP_NUMB_BITS;
}
MPFR_ASSERTD( shift_c >= 0 && shift_c < GMP_NUMB_BITS);
if (MPFR_UNLIKELY(shift_c == 0))
{
cp = MPFR_MANT(c);
/* Ensure ap != cp */
if (ap == cp)
{
cp = MPFR_TMP_LIMBS_ALLOC (cn);
MPN_COPY(cp, ap, cn);
}
}
else
{
cp = MPFR_TMP_LIMBS_ALLOC (cn + 1);
cp[0] = mpn_rshift (cp + 1, MPFR_MANT(c), cn++, shift_c);
}
#ifdef DEBUG
printf ("rnd=%s shift_b=%d shift_c=%d diffexp=%lu\n",
mpfr_print_rnd_mode (rnd_mode), shift_b, shift_c,
(unsigned long) diff_exp);
#endif
MPFR_ASSERTD (ap != cp);
MPFR_ASSERTD (bp != cp);
/* here we have shift_c = (diff_exp - cancel) % GMP_NUMB_BITS,
0 <= shift_c < GMP_NUMB_BITS
thus we want cancel2 = ceil((cancel - diff_exp) / GMP_NUMB_BITS) */
/* Possible optimization with a C99 compiler (i.e. well-defined
integer division): if MPFR_PREC_MAX is reduced to
((mpfr_prec_t)((mpfr_uprec_t)(~(mpfr_uprec_t)0)>>1) - GMP_NUMB_BITS + 1)
and diff_exp is of type mpfr_exp_t (no need for mpfr_uexp_t, since
the sum or difference of 2 exponents must be representable, as used
by the multiplication code), then the computation of cancel2 could
be simplified to
cancel2 = (cancel - (diff_exp - shift_c)) / GMP_NUMB_BITS;
because cancel, diff_exp and shift_c are all non-negative and
these variables are signed. */
MPFR_ASSERTD (cancel >= 0);
if (cancel >= diff_exp)
/* Note that cancel is signed and will be converted to mpfr_uexp_t
(type of diff_exp) in the expression below, so that this will
work even if cancel is very large and diff_exp = 0. */
cancel2 = (cancel - diff_exp + (GMP_NUMB_BITS - 1)) / GMP_NUMB_BITS;
else
cancel2 = - (mp_size_t) ((diff_exp - cancel) / GMP_NUMB_BITS);
/* the high cancel2 limbs from b should not be taken into account */
#ifdef DEBUG
printf ("cancel=%lu cancel1=%lu cancel2=%ld\n",
(unsigned long) cancel, (unsigned long) cancel1, (long) cancel2);
#endif
/* ap[an-1] ap[0]
<----------------+-----------|---->
<----------PREC(a)----------><-sh->
cancel1
limbs bp[bn-cancel1-1]
<--...-----><----------------+-----------+----------->
cancel2
limbs cp[cn-cancel2-1] cancel2 >= 0
<--...--><----------------+----------------+---------------->
(-cancel2) cancel2 < 0
limbs <----------------+---------------->
*/
/* first part: put in ap[0..an-1] the value of high(b) - high(c),
where high(b) consists of the high an+cancel1 limbs of b,
and high(c) consists of the high an+cancel2 limbs of c.
*/
/* copy high(b) into a */
if (MPFR_LIKELY(an + (mp_size_t) cancel1 <= bn))
/* a: <----------------+-----------|---->
b: <-----------------------------------------> */
MPN_COPY (ap, bp + bn - (an + cancel1), an);
else
/* a: <----------------+-----------|---->
b: <-------------------------> */
if ((mp_size_t) cancel1 < bn) /* otherwise b does not overlap with a */
{
MPN_ZERO (ap, an + cancel1 - bn);
MPN_COPY (ap + (an + cancel1 - bn), bp, bn - cancel1);
}
else
MPN_ZERO (ap, an);
#ifdef DEBUG
printf("after copying high(b), a="); mpfr_print_binary(a); putchar('\n');
#endif
/* subtract high(c) */
if (MPFR_LIKELY(an + cancel2 > 0)) /* otherwise c does not overlap with a */
{
mp_limb_t *ap2;
if (cancel2 >= 0)
{
if (an + cancel2 <= cn)
/* a: <----------------------------->
c: <-----------------------------------------> */
mpn_sub_n (ap, ap, cp + cn - (an + cancel2), an);
else
/* a: <---------------------------->
c: <-------------------------> */
{
ap2 = ap + an + (cancel2 - cn);
if (cn > cancel2)
mpn_sub_n (ap2, ap2, cp, cn - cancel2);
}
}
else /* cancel2 < 0 */
{
mp_limb_t borrow;
if (an + cancel2 <= cn)
/* a: <----------------------------->
c: <-----------------------------> */
borrow = mpn_sub_n (ap, ap, cp + cn - (an + cancel2),
an + cancel2);
else
/* a: <---------------------------->
c: <----------------> */
{
ap2 = ap + an + cancel2 - cn;
borrow = mpn_sub_n (ap2, ap2, cp, cn);
}
ap2 = ap + an + cancel2;
mpn_sub_1 (ap2, ap2, -cancel2, borrow);
}
}
#ifdef DEBUG
printf("after subtracting high(c), a=");
mpfr_print_binary(a);
putchar('\n');
#endif
/* now perform rounding */
sh = (mpfr_prec_t) an * GMP_NUMB_BITS - MPFR_PREC(a);
/* last unused bits from a */
carry = ap[0] & MPFR_LIMB_MASK (sh);
ap[0] -= carry;
if (MPFR_LIKELY(rnd_mode == MPFR_RNDN))
{
if (MPFR_LIKELY(sh))
{
/* can decide except when carry = 2^(sh-1) [middle]
or carry = 0 [truncate, but cannot decide inexact flag] */
if (carry > (MPFR_LIMB_ONE << (sh - 1)))
goto add_one_ulp;
else if ((0 < carry) && (carry < (MPFR_LIMB_ONE << (sh - 1))))
{
inexact = -1; /* result if smaller than exact value */
goto truncate;
}
/* now carry = 2^(sh-1), in which case cmp_low=2,
or carry = 0, in which case cmp_low=0 */
cmp_low = (carry == 0) ? 0 : 2;
}
}
else /* directed rounding: set rnd_mode to RNDZ iff toward zero */
{
if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd_mode, MPFR_IS_NEG(a)))
rnd_mode = MPFR_RNDZ;
if (carry)
{
if (rnd_mode == MPFR_RNDZ)
{
inexact = -1;
goto truncate;
}
else /* round away */
goto add_one_ulp;
}
}
/* we have to consider the low (bn - (an+cancel1)) limbs from b,
and the (cn - (an+cancel2)) limbs from c. */
bn -= an + cancel1;
cn0 = cn;
cn -= an + cancel2;
#ifdef DEBUG
printf ("last sh=%d bits from a are %lu, bn=%ld, cn=%ld\n",
sh, (unsigned long) carry, (long) bn, (long) cn);
#endif
/* for rounding to nearest, we couldn't conclude up to here in the following
cases:
1. sh = 0, then cmp_low=0: we can either truncate, subtract one ulp
or add one ulp: -1 ulp < low(b)-low(c) < 1 ulp
2. sh > 0 but the low sh bits from high(b)-high(c) equal 2^(sh-1):
-0.5 ulp <= -1/2^sh < low(b)-low(c)-0.5 < 1/2^sh <= 0.5 ulp
we can't decide the rounding, in that case cmp_low=2:
either we truncate and flag=-1, or we add one ulp and flag=1
3. the low sh>0 bits from high(b)-high(c) equal 0: we know we have to
truncate but we can't decide the ternary value, here cmp_low=0:
-0.5 ulp <= -1/2^sh < low(b)-low(c) < 1/2^sh <= 0.5 ulp
we always truncate and inexact can be any of -1,0,1
*/
/* note: here cn might exceed cn0, in which case we consider a zero limb */
for (k = 0; (bn > 0) || (cn > 0); k = 1)
{
/* if cmp_low < 0, we know low(b) - low(c) < 0
if cmp_low > 0, we know low(b) - low(c) > 0
(more precisely if cmp_low = 2, low(b) - low(c) = 0.5 ulp so far)
if cmp_low = 0, so far low(b) - low(c) = 0 */
/* get next limbs */
bb = (bn > 0) ? bp[--bn] : 0;
if ((cn > 0) && (cn-- <= cn0))
cc = cp[cn];
else
cc = 0;
/* cmp_low compares low(b) and low(c) */
if (cmp_low == 0) /* case 1 or 3 */
cmp_low = (bb < cc) ? -2+k : (bb > cc) ? 1 : 0;
/* Case 1 for k=0 splits into 7 subcases:
1a: bb > cc + half
1b: bb = cc + half
1c: 0 < bb - cc < half
1d: bb = cc
1e: -half < bb - cc < 0
1f: bb - cc = -half
1g: bb - cc < -half
Case 2 splits into 3 subcases:
2a: bb > cc
2b: bb = cc
2c: bb < cc
Case 3 splits into 3 subcases:
3a: bb > cc
3b: bb = cc
3c: bb < cc
*/
/* the case rounding to nearest with sh=0 is special since one couldn't
subtract above 1/2 ulp in the trailing limb of the result */
if (rnd_mode == MPFR_RNDN && sh == 0 && k == 0) /* case 1 for k=0 */
{
mp_limb_t half = MPFR_LIMB_HIGHBIT;
/* add one ulp if bb > cc + half
truncate if cc - half < bb < cc + half
sub one ulp if bb < cc - half
*/
if (cmp_low < 0) /* bb < cc: -1 ulp < low(b) - low(c) < 0,
cases 1e, 1f and 1g */
{
if (cc >= half)
cc -= half;
else /* since bb < cc < half, bb+half < 2*half */
bb += half;
/* now we have bb < cc + half:
we have to subtract one ulp if bb < cc,
and truncate if bb > cc */
}
else if (cmp_low >= 0) /* bb >= cc, cases 1a to 1d */
{
if (cc < half)
cc += half;
else /* since bb >= cc >= half, bb - half >= 0 */
bb -= half;
/* now we have bb > cc - half: we have to add one ulp if bb > cc,
and truncate if bb < cc */
if (cmp_low > 0)
cmp_low = 2;
}
}
#ifdef DEBUG
printf ("k=%u bb=%lu cc=%lu cmp_low=%d\n", k,
(unsigned long) bb, (unsigned long) cc, cmp_low);
#endif
if (cmp_low < 0) /* low(b) - low(c) < 0: either truncate or subtract
one ulp */
{
if (rnd_mode == MPFR_RNDZ)
goto sub_one_ulp; /* set inexact=-1 */
else if (rnd_mode != MPFR_RNDN) /* round away */
{
inexact = 1;
goto truncate;
}
else /* round to nearest */
{
/* If cmp_low < 0 and bb > cc, then -0.5 ulp < low(b)-low(c) < 0,
whatever the value of sh.
If sh>0, then cmp_low < 0 implies that the initial neglected
sh bits were 0 (otherwise cmp_low=2 initially), thus the
weight of the new bits is less than 0.5 ulp too.
If k > 0 (and sh=0) this means that either the first neglected
limbs bb and cc were equal (thus cmp_low was 0 for k=0),
or we had bb - cc = -0.5 ulp or 0.5 ulp.
The last case is not possible here since we would have
cmp_low > 0 which is sticky.
In the first case (where we have cmp_low = -1), we truncate,
whereas in the 2nd case we have cmp_low = -2 and we subtract
one ulp.
*/
if (bb > cc || sh > 0 || cmp_low == -1)
{ /* -0.5 ulp < low(b)-low(c) < 0,
bb > cc corresponds to cases 1e and 1f1
sh > 0 corresponds to cases 3c and 3b3
cmp_low = -1 corresponds to case 1d3 (also 3b3) */
inexact = 1;
goto truncate;
}
else if (bb < cc) /* here sh = 0 and low(b)-low(c) < -0.5 ulp,
this corresponds to cases 1g and 1f3 */
goto sub_one_ulp;
/* the only case where we can't conclude is sh=0 and bb=cc,
i.e., we have low(b) - low(c) = -0.5 ulp (up to now), thus
we don't know if we must truncate or subtract one ulp.
Note: for sh=0 we can't have low(b) - low(c) = -0.5 ulp up to
now, since low(b) - low(c) > 1/2^sh */
}
}
else if (cmp_low > 0) /* 0 < low(b) - low(c): either truncate or
add one ulp */
{
if (rnd_mode == MPFR_RNDZ)
{
inexact = -1;
goto truncate;
}
else if (rnd_mode != MPFR_RNDN) /* round away */
goto add_one_ulp;
else /* round to nearest */
{
if (bb > cc)
{
/* if sh=0, then bb>cc means that low(b)-low(c) > 0.5 ulp,
and similarly when cmp_low=2 */
if (cmp_low == 2) /* cases 1a, 1b1, 2a and 2b1 */
goto add_one_ulp;
/* sh > 0 and cmp_low > 0: this implies that the sh initial
neglected bits were 0, and the remaining low(b)-low(c)>0,
but its weight is less than 0.5 ulp */
else /* 0 < low(b) - low(c) < 0.5 ulp, this corresponds to
cases 3a, 1d1 and 3b1 */
{
inexact = -1;
goto truncate;
}
}
else if (bb < cc) /* 0 < low(b) - low(c) < 0.5 ulp, cases 1c,
1b3, 2b3 and 2c */
{
inexact = -1;
goto truncate;
}
/* the only case where we can't conclude is bb=cc, i.e.,
low(b) - low(c) = 0.5 ulp (up to now), thus we don't know
if we must truncate or add one ulp. */
}
}
/* after k=0, we cannot conclude in the following cases, we split them
according to the values of bb and cc for k=1:
1b. sh=0 and cmp_low = 1 and bb-cc = half [around 0.5 ulp]
1b1. bb > cc: add one ulp, inex = 1
1b2: bb = cc: cannot conclude
1b3: bb < cc: truncate, inex = -1
1d. sh=0 and cmp_low = 0 and bb-cc = 0 [around 0]
1d1: bb > cc: truncate, inex = -1
1d2: bb = cc: cannot conclude
1d3: bb < cc: truncate, inex = +1
1f. sh=0 and cmp_low = -1 and bb-cc = -half [around -0.5 ulp]
1f1: bb > cc: truncate, inex = +1
1f2: bb = cc: cannot conclude
1f3: bb < cc: sub one ulp, inex = -1
2b. sh > 0 and cmp_low = 2 and bb=cc [around 0.5 ulp]
2b1. bb > cc: add one ulp, inex = 1
2b2: bb = cc: cannot conclude
2b3: bb < cc: truncate, inex = -1
3b. sh > 0 and cmp_low = 0 [around 0]
3b1. bb > cc: truncate, inex = -1
3b2: bb = cc: cannot conclude
3b3: bb < cc: truncate, inex = +1
*/
}
if ((rnd_mode == MPFR_RNDN) && cmp_low != 0)
{
/* even rounding rule */
if ((ap[0] >> sh) & 1)
{
if (cmp_low < 0)
goto sub_one_ulp;
else
goto add_one_ulp;
}
else
inexact = (cmp_low > 0) ? -1 : 1;
}
else
inexact = 0;
goto truncate;
sub_one_ulp: /* sub one unit in last place to a */
mpn_sub_1 (ap, ap, an, MPFR_LIMB_ONE << sh);
inexact = -1;
goto end_of_sub;
add_one_ulp: /* add one unit in last place to a */
if (MPFR_UNLIKELY(mpn_add_1 (ap, ap, an, MPFR_LIMB_ONE << sh)))
/* result is a power of 2: 11111111111111 + 1 = 1000000000000000 */
{
ap[an-1] = MPFR_LIMB_HIGHBIT;
add_exp = 1;
}
inexact = 1; /* result larger than exact value */
truncate:
if (MPFR_UNLIKELY((ap[an-1] >> (GMP_NUMB_BITS - 1)) == 0))
/* case 1 - epsilon */
{
ap[an-1] = MPFR_LIMB_HIGHBIT;
add_exp = 1;
}
end_of_sub:
/* we have to set MPFR_EXP(a) to MPFR_EXP(b) - cancel + add_exp, taking
care of underflows/overflows in that computation, and of the allowed
exponent range */
if (MPFR_LIKELY(cancel))
{
mpfr_exp_t exp_a;
cancel -= add_exp; /* OK: add_exp is an int equal to 0 or 1 */
exp_a = MPFR_GET_EXP (b) - cancel;
if (MPFR_UNLIKELY(exp_a < __gmpfr_emin))
{
MPFR_TMP_FREE(marker);
if (rnd_mode == MPFR_RNDN &&
(exp_a < __gmpfr_emin - 1 ||
(inexact >= 0 && mpfr_powerof2_raw (a))))
rnd_mode = MPFR_RNDZ;
return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
}
MPFR_SET_EXP (a, exp_a);
}
else /* cancel = 0: MPFR_EXP(a) <- MPFR_EXP(b) + add_exp */
{
/* in case cancel = 0, add_exp can still be 1, in case b is just
below a power of two, c is very small, prec(a) < prec(b),
and rnd=away or nearest */
mpfr_exp_t exp_b;
exp_b = MPFR_GET_EXP (b);
if (MPFR_UNLIKELY(add_exp && exp_b == __gmpfr_emax))
{
MPFR_TMP_FREE(marker);
return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
}
MPFR_SET_EXP (a, exp_b + add_exp);
}
MPFR_TMP_FREE(marker);
#ifdef DEBUG
printf ("result is a="); mpfr_print_binary(a); putchar('\n');
#endif
/* check that result is msb-normalized */
MPFR_ASSERTD(ap[an-1] > ~ap[an-1]);
MPFR_RET (inexact * MPFR_INT_SIGN (a));
}