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