xref: /netbsd-src/external/lgpl3/mpfr/dist/src/add1sp.c (revision ba125506a622fe649968631a56eba5d42ff57863)
1 /* mpfr_add1sp -- internal function to perform a "real" addition
2    All the op must have the same precision
3 
4 Copyright 2004-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6 
7 This file is part of the GNU MPFR Library.
8 
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.
13 
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.
18 
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 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
26 
27 #if MPFR_WANT_ASSERT >= 2
28 /* Check the result of mpfr_add1sp with mpfr_add1.
29 
30    Note: mpfr_add1sp itself has two algorithms: one always valid and one
31    faster for small precisions (up to 3 limbs). The latter one is disabled
32    if MPFR_GENERIC_ABI is defined. When MPFR_WANT_ASSERT >= 2, it could be
33    interesting to compare the results of these different algorithms. For
34    the time being, this is currently done by running the same code on the
35    same data with and without MPFR_GENERIC_ABI defined, where we have the
36    following comparisons in small precisions:
37      mpfr_add1sp slow <-> mpfr_add1 when MPFR_GENERIC_ABI is defined;
38      mpfr_add1sp fast <-> mpfr_add1 when MPFR_GENERIC_ABI is not defined.
39    By transitivity, the absence of failures implies that the 3 results are
40    the same.
41 */
42 
43 int mpfr_add1sp_ref (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)44 int mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
45 {
46   mpfr_t tmpa, tmpb, tmpc, tmpd;
47   mpfr_flags_t old_flags, flags, flags2;
48   int inexb, inexc, inexact, inexact2;
49 
50   if (rnd_mode == MPFR_RNDF)
51     return mpfr_add1sp_ref (a, b, c, rnd_mode);
52 
53   old_flags = __gmpfr_flags;
54 
55   mpfr_init2 (tmpa, MPFR_PREC (a));
56   mpfr_init2 (tmpb, MPFR_PREC (b));
57   mpfr_init2 (tmpc, MPFR_PREC (c));
58 
59   inexb = mpfr_set (tmpb, b, MPFR_RNDN);
60   MPFR_ASSERTN (inexb == 0);
61 
62   inexc = mpfr_set (tmpc, c, MPFR_RNDN);
63   MPFR_ASSERTN (inexc == 0);
64 
65   MPFR_ASSERTN (__gmpfr_flags == old_flags);
66 
67   if (MPFR_GET_EXP (tmpb) < MPFR_GET_EXP (tmpc))
68     {
69       /* The sign for the result will be taken from the second argument
70          (= first input value, which is tmpb). */
71       MPFR_ALIAS (tmpd, tmpc, MPFR_SIGN (tmpb), MPFR_EXP (tmpc));
72       inexact2 = mpfr_add1 (tmpa, tmpd, tmpb, rnd_mode);
73     }
74   else
75     {
76       inexact2 = mpfr_add1 (tmpa, tmpb, tmpc, rnd_mode);
77     }
78   flags2 = __gmpfr_flags;
79 
80   __gmpfr_flags = old_flags;
81   inexact = mpfr_add1sp_ref (a, b, c, rnd_mode);
82   flags = __gmpfr_flags;
83 
84   /* Convert the ternary values to (-1,0,1). */
85   inexact2 = VSIGN (inexact2);
86   inexact = VSIGN (inexact);
87 
88   if (! mpfr_equal_p (tmpa, a) || inexact != inexact2 || flags != flags2)
89     {
90       fprintf (stderr, "add1 & add1sp return different values for %s\n"
91                "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
92                mpfr_print_rnd_mode (rnd_mode),
93                (unsigned long) MPFR_PREC (a),
94                (unsigned long) MPFR_PREC (b),
95                (unsigned long) MPFR_PREC (c));
96       mpfr_fdump (stderr, tmpb);
97       fprintf (stderr, "C = ");
98       mpfr_fdump (stderr, tmpc);
99       fprintf (stderr, "\nadd1  : ");
100       mpfr_fdump (stderr, tmpa);
101       fprintf (stderr, "add1sp: ");
102       mpfr_fdump (stderr, a);
103       fprintf (stderr, "add1  : ternary = %2d, flags =", inexact2);
104       flags_fout (stderr, flags2);
105       fprintf (stderr, "add1sp: ternary = %2d, flags =", inexact);
106       flags_fout (stderr, flags);
107       MPFR_ASSERTN (0);
108     }
109   mpfr_clears (tmpa, tmpb, tmpc, (mpfr_ptr) 0);
110   return inexact;
111 }
112 # define mpfr_add1sp mpfr_add1sp_ref
113 #endif  /* MPFR_WANT_ASSERT >= 2 */
114 
115 #if !defined(MPFR_GENERIC_ABI)
116 
117 #if defined(MPFR_WANT_PROVEN_CODE) && GMP_NUMB_BITS == 64 && \
118   UINT_MAX == 0xffffffff && MPFR_PREC_BITS == 64 && \
119   _MPFR_PREC_FORMAT == 3 && _MPFR_EXP_FORMAT == _MPFR_PREC_FORMAT
120 
121 /* The code assumes that mp_limb_t has 64 bits exactly, unsigned int
122    has 32 bits exactly, mpfr_prec_t and mpfr_exp_t are of type long,
123    which has 64 bits exactly. */
124 
125 #include "add1sp1_extracted.c"
126 
127 #else
128 
129 /* same as mpfr_add1sp, but for p < GMP_NUMB_BITS */
130 static int
mpfr_add1sp1(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode,mpfr_prec_t p)131 mpfr_add1sp1 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
132               mpfr_prec_t p)
133 {
134   mpfr_exp_t bx = MPFR_GET_EXP (b);
135   mpfr_exp_t cx = MPFR_GET_EXP (c);
136   mp_limb_t *ap = MPFR_MANT(a);
137   mp_limb_t *bp = MPFR_MANT(b);
138   mp_limb_t *cp = MPFR_MANT(c);
139   mpfr_prec_t sh = GMP_NUMB_BITS - p;
140   mp_limb_t rb; /* round bit */
141   mp_limb_t sb; /* sticky bit */
142   mp_limb_t a0; /* to store the result */
143   mp_limb_t mask;
144   mpfr_uexp_t d;
145 
146   MPFR_ASSERTD(p < GMP_NUMB_BITS);
147 
148   if (bx == cx)
149     {
150       /* The following line is probably better than
151            a0 = MPFR_LIMB_HIGHBIT | ((bp[0] + cp[0]) >> 1);
152          as it has less dependency and doesn't need a long constant on some
153          processors. On ARM, it can also probably benefit from shift-and-op
154          in a better way. Timings cannot be conclusive. */
155       a0 = (bp[0] >> 1) + (cp[0] >> 1);
156       bx ++;
157       rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
158       ap[0] = a0 ^ rb;
159       sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */
160     }
161   else
162     {
163       if (bx < cx)  /* swap b and c */
164         {
165           mpfr_exp_t tx;
166           mp_limb_t *tp;
167           tx = bx; bx = cx; cx = tx;
168           tp = bp; bp = cp; cp = tp;
169         }
170       MPFR_ASSERTD (bx > cx);
171       d = (mpfr_uexp_t) bx - cx;
172       mask = MPFR_LIMB_MASK(sh);
173       /* TODO: Should the case d < sh be removed, i.e. seen as a particular
174          case of d < GMP_NUMB_BITS? This case would do a bit more operations
175          but a test would be removed, avoiding pipeline stall issues. */
176       if (d < sh)
177         {
178           /* we can shift c by d bits to the right without losing any bit,
179              moreover we can shift one more if there is an exponent increase */
180           a0 = bp[0] + (cp[0] >> d);
181           if (a0 < bp[0]) /* carry */
182             {
183               MPFR_ASSERTD ((a0 & MPFR_LIMB_ONE) == 0);
184               a0 = MPFR_LIMB_HIGHBIT | (a0 >> 1);
185               bx ++;
186             }
187           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
188           sb = (a0 & mask) ^ rb;
189           ap[0] = a0 & ~mask;
190         }
191       else if (d < GMP_NUMB_BITS) /* sh <= d < GMP_NUMB_BITS */
192         {
193           sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */
194           a0 = bp[0] + (cp[0] >> d);
195           if (a0 < bp[0]) /* carry */
196             {
197               sb |= a0 & MPFR_LIMB_ONE;
198               a0 = MPFR_LIMB_HIGHBIT | (a0 >> 1);
199               bx ++;
200             }
201           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
202           sb |= (a0 & mask) ^ rb;
203           ap[0] = a0 & ~mask;
204         }
205       else /* d >= GMP_NUMB_BITS */
206         {
207           ap[0] = bp[0];
208           rb = 0; /* since p < GMP_NUMB_BITS */
209           sb = 1; /* since c <> 0 */
210         }
211     }
212 
213   /* Note: we could keep the output significand in a0 for the rounding,
214      and only store it in ap[0] at the very end, but this seems slower
215      on average (but better for the worst case). */
216 
217   /* now perform rounding */
218   if (MPFR_UNLIKELY(bx > __gmpfr_emax))
219     return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
220 
221   MPFR_SET_EXP (a, bx);
222   if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
223     MPFR_RET(0);
224   else if (rnd_mode == MPFR_RNDN)
225     {
226       /* the condition below should be rb == 0 || (rb != 0 && ...), but this
227          is equivalent to rb == 0 || (...) */
228       if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
229         goto truncate;
230       else
231         goto add_one_ulp;
232     }
233   else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
234     {
235     truncate:
236       MPFR_RET(-MPFR_SIGN(a));
237     }
238   else /* round away from zero */
239     {
240     add_one_ulp:
241       ap[0] += MPFR_LIMB_ONE << sh;
242       if (MPFR_UNLIKELY(ap[0] == 0))
243         {
244           ap[0] = MPFR_LIMB_HIGHBIT;
245           /* no need to have MPFR_LIKELY here, since we are in a rare branch */
246           if (bx + 1 <= __gmpfr_emax)
247             MPFR_SET_EXP (a, bx + 1);
248           else /* overflow */
249             return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
250         }
251       MPFR_RET(MPFR_SIGN(a));
252     }
253 }
254 
255 #endif /* MPFR_WANT_PROVEN_CODE */
256 
257 /* same as mpfr_add1sp, but for p = GMP_NUMB_BITS */
258 static int
mpfr_add1sp1n(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)259 mpfr_add1sp1n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
260 {
261   mpfr_exp_t bx = MPFR_GET_EXP (b);
262   mpfr_exp_t cx = MPFR_GET_EXP (c);
263   mp_limb_t *ap = MPFR_MANT(a);
264   mp_limb_t *bp = MPFR_MANT(b);
265   mp_limb_t *cp = MPFR_MANT(c);
266   mp_limb_t rb; /* round bit */
267   mp_limb_t sb; /* sticky bit */
268   mp_limb_t a0; /* to store the result */
269   mpfr_uexp_t d;
270 
271   MPFR_ASSERTD(MPFR_PREC (a) == GMP_NUMB_BITS);
272   MPFR_ASSERTD(MPFR_PREC (b) == GMP_NUMB_BITS);
273   MPFR_ASSERTD(MPFR_PREC (c) == GMP_NUMB_BITS);
274 
275   if (bx == cx)
276     {
277       a0 = bp[0] + cp[0];
278       rb = a0 & MPFR_LIMB_ONE;
279       ap[0] = MPFR_LIMB_HIGHBIT | (a0 >> 1);
280       bx ++;
281       sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */
282     }
283   else
284     {
285       if (bx < cx)  /* swap b and c */
286         {
287           mpfr_exp_t tx;
288           mp_limb_t *tp;
289           tx = bx; bx = cx; cx = tx;
290           tp = bp; bp = cp; cp = tp;
291         }
292       MPFR_ASSERTD (bx > cx);
293       d = (mpfr_uexp_t) bx - cx;
294       if (d < GMP_NUMB_BITS) /* 1 <= d < GMP_NUMB_BITS */
295         {
296           a0 = bp[0] + (cp[0] >> d);
297           sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */
298           if (a0 < bp[0]) /* carry */
299             {
300               ap[0] = MPFR_LIMB_HIGHBIT | (a0 >> 1);
301               rb = a0 & MPFR_LIMB_ONE;
302               bx ++;
303             }
304           else /* no carry */
305             {
306               ap[0] = a0;
307               rb = sb & MPFR_LIMB_HIGHBIT;
308               sb &= ~MPFR_LIMB_HIGHBIT;
309             }
310         }
311       else /* d >= GMP_NUMB_BITS */
312         {
313           sb = d != GMP_NUMB_BITS || cp[0] != MPFR_LIMB_HIGHBIT;
314           ap[0] = bp[0];
315           rb = d == GMP_NUMB_BITS;
316         }
317     }
318 
319   /* Note: we could keep the output significand in a0 for the rounding,
320      and only store it in ap[0] at the very end, but this seems slower
321      on average (but better for the worst case). */
322 
323   /* now perform rounding */
324   if (MPFR_UNLIKELY(bx > __gmpfr_emax))
325     return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
326 
327   MPFR_SET_EXP (a, bx);
328   if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
329     MPFR_RET(0);
330   else if (rnd_mode == MPFR_RNDN)
331     {
332       /* the condition below should be rb == 0 || (rb != 0 && ...), but this
333          is equivalent to rb == 0 || (...) */
334       if (rb == 0 || (sb == 0 && (ap[0] & MPFR_LIMB_ONE) == 0))
335         goto truncate;
336       else
337         goto add_one_ulp;
338     }
339   else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
340     {
341     truncate:
342       MPFR_RET(-MPFR_SIGN(a));
343     }
344   else /* round away from zero */
345     {
346     add_one_ulp:
347       ap[0] += MPFR_LIMB_ONE;
348       if (MPFR_UNLIKELY(ap[0] == 0))
349         {
350           ap[0] = MPFR_LIMB_HIGHBIT;
351           /* no need to have MPFR_LIKELY here, since we are in a rare branch */
352           if (bx + 1 <= __gmpfr_emax)
353             MPFR_SET_EXP (a, bx + 1);
354           else /* overflow */
355             return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
356         }
357       MPFR_RET(MPFR_SIGN(a));
358     }
359 }
360 
361 /* same as mpfr_add1sp, but for GMP_NUMB_BITS < p < 2*GMP_NUMB_BITS */
362 static int
mpfr_add1sp2(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode,mpfr_prec_t p)363 mpfr_add1sp2 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
364               mpfr_prec_t p)
365 {
366   mpfr_exp_t bx = MPFR_GET_EXP (b);
367   mpfr_exp_t cx = MPFR_GET_EXP (c);
368   mp_limb_t *ap = MPFR_MANT(a);
369   mp_limb_t *bp = MPFR_MANT(b);
370   mp_limb_t *cp = MPFR_MANT(c);
371   mpfr_prec_t sh = 2*GMP_NUMB_BITS - p;
372   mp_limb_t rb; /* round bit */
373   mp_limb_t sb; /* sticky bit */
374   mp_limb_t a1, a0;
375   mp_limb_t mask;
376   mpfr_uexp_t d;
377 
378   MPFR_ASSERTD(GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS);
379 
380   if (bx == cx)
381     {
382       /* since bp[1], cp[1] >= MPFR_LIMB_HIGHBIT, a carry always occurs */
383       a0 = bp[0] + cp[0];
384       a1 = bp[1] + cp[1] + (a0 < bp[0]);
385       a0 = (a0 >> 1) | (a1 << (GMP_NUMB_BITS - 1));
386       bx ++;
387       rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
388       ap[1] = MPFR_LIMB_HIGHBIT | (a1 >> 1);
389       ap[0] = a0 ^ rb;
390       sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */
391     }
392   else
393     {
394       if (bx < cx)  /* swap b and c */
395         {
396           mpfr_exp_t tx;
397           mp_limb_t *tp;
398           tx = bx; bx = cx; cx = tx;
399           tp = bp; bp = cp; cp = tp;
400         }
401       MPFR_ASSERTD (bx > cx);
402       d = (mpfr_uexp_t) bx - cx;
403       mask = MPFR_LIMB_MASK(sh);
404       if (d < GMP_NUMB_BITS) /* 0 < d < GMP_NUMB_BITS */
405         {
406           sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */
407           a0 = bp[0] + ((cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d));
408           a1 = bp[1] + (cp[1] >> d) + (a0 < bp[0]);
409           if (a1 < bp[1]) /* carry in high word */
410             {
411             exponent_shift:
412               sb |= a0 & MPFR_LIMB_ONE;
413               /* shift a by 1 */
414               a0 = (a1 << (GMP_NUMB_BITS - 1)) | (a0 >> 1);
415               ap[1] = MPFR_LIMB_HIGHBIT | (a1 >> 1);
416               bx ++;
417             }
418           else
419             ap[1] = a1;
420           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
421           sb |= (a0 & mask) ^ rb;
422           ap[0] = a0 & ~mask;
423         }
424       else if (d < 2*GMP_NUMB_BITS) /* GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS */
425         {
426           sb = (d == GMP_NUMB_BITS) ? cp[0]
427             : cp[0] | (cp[1] << (2*GMP_NUMB_BITS-d));
428           a0 = bp[0] + (cp[1] >> (d - GMP_NUMB_BITS));
429           a1 = bp[1] + (a0 < bp[0]);
430           if (a1 == 0)
431             goto exponent_shift;
432           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
433           sb |= (a0 & mask) ^ rb;
434           ap[0] = a0 & ~mask;
435           ap[1] = a1;
436         }
437       else /* d >= 2*GMP_NUMB_BITS */
438         {
439           ap[0] = bp[0];
440           ap[1] = bp[1];
441           rb = 0; /* since p < 2*GMP_NUMB_BITS */
442           sb = 1; /* since c <> 0 */
443         }
444     }
445 
446   /* now perform rounding */
447   if (MPFR_UNLIKELY(bx > __gmpfr_emax))
448     return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
449 
450   MPFR_SET_EXP (a, bx);
451   if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
452     MPFR_RET(0);
453   else if (rnd_mode == MPFR_RNDN)
454     {
455       if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
456         goto truncate;
457       else
458         goto add_one_ulp;
459     }
460   else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
461     {
462     truncate:
463       MPFR_RET(-MPFR_SIGN(a));
464     }
465   else /* round away from zero */
466     {
467     add_one_ulp:
468       ap[0] += MPFR_LIMB_ONE << sh;
469       ap[1] += (ap[0] == 0);
470       if (MPFR_UNLIKELY(ap[1] == 0))
471         {
472           ap[1] = MPFR_LIMB_HIGHBIT;
473           /* no need to have MPFR_LIKELY here, since we are in a rare branch */
474           if (bx + 1 <= __gmpfr_emax)
475             MPFR_SET_EXP (a, bx + 1);
476           else /* overflow */
477             return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
478         }
479       MPFR_RET(MPFR_SIGN(a));
480     }
481 }
482 
483 /* same as mpfr_add1sp, but for p = 2*GMP_NUMB_BITS */
484 static int
mpfr_add1sp2n(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)485 mpfr_add1sp2n (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
486 {
487   mpfr_exp_t bx = MPFR_GET_EXP (b);
488   mpfr_exp_t cx = MPFR_GET_EXP (c);
489   mp_limb_t *ap = MPFR_MANT(a);
490   mp_limb_t *bp = MPFR_MANT(b);
491   mp_limb_t *cp = MPFR_MANT(c);
492   mp_limb_t rb; /* round bit */
493   mp_limb_t sb; /* sticky bit */
494   mp_limb_t a1, a0;
495   mpfr_uexp_t d;
496 
497   if (bx == cx)
498     {
499       /* since bp[1], cp[1] >= MPFR_LIMB_HIGHBIT, a carry always occurs */
500       a0 = bp[0] + cp[0];
501       a1 = bp[1] + cp[1] + (a0 < bp[0]);
502       rb = a0 & MPFR_LIMB_ONE;
503       sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */
504       ap[0] = (a1 << (GMP_NUMB_BITS - 1)) | (a0 >> 1);
505       ap[1] = MPFR_LIMB_HIGHBIT | (a1 >> 1);
506       bx ++;
507     }
508   else
509     {
510       if (bx < cx)  /* swap b and c */
511         {
512           mpfr_exp_t tx;
513           mp_limb_t *tp;
514           tx = bx; bx = cx; cx = tx;
515           tp = bp; bp = cp; cp = tp;
516         }
517       MPFR_ASSERTD (bx > cx);
518       d = (mpfr_uexp_t) bx - cx;
519       if (d >= 2 * GMP_NUMB_BITS)
520         {
521           if (d == 2 * GMP_NUMB_BITS)
522             {
523               rb = 1;
524               sb = (cp[0] != MPFR_LIMB_ZERO ||
525                     cp[1] > MPFR_LIMB_HIGHBIT);
526             }
527           else
528             {
529               rb = 0;
530               sb = 1;
531             }
532           ap[0] = bp[0];
533           ap[1] = bp[1];
534         }
535       else
536         {
537           /* First, compute (a0,a1) = b + (c >> d), and determine sb from
538              the bits shifted out such that (MSB, other bits) is regarded
539              as (rounding bit, sticky bit), assuming no carry. */
540           if (d < GMP_NUMB_BITS) /* 0 < d < GMP_NUMB_BITS */
541             {
542               sb = cp[0] << (GMP_NUMB_BITS - d);
543               a0 = bp[0] + ((cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d));
544               a1 = bp[1] + (cp[1] >> d) + (a0 < bp[0]);
545             }
546           else /* GMP_NUMB_BITS <= d < 2 * GMP_NUMB_BITS */
547             {
548               /* The most significant bit of sb should be the rounding bit,
549                  while the other bits represent the sticky bit:
550                  * if d = GMP_NUMB_BITS, we get cp[0];
551                  * if d > GMP_NUMB_BITS: we get the least d-GMP_NUMB_BITS bits
552                    of cp[1], and those from cp[0] as the LSB of sb. */
553               sb = (d == GMP_NUMB_BITS) ? cp[0]
554                 : (cp[1] << (2*GMP_NUMB_BITS-d)) | (cp[0] != 0);
555               a0 = bp[0] + (cp[1] >> (d - GMP_NUMB_BITS));
556               a1 = bp[1] + (a0 < bp[0]);
557             }
558           if (a1 < bp[1]) /* carry in high word */
559             {
560               rb = a0 << (GMP_NUMB_BITS - 1);
561               /* and sb is the real sticky bit. */
562               /* Shift the result by 1 to the right. */
563               ap[0] = (a1 << (GMP_NUMB_BITS - 1)) | (a0 >> 1);
564               ap[1] = MPFR_LIMB_HIGHBIT | (a1 >> 1);
565               bx ++;
566             }
567           else
568             {
569               rb = MPFR_LIMB_MSB (sb);
570               sb <<= 1;
571               ap[0] = a0;
572               ap[1] = a1;
573             }
574         }
575     }
576 
577   /* now perform rounding */
578   if (MPFR_UNLIKELY(bx > __gmpfr_emax))
579     return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
580 
581   MPFR_SET_EXP (a, bx);
582   if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
583     MPFR_RET(0);
584   else if (rnd_mode == MPFR_RNDN)
585     {
586       if (rb == 0 || (sb == 0 && (ap[0] & MPFR_LIMB_ONE) == 0))
587         goto truncate;
588       else
589         goto add_one_ulp;
590     }
591   else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
592     {
593     truncate:
594       MPFR_RET(-MPFR_SIGN(a));
595     }
596   else /* round away from zero */
597     {
598     add_one_ulp:
599       ap[0] += MPFR_LIMB_ONE;
600       ap[1] += (ap[0] == 0);
601       if (MPFR_UNLIKELY(ap[1] == 0))
602         {
603           ap[1] = MPFR_LIMB_HIGHBIT;
604           /* no need to have MPFR_LIKELY here, since we are in a rare branch */
605           if (bx + 1 <= __gmpfr_emax)
606             MPFR_SET_EXP (a, bx + 1);
607           else /* overflow */
608             return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
609         }
610       MPFR_RET(MPFR_SIGN(a));
611     }
612 }
613 
614 /* same as mpfr_add1sp, but for 2*GMP_NUMB_BITS < p < 3*GMP_NUMB_BITS */
615 static int
mpfr_add1sp3(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode,mpfr_prec_t p)616 mpfr_add1sp3 (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode,
617               mpfr_prec_t p)
618 {
619   mpfr_exp_t bx = MPFR_GET_EXP (b);
620   mpfr_exp_t cx = MPFR_GET_EXP (c);
621   mp_limb_t *ap = MPFR_MANT(a);
622   mp_limb_t *bp = MPFR_MANT(b);
623   mp_limb_t *cp = MPFR_MANT(c);
624   mpfr_prec_t sh = 3*GMP_NUMB_BITS - p;
625   mp_limb_t rb; /* round bit */
626   mp_limb_t sb; /* sticky bit */
627   mp_limb_t a2, a1, a0;
628   mp_limb_t mask;
629   mpfr_uexp_t d;
630 
631   MPFR_ASSERTD(2 * GMP_NUMB_BITS < p && p < 3 * GMP_NUMB_BITS);
632 
633   if (bx == cx)
634     {
635       /* since bp[2], cp[2] >= MPFR_LIMB_HIGHBIT, a carry always occurs */
636       a0 = bp[0] + cp[0];
637       a1 = bp[1] + cp[1] + (a0 < bp[0]);
638       a2 = bp[2] + cp[2] + (a1 < bp[1] || (a1 == bp[1] && a0 < bp[0]));
639       /* since p < 3 * GMP_NUMB_BITS, we lose no bit in a0 >> 1 */
640       a0 = (a1 << (GMP_NUMB_BITS - 1)) | (a0 >> 1);
641       bx ++;
642       rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
643       ap[2] = MPFR_LIMB_HIGHBIT | (a2 >> 1);
644       ap[1] = (a2 << (GMP_NUMB_BITS - 1)) | (a1 >> 1);
645       ap[0] = a0 ^ rb;
646       sb = 0; /* since b + c fits on p+1 bits, the sticky bit is zero */
647     }
648   else
649     {
650       if (bx < cx)  /* swap b and c */
651         {
652           mpfr_exp_t tx;
653           mp_limb_t *tp;
654           tx = bx; bx = cx; cx = tx;
655           tp = bp; bp = cp; cp = tp;
656         }
657       MPFR_ASSERTD (bx > cx);
658       d = (mpfr_uexp_t) bx - cx;
659       mask = MPFR_LIMB_MASK(sh);
660       if (d < GMP_NUMB_BITS) /* 0 < d < GMP_NUMB_BITS */
661         {
662           mp_limb_t cy;
663           sb = cp[0] << (GMP_NUMB_BITS - d); /* bits from cp[-1] after shift */
664           a0 = bp[0] + ((cp[1] << (GMP_NUMB_BITS - d)) | (cp[0] >> d));
665           a1 = bp[1] + ((cp[2] << (GMP_NUMB_BITS - d)) | (cp[1] >> d))
666             + (a0 < bp[0]);
667           cy = a1 < bp[1] || (a1 == bp[1] && a0 < bp[0]); /* carry in a1 */
668           a2 = bp[2] + (cp[2] >> d) + cy;
669           if (a2 < bp[2] || (a2 == bp[2] && cy)) /* carry in high word */
670             {
671             exponent_shift:
672               sb |= a0 & MPFR_LIMB_ONE;
673               /* shift a by 1 */
674               a0 = (a1 << (GMP_NUMB_BITS - 1)) | (a0 >> 1);
675               ap[1] = (a2 << (GMP_NUMB_BITS - 1)) | (a1 >> 1);
676               ap[2] = MPFR_LIMB_HIGHBIT | (a2 >> 1);
677               bx ++;
678             }
679           else
680             {
681               ap[1] = a1;
682               ap[2] = a2;
683             }
684           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
685           sb |= (a0 & mask) ^ rb;
686           ap[0] = a0 & ~mask;
687         }
688       else if (d < 2*GMP_NUMB_BITS) /* GMP_NUMB_BITS <= d < 2*GMP_NUMB_BITS */
689         {
690           mp_limb_t c0shifted;
691           sb = (d == GMP_NUMB_BITS) ? cp[0]
692             : (cp[1] << (2*GMP_NUMB_BITS - d)) | cp[0];
693           c0shifted = (d == GMP_NUMB_BITS) ? cp[1]
694             : (cp[2] << (2*GMP_NUMB_BITS-d)) | (cp[1] >> (d - GMP_NUMB_BITS));
695           a0 = bp[0] + c0shifted;
696           a1 = bp[1] + (cp[2] >> (d - GMP_NUMB_BITS)) + (a0 < bp[0]);
697           /* if a1 < bp[1], there was a carry in the above addition,
698              or when a1 = bp[1] and one of the added terms is non-zero
699              (the sum of cp[2] >> (d - GMP_NUMB_BITS) and a0 < bp[0]
700              is at most 2^GMP_NUMB_BITS-d) */
701           a2 = bp[2] + ((a1 < bp[1]) || (a1 == bp[1] && a0 < bp[0]));
702           if (a2 == 0)
703             goto exponent_shift;
704           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
705           sb |= (a0 & mask) ^ rb;
706           ap[0] = a0 & ~mask;
707           ap[1] = a1;
708           ap[2] = a2;
709         }
710       else if (d < 3*GMP_NUMB_BITS) /* 2*GMP_NUMB_BITS<=d<3*GMP_NUMB_BITS */
711         {
712           MPFR_ASSERTD (2*GMP_NUMB_BITS <= d && d < 3*GMP_NUMB_BITS);
713           sb = (d == 2*GMP_NUMB_BITS ? 0 : cp[2] << (3*GMP_NUMB_BITS - d))
714             | cp[1] | cp[0];
715           a0 = bp[0] + (cp[2] >> (d - 2*GMP_NUMB_BITS));
716           a1 = bp[1] + (a0 < bp[0]);
717           a2 = bp[2] + (a1 < bp[1]);
718           if (a2 == 0)
719             goto exponent_shift;
720           rb = a0 & (MPFR_LIMB_ONE << (sh - 1));
721           sb |= (a0 & mask) ^ rb;
722           ap[0] = a0 & ~mask;
723           ap[1] = a1;
724           ap[2] = a2;
725         }
726       else /* d >= 2*GMP_NUMB_BITS */
727         {
728           ap[0] = bp[0];
729           ap[1] = bp[1];
730           ap[2] = bp[2];
731           rb = 0; /* since p < 3*GMP_NUMB_BITS */
732           sb = 1; /* since c <> 0 */
733         }
734     }
735 
736   /* now perform rounding */
737   if (MPFR_UNLIKELY(bx > __gmpfr_emax))
738     return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
739 
740   MPFR_SET_EXP (a, bx);
741   if ((rb == 0 && sb == 0) || rnd_mode == MPFR_RNDF)
742     MPFR_RET(0);
743   else if (rnd_mode == MPFR_RNDN)
744     {
745       if (rb == 0 || (sb == 0 && (ap[0] & (MPFR_LIMB_ONE << sh)) == 0))
746         goto truncate;
747       else
748         goto add_one_ulp;
749     }
750   else if (MPFR_IS_LIKE_RNDZ(rnd_mode, MPFR_IS_NEG(a)))
751     {
752     truncate:
753       MPFR_RET(-MPFR_SIGN(a));
754     }
755   else /* round away from zero */
756     {
757     add_one_ulp:
758       ap[0] += MPFR_LIMB_ONE << sh;
759       ap[1] += (ap[0] == 0);
760       ap[2] += (ap[1] == 0 && ap[0] == 0);
761       if (MPFR_UNLIKELY(ap[2] == 0))
762         {
763           ap[2] = MPFR_LIMB_HIGHBIT;
764           /* no need to have MPFR_LIKELY here, since we are in a rare branch */
765           if (bx + 1 <= __gmpfr_emax)
766             MPFR_SET_EXP (a, bx + 1);
767           else /* overflow */
768             return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
769         }
770       MPFR_RET(MPFR_SIGN(a));
771     }
772 }
773 
774 #endif /* !defined(MPFR_GENERIC_ABI) */
775 
776 /* {ap, n} <- {bp, n} + {cp + q, n - q} >> r where d = q * GMP_NUMB_BITS + r.
777    Return the carry at ap[n+1] (0 or 1) and set *low so that:
778    * the most significant bit of *low would be that of ap[-1] if we would
779      compute one more limb of the (infinite) addition
780    * the GMP_NUMB_BITS-1 least significant bits of *low are zero iff all bits
781      of ap[-1], ap[-2], ... would be zero (except the most significant bit
782      of ap[-1).
783    Assume 0 < d < GMP_NUMB_BITS*n. */
784 static mp_limb_t
mpfr_addrsh(mp_limb_t * ap,mp_limb_t * bp,mp_limb_t * cp,mp_size_t n,mp_size_t d,mp_limb_t * low)785 mpfr_addrsh (mp_limb_t *ap, mp_limb_t *bp, mp_limb_t *cp, mp_size_t n,
786              mp_size_t d, mp_limb_t *low)
787 {
788   mp_limb_t cy, cy2, c_shifted;
789   mp_size_t i;
790 
791   if (d < GMP_NUMB_BITS)
792     {
793       /* {ap, n} <- {bp, n} + {cp, n} >> d */
794       MPFR_ASSERTD (d > 0);
795       /* thus 0 < GMP_NUMB_BITS - d < GMP_NUMB_BITS */
796       *low = cp[0] << (GMP_NUMB_BITS - d);
797       for (i = 0, cy = 0; i < n - 1; i++)
798         {
799           c_shifted = (cp[i+1] << (GMP_NUMB_BITS - d)) | (cp[i] >> d);
800           ap[i] = bp[i] + c_shifted;
801           cy2 = ap[i] < c_shifted;
802           ap[i] += cy;
803           cy = cy2 + (ap[i] < cy);
804         }
805       /* most significant limb is special */
806       c_shifted = cp[i] >> d;
807       ap[i] = bp[i] + c_shifted;
808       cy2 = ap[i] < c_shifted;
809       ap[i] += cy;
810       cy = cy2 + (ap[i] < cy);
811     }
812   else /* d >= GMP_NUMB_BITS */
813     {
814       mp_size_t q = d / GMP_NUMB_BITS;
815       mpfr_uexp_t r = d % GMP_NUMB_BITS;
816       if (r == 0)
817         {
818           MPFR_ASSERTD(q > 0);
819           *low = cp[q-1];
820           for (i = 0; i < q-1; i++)
821             *low |= !!cp[i];
822           cy = mpn_add_n (ap, bp, cp + q, n - q);
823           cy = mpn_add_1 (ap + n - q, bp + n - q, q, cy);
824         }
825       else /* 0 < r < GMP_NUMB_BITS */
826         {
827           *low = cp[q] << (GMP_NUMB_BITS - r);
828           for (i = 0; i < q; i++)
829             *low |= !!cp[i];
830           for (i = 0, cy = 0; i < n - q - 1; i++)
831             {
832               c_shifted = (cp[q+i+1] << (GMP_NUMB_BITS - r)) | (cp[q+i] >> r);
833               ap[i] = bp[i] + c_shifted;
834               cy2 = ap[i] < c_shifted;
835               ap[i] += cy;
836               cy = cy2 + (ap[i] < cy);
837             }
838           /* most significant limb of c is special */
839           MPFR_ASSERTD(i == n - q - 1);
840           c_shifted = cp[n-1] >> r;
841           ap[i] = bp[i] + c_shifted;
842           cy2 = ap[i] < c_shifted;
843           ap[i] += cy;
844           cy = cy2 + (ap[i] < cy);
845           /* upper limbs are copied */
846           cy = mpn_add_1 (ap + n - q, bp + n - q, q, cy);
847         }
848     }
849   return cy;
850 }
851 
852 /* compute sign(b) * (|b| + |c|).
853    Returns 0 iff result is exact,
854    a negative value when the result is less than the exact value,
855    a positive value otherwise. */
856 MPFR_HOT_FUNCTION_ATTR int
mpfr_add1sp(mpfr_ptr a,mpfr_srcptr b,mpfr_srcptr c,mpfr_rnd_t rnd_mode)857 mpfr_add1sp (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
858 {
859   mpfr_uexp_t d;
860   mpfr_prec_t p;
861   unsigned int sh;
862   mp_size_t n;
863   mp_limb_t *ap = MPFR_MANT(a);
864   mpfr_exp_t bx;
865   mp_limb_t limb, rb, sb;
866   int inexact;
867   int neg;
868 
869   MPFR_ASSERTD(MPFR_PREC(a) == MPFR_PREC(b) && MPFR_PREC(b) == MPFR_PREC(c));
870   MPFR_ASSERTD(MPFR_IS_PURE_FP(b));
871   MPFR_ASSERTD(MPFR_IS_PURE_FP(c));
872 
873   MPFR_SET_SAME_SIGN (a, b);
874 
875   /* Read prec and num of limbs */
876   p = MPFR_GET_PREC (b);
877 
878 #if !defined(MPFR_GENERIC_ABI)
879   if (p < GMP_NUMB_BITS)
880     return mpfr_add1sp1 (a, b, c, rnd_mode, p);
881 
882   if (GMP_NUMB_BITS < p && p < 2 * GMP_NUMB_BITS)
883     return mpfr_add1sp2 (a, b, c, rnd_mode, p);
884 
885   /* we put this test after the mpfr_add1sp2() call, to avoid slowing down
886      the more frequent case GMP_NUMB_BITS < p < 2 * GMP_NUMB_BITS */
887   if (p == GMP_NUMB_BITS)
888     return mpfr_add1sp1n (a, b, c, rnd_mode);
889 
890   if (2 * GMP_NUMB_BITS < p && p < 3 * GMP_NUMB_BITS)
891     return mpfr_add1sp3 (a, b, c, rnd_mode, p);
892 
893   if (p == 2 * GMP_NUMB_BITS)
894     return mpfr_add1sp2n (a, b, c, rnd_mode);
895 #endif
896 
897   /* We need to get the sign before the possible exchange. */
898   neg = MPFR_IS_NEG (b);
899 
900   bx = MPFR_GET_EXP(b);
901   if (bx < MPFR_GET_EXP(c))
902     {
903       mpfr_srcptr t = b;
904       bx = MPFR_GET_EXP(c);
905       b = c;
906       c = t;
907     }
908   MPFR_ASSERTD(bx >= MPFR_GET_EXP(c));
909 
910   n = MPFR_PREC2LIMBS (p);
911   MPFR_UNSIGNED_MINUS_MODULO(sh, p);
912   d = (mpfr_uexp_t) (bx - MPFR_GET_EXP(c));
913 
914   /* printf ("New add1sp with diff=%lu\n", (unsigned long) d); */
915 
916   if (d == 0)
917     {
918       /* d==0 */
919       /* mpfr_print_mant_binary("C= ", MPFR_MANT(c), p); */
920       /* mpfr_print_mant_binary("B= ", MPFR_MANT(b), p); */
921       bx++;                                /* exp + 1 */
922       limb = mpn_add_n (ap, MPFR_MANT(b), MPFR_MANT(c), n);
923       /* mpfr_print_mant_binary("A= ", ap, p); */
924       MPFR_ASSERTD(limb != 0);             /* There must be a carry */
925       rb = ap[0] & (MPFR_LIMB_ONE << sh);  /* Get round bit (sb=0) */
926       mpn_rshift (ap, ap, n, 1);           /* Shift mantissa A */
927       ap[n-1] |= MPFR_LIMB_HIGHBIT;        /* Set MSB */
928       ap[0]   &= ~MPFR_LIMB_MASK(sh);      /* Clear round bit */
929 
930 
931       /* Fast track for faithful rounding: since b and c have the same
932          precision and the same exponent, the neglected value is either
933          0 or 1/2 ulp(a), thus returning a is fine. */
934       if (rnd_mode == MPFR_RNDF)
935         { inexact = 0; goto set_exponent; }
936 
937       if (rb == 0) /* Check exact case */
938         { inexact = 0; goto set_exponent; }
939 
940       /* Zero: Truncate
941          Nearest: Even Rule => truncate or add 1
942          Away: Add 1 */
943       if (MPFR_LIKELY (rnd_mode == MPFR_RNDN))
944         {
945           if ((ap[0] & (MPFR_LIMB_ONE << sh)) == 0)
946             { inexact = -1; goto set_exponent; }
947           else
948             goto add_one_ulp;
949         }
950       MPFR_UPDATE_RND_MODE(rnd_mode, neg);
951       if (rnd_mode == MPFR_RNDZ)
952         { inexact = -1; goto set_exponent; }
953       else
954         goto add_one_ulp;
955     }
956   else if (MPFR_UNLIKELY (d >= p))
957     {
958       /* fast track for RNDF */
959       if (MPFR_LIKELY(rnd_mode == MPFR_RNDF))
960         goto copy_set_exponent;
961 
962       if (MPFR_LIKELY (d > p))
963         {
964           /* d > p : Copy B in A */
965           /* Away:    Add 1
966              Nearest: Trunc
967              Zero:    Trunc */
968           if (rnd_mode == MPFR_RNDN || MPFR_IS_LIKE_RNDZ (rnd_mode, neg))
969             {
970             copy_set_exponent:
971               if (a != b)
972                 MPN_COPY (ap, MPFR_MANT(b), n);
973               inexact = -1;
974               goto set_exponent;
975             }
976           else
977             {
978             copy_add_one_ulp:
979               if (a != b)
980                 MPN_COPY (ap, MPFR_MANT(b), n);
981               goto add_one_ulp;
982             }
983         }
984       else
985         {
986           /* d==p : Copy B in A */
987           /* Away:     Add 1
988              Nearest:  Even Rule if C is a power of 2, else Add 1
989              Zero:     Trunc */
990           if (rnd_mode == MPFR_RNDN)
991             {
992               /* Check if C was a power of 2 */
993               if (mpfr_powerof2_raw (c) &&
994                   ((MPFR_MANT (b))[0] & (MPFR_LIMB_ONE << sh)) == 0)
995                 goto copy_set_exponent;
996               /* Not a Power of 2 */
997               goto copy_add_one_ulp;
998             }
999           else if (MPFR_IS_LIKE_RNDZ (rnd_mode, neg))
1000             goto copy_set_exponent;
1001           else
1002             goto copy_add_one_ulp;
1003         }
1004     }
1005   else /* 0 < d < p */
1006     {
1007       mp_limb_t mask = ~MPFR_LIMB_MASK(sh);
1008 
1009       /* General case: 1 <= d < p */
1010 
1011       limb = mpfr_addrsh (ap, MPFR_MANT(b), MPFR_MANT(c), n, d, &sb);
1012       /* the most significant bit of sb contains what would be the most
1013          significant bit of ap[-1], and the remaining bits of sb are 0
1014          iff the remaining bits of ap[-1], ap[-2], ... are all zero */
1015 
1016       if (sh > 0)
1017         {
1018           /* The round bit and a part of the sticky bit are in ap[0]. */
1019           rb = (ap[0] & (MPFR_LIMB_ONE << (sh - 1)));
1020           sb |= ap[0] & MPFR_LIMB_MASK (sh - 1);
1021         }
1022       else
1023         {
1024           /* The round bit and possibly a part of the sticky bit are
1025              in sb. */
1026           rb = sb & MPFR_LIMB_HIGHBIT;
1027           sb &= ~MPFR_LIMB_HIGHBIT;
1028         }
1029 
1030       ap[0] &= mask;
1031 
1032       /* Check for carry out */
1033       if (MPFR_UNLIKELY (limb != 0))
1034         {
1035           limb = ap[0] & (MPFR_LIMB_ONE << sh); /* Get LSB (will be new rb) */
1036           mpn_rshift (ap, ap, n, 1);            /* Shift significand */
1037           bx++;                                 /* Increase exponent */
1038           ap[n-1] |= MPFR_LIMB_HIGHBIT;         /* Set MSB */
1039           ap[0]   &= mask;                      /* Clear LSB */
1040           sb      |= rb;                        /* Update sb */
1041           rb      = limb;                       /* New rb */
1042           /* printf ("(Overflow) rb=%lu sb=%lu\n",
1043              (unsigned long) rb, (unsigned long) sb);
1044              mpfr_print_mant_binary ("Add=  ", ap, p); */
1045         }
1046 
1047       /* Round:
1048           Zero: Truncate but could be exact.
1049           Away: Add 1 if rb or sb !=0
1050           Nearest: Truncate but could be exact if sb==0
1051                    Add 1 if rb !=0,
1052                    Even rule else */
1053       if (MPFR_LIKELY(rnd_mode == MPFR_RNDF))
1054         { inexact = 0; goto set_exponent; }
1055       else if (rnd_mode == MPFR_RNDN)
1056         {
1057           inexact = - (sb != 0);
1058           if (rb == 0)
1059             goto set_exponent;
1060           else if (MPFR_UNLIKELY (sb == 0) &&
1061                    (ap[0] & (MPFR_LIMB_ONE << sh)) == 0)
1062             { inexact = -1; goto set_exponent; }
1063           else
1064             goto add_one_ulp;
1065         }
1066       MPFR_UPDATE_RND_MODE(rnd_mode, neg);
1067       inexact = - (rb != 0 || sb != 0);
1068       if (rnd_mode == MPFR_RNDZ || inexact == 0)
1069         goto set_exponent;
1070       else
1071         goto add_one_ulp;
1072     }
1073   MPFR_RET_NEVER_GO_HERE();
1074 
1075  add_one_ulp:
1076   /* add one unit in last place to a */
1077   /* printf("AddOneUlp\n"); */
1078   if (MPFR_UNLIKELY (mpn_add_1 (ap, ap, n, MPFR_LIMB_ONE << sh)))
1079     {
1080       /* Case 100000x0 = 0x1111x1 + 1*/
1081       /* printf("Pow of 2\n"); */
1082       bx++;
1083       ap[n-1] = MPFR_LIMB_HIGHBIT;
1084     }
1085   inexact = 1;
1086 
1087  set_exponent:
1088   if (MPFR_UNLIKELY(bx > __gmpfr_emax)) /* Check for overflow */
1089     {
1090       /* printf("Overflow\n"); */
1091       return mpfr_overflow (a, rnd_mode, MPFR_SIGN(a));
1092     }
1093   MPFR_SET_EXP (a, bx);
1094 
1095   MPFR_RET (inexact * MPFR_INT_SIGN (a));
1096 }
1097