1 /* mpfr_sub1sp -- internal function to perform a "real" substraction
2 All the op must have the same precision
4 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
27 /* Check if we have to check the result of mpfr_sub1sp with mpfr_sub1 */
31 int mpfr_sub1sp2 (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mpfr_rnd_t rnd_mode
);
32 int mpfr_sub1sp (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mpfr_rnd_t rnd_mode
)
34 mpfr_t tmpa
, tmpb
, tmpc
;
35 int inexb
, inexc
, inexact
, inexact2
;
37 mpfr_init2 (tmpa
, MPFR_PREC (a
));
38 mpfr_init2 (tmpb
, MPFR_PREC (b
));
39 mpfr_init2 (tmpc
, MPFR_PREC (c
));
41 inexb
= mpfr_set (tmpb
, b
, MPFR_RNDN
);
42 MPFR_ASSERTN (inexb
== 0);
44 inexc
= mpfr_set (tmpc
, c
, MPFR_RNDN
);
45 MPFR_ASSERTN (inexc
== 0);
47 inexact2
= mpfr_sub1 (tmpa
, tmpb
, tmpc
, rnd_mode
);
48 inexact
= mpfr_sub1sp2(a
, b
, c
, rnd_mode
);
50 if (mpfr_cmp (tmpa
, a
) || inexact
!= inexact2
)
52 fprintf (stderr
, "sub1 & sub1sp return different values for %s\n"
53 "Prec_a = %lu, Prec_b = %lu, Prec_c = %lu\nB = ",
54 mpfr_print_rnd_mode (rnd_mode
), (unsigned long) MPFR_PREC (a
),
55 (unsigned long) MPFR_PREC (b
), (unsigned long) MPFR_PREC (c
));
56 mpfr_fprint_binary (stderr
, tmpb
);
57 fprintf (stderr
, "\nC = ");
58 mpfr_fprint_binary (stderr
, tmpc
);
59 fprintf (stderr
, "\nSub1 : ");
60 mpfr_fprint_binary (stderr
, tmpa
);
61 fprintf (stderr
, "\nSub1sp: ");
62 mpfr_fprint_binary (stderr
, a
);
63 fprintf (stderr
, "\nInexact sp = %d | Inexact = %d\n",
67 mpfr_clears (tmpa
, tmpb
, tmpc
, (mpfr_ptr
) 0);
70 # define mpfr_sub1sp mpfr_sub1sp2
74 /* Debugging support */
79 # define DEBUG(x) /**/
85 compute sgn(b)*(|b| - |c|) if |b|>|c| else -sgn(b)*(|c| -|b|)
86 Returns 0 iff result is exact,
87 a negative value when the result is less than the exact value,
88 a positive value otherwise.
94 * Cp = -1 if calculated from c mantissa
95 * Cp = 0 if 0 from a or c
96 * Cp = 1 if calculated from a.
97 * C'p+1 = First bit not null or 0 if there isn't one
99 * Can't have Cp=-1 and C'p+1=1*/
102 * + if Cp=0 and C'p+1=0,1, Truncate.
103 * + if Cp=0 and C'p+1=-1, SubOneUlp
104 * + if Cp=-1, SubOneUlp
105 * + if Cp=1, AddOneUlp
106 * RND = MPFR_RNDA (Away)
107 * + if Cp=0 and C'p+1=0,-1, Truncate
108 * + if Cp=0 and C'p+1=1, AddOneUlp
109 * + if Cp=1, AddOneUlp
110 * + if Cp=-1, Truncate
112 * + if Cp=0, Truncate
113 * + if Cp=1 and C'p+1=1, AddOneUlp
114 * + if Cp=1 and C'p+1=-1, Truncate
115 * + if Cp=1 and C'p+1=0, Truncate if Ap-1=0, AddOneUlp else
116 * + if Cp=-1 and C'p+1=-1, SubOneUlp
117 * + if Cp=-1 and C'p+1=0, Truncate if Ap-1=0, SubOneUlp else
120 * If carry, then it is 11111111111 + 1 = 10000000000000
121 * ap[n-1]=MPFR_HIGHT_BIT
123 * If we lose one bit, it is 1000000000 - 1 = 0111111111111
124 * Then shift, and put as last bit x which is calculated
125 * according Cp, Cp-1 and rnd_mode.
127 * If it is a power of 2,
128 * we may have to suboneulp in some special cases.
130 * To simplify, we don't use Cp = 1.
135 mpfr_sub1sp (mpfr_ptr a
, mpfr_srcptr b
, mpfr_srcptr c
, mpfr_rnd_t rnd_mode
)
139 mpfr_prec_t p
, sh
, cnt
;
141 mp_limb_t
*ap
, *bp
, *cp
;
144 mp_limb_t bcp
,bcp1
; /* Cp and C'p+1 */
145 mp_limb_t bbcp
= (mp_limb_t
) -1, bbcp1
= (mp_limb_t
) -1; /* Cp+1 and C'p+2,
146 gcc claims that they might be used uninitialized. We fill them with invalid
147 values, which should produce a failure if so. See README.dev file. */
149 MPFR_TMP_DECL(marker
);
151 MPFR_TMP_MARK(marker
);
153 MPFR_ASSERTD(MPFR_PREC(a
) == MPFR_PREC(b
) && MPFR_PREC(b
) == MPFR_PREC(c
));
154 MPFR_ASSERTD(MPFR_IS_PURE_FP(b
));
155 MPFR_ASSERTD(MPFR_IS_PURE_FP(c
));
157 /* Read prec and num of limbs */
159 n
= MPFR_PREC2LIMBS (p
);
161 /* Fast cmp of |b| and |c|*/
162 bx
= MPFR_GET_EXP (b
);
163 cx
= MPFR_GET_EXP (c
);
164 if (MPFR_UNLIKELY(bx
== cx
))
167 /* Check mantissa since exponent are equals */
170 while (k
>=0 && MPFR_UNLIKELY(bp
[k
] == cp
[k
]))
172 if (MPFR_UNLIKELY(k
< 0))
175 /* Return exact number 0 */
176 if (rnd_mode
== MPFR_RNDD
)
183 else if (bp
[k
] > cp
[k
])
187 MPFR_ASSERTD(bp
[k
]<cp
[k
]);
191 else if (MPFR_UNLIKELY(bx
< cx
))
193 /* Swap b and c and set sign */
197 MPFR_SET_OPPOSITE_SIGN(a
,b
);
199 tx
= bx
; bx
= cx
; cx
= tx
;
205 MPFR_SET_SAME_SIGN(a
,b
);
209 MPFR_ASSERTD(bx
>= cx
);
210 d
= (mpfr_uexp_t
) bx
- cx
;
211 DEBUG (printf ("New with diff=%lu\n", (unsigned long) d
));
213 if (MPFR_UNLIKELY(d
<= 1))
215 if (MPFR_LIKELY(d
< 1))
218 <-- c --> : exact sub */
220 mpn_sub_n (ap
, MPFR_MANT(b
), MPFR_MANT(c
), n
);
224 if (MPFR_LIKELY(limb
))
226 /* First limb is not zero. */
227 count_leading_zeros(cnt
, limb
);
228 /* cnt could be == 0 <= SubD1Lose */
229 if (MPFR_LIKELY(cnt
))
231 mpn_lshift(ap
, ap
, n
, cnt
); /* Normalize number */
232 bx
-= cnt
; /* Update final expo */
234 /* Last limb should be ok */
235 MPFR_ASSERTD(!(ap
[0] & MPFR_LIMB_MASK((unsigned int) (-p
)
240 /* First limb is zero */
241 mp_size_t k
= n
-1, len
;
242 /* Find the first limb not equal to zero.
243 FIXME:It is assume it exists (since |b| > |c| and same prec)*/
246 MPFR_ASSERTD( k
> 0 );
250 MPFR_ASSERTD(limb
!= 0);
251 count_leading_zeros(cnt
, limb
);
253 len
= n
- k
; /* Number of last limb */
254 MPFR_ASSERTD(k
>= 0);
255 if (MPFR_LIKELY(cnt
))
256 mpn_lshift(ap
+len
, ap
, k
, cnt
); /* Normalize the High Limb*/
259 /* Must use DECR since src and dest may overlap & dest>=src*/
260 MPN_COPY_DECR(ap
+len
, ap
, k
);
262 MPN_ZERO(ap
, len
); /* Zeroing the last limbs */
263 bx
-= cnt
+ len
*GMP_NUMB_BITS
; /* Update Expo */
264 /* Last limb should be ok */
265 MPFR_ASSERTD(!(ap
[len
]&MPFR_LIMB_MASK((unsigned int) (-p
)
268 /* Check expo underflow */
269 if (MPFR_UNLIKELY(bx
< __gmpfr_emin
))
271 MPFR_TMP_FREE(marker
);
273 DEBUG( printf("(D==0 Underflow)\n") );
274 if (rnd_mode
== MPFR_RNDN
&&
275 (bx
< __gmpfr_emin
- 1 ||
276 (/*inexact >= 0 &&*/ mpfr_powerof2_raw (a
))))
277 rnd_mode
= MPFR_RNDZ
;
278 return mpfr_underflow (a
, rnd_mode
, MPFR_SIGN(a
));
280 MPFR_SET_EXP (a
, bx
);
281 /* No rounding is necessary since the result is exact */
282 MPFR_ASSERTD(ap
[n
-1] > ~ap
[n
-1]);
283 MPFR_TMP_FREE(marker
);
286 else /* if (d == 1) */
292 MPFR_UNSIGNED_MINUS_MODULO(sh
, p
);
293 /* If we lose at least one bit, compute 2*b-c (Exact)
294 * else compute b-c/2 */
298 limb
= bp
[k
] - cp
[k
]/2;
299 if (limb
> MPFR_LIMB_HIGHBIT
)
301 /* We can't lose precision: compute b-c/2 */
302 /* Shift c in the allocated temporary block */
304 c0
= cp
[0] & (MPFR_LIMB_ONE
<<sh
);
305 cp
= MPFR_TMP_LIMBS_ALLOC (n
);
306 mpn_rshift(cp
, MPFR_MANT(c
), n
, 1);
307 if (MPFR_LIKELY(c0
== 0))
309 /* Result is exact: no need of rounding! */
311 mpn_sub_n (ap
, bp
, cp
, n
);
312 MPFR_SET_EXP(a
, bx
); /* No expo overflow! */
313 /* No truncate or normalize is needed */
314 MPFR_ASSERTD(ap
[n
-1] > ~ap
[n
-1]);
315 /* No rounding is necessary since the result is exact */
316 MPFR_TMP_FREE(marker
);
320 mask
= ~MPFR_LIMB_MASK(sh
);
321 cp
[0] &= mask
; /* Delete last bit of c */
322 mpn_sub_n (ap
, bp
, cp
, n
);
323 MPFR_SET_EXP(a
, bx
); /* No expo overflow! */
324 MPFR_ASSERTD( !(ap
[0] & ~mask
) ); /* Check last bits */
325 /* No normalize is needed */
326 MPFR_ASSERTD(ap
[n
-1] > ~ap
[n
-1]);
327 /* Rounding is necessary since c0 = 1*/
328 /* Cp =-1 and C'p+1=0 */
330 if (MPFR_LIKELY(rnd_mode
== MPFR_RNDN
))
332 /* Even Rule apply: Check Ap-1 */
333 if (MPFR_LIKELY( (ap
[0] & (MPFR_LIMB_ONE
<<sh
)) == 0) )
338 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(a
));
339 if (rnd_mode
== MPFR_RNDZ
)
344 else if (MPFR_LIKELY(limb
< MPFR_LIMB_HIGHBIT
))
346 /* We lose at least one bit of prec */
347 /* Calcul of 2*b-c (Exact) */
348 /* Shift b in the allocated temporary block */
350 bp
= MPFR_TMP_LIMBS_ALLOC (n
);
351 mpn_lshift (bp
, MPFR_MANT(b
), n
, 1);
353 mpn_sub_n (ap
, bp
, cp
, n
);
359 /* Case: limb = 100000000000 */
360 /* Check while b[k] == c'[k] (C' is C shifted by 1) */
361 /* If b[k]<c'[k] => We lose at least one bit*/
362 /* If b[k]>c'[k] => We don't lose any bit */
363 /* If k==-1 => We don't lose any bit
364 AND the result is 100000000000 0000000000 00000000000 */
367 carry
= cp
[k
]&MPFR_LIMB_ONE
;
370 bp
[k
]==(carry
=cp
[k
]/2+(carry
<<(GMP_NUMB_BITS
-1))));
371 if (MPFR_UNLIKELY(k
<0))
373 /*If carry then (sh==0 and Virtual c'[-1] > Virtual b[-1]) */
374 if (MPFR_UNLIKELY(carry
)) /* carry = cp[0]&MPFR_LIMB_ONE */
376 /* FIXME: Can be faster? */
377 MPFR_ASSERTD(sh
== 0);
380 /* Result is a power of 2 */
383 ap
[n
-1] = MPFR_LIMB_HIGHBIT
;
384 MPFR_SET_EXP (a
, bx
); /* No expo overflow! */
385 /* No Normalize is needed*/
386 /* No Rounding is needed */
387 MPFR_TMP_FREE (marker
);
390 /* carry = cp[k]/2+(cp[k-1]&1)<<(GMP_NUMB_BITS-1) = c'[k]*/
391 else if (bp
[k
] > carry
)
395 MPFR_ASSERTD(bp
[k
]<carry
);
401 else if (MPFR_UNLIKELY(d
>= p
))
404 MPFR_UNSIGNED_MINUS_MODULO(sh
, p
);
405 /* We can't set A before since we use cp for rounding... */
406 /* Perform rounding: check if a=b or a=b-ulp(b) */
407 if (MPFR_UNLIKELY(d
== p
))
409 /* cp == -1 and c'p+1 = ? */
411 /* We need Cp+1 later for a very improbable case. */
412 bbcp
= (MPFR_MANT(c
)[n
-1] & (MPFR_LIMB_ONE
<<(GMP_NUMB_BITS
-2)));
413 /* We need also C'p+1 for an even more unprobable case... */
414 if (MPFR_LIKELY( bbcp
))
419 if (MPFR_UNLIKELY(cp
[n
-1] == MPFR_LIMB_HIGHBIT
))
424 } while (k
>=0 && cp
[k
]==0);
430 DEBUG( printf("(D=P) Cp=-1 Cp+1=%d C'p+1=%d \n", bbcp
!=0, bcp1
!=0) );
433 /* Even if src and dest overlap, it is ok using MPN_COPY */
434 if (MPFR_LIKELY(rnd_mode
== MPFR_RNDN
))
436 if (MPFR_UNLIKELY( bcp
&& bcp1
==0 ))
437 /* Cp=-1 and C'p+1=0: Even rule Apply! */
438 /* Check Ap-1 = Bp-1 */
439 if ((bp
[0] & (MPFR_LIMB_ONE
<<sh
)) == 0)
447 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(a
));
448 if (rnd_mode
== MPFR_RNDZ
)
461 /* Cp=0, Cp+1=-1 if d==p+1, C'p+1=-1 */
462 bcp
= 0; bbcp
= (d
==p
+1); bcp1
= 1;
463 DEBUG( printf("(D>P) Cp=%d Cp+1=%d C'p+1=%d\n", bcp
!=0,bbcp
!=0,bcp1
!=0) );
464 /* Need to compute C'p+2 if d==p+1 and if rnd_mode=NEAREST
465 (Because of a very improbable case) */
466 if (MPFR_UNLIKELY(d
==p
+1 && rnd_mode
==MPFR_RNDN
))
469 if (MPFR_UNLIKELY(cp
[n
-1] == MPFR_LIMB_HIGHBIT
))
474 } while (k
>=0 && cp
[k
]==0);
479 DEBUG( printf("(D>P) C'p+2=%d\n", bbcp1
!=0) );
481 /* Copy mantissa B in A */
482 MPN_COPY(ap
, MPFR_MANT(b
), n
);
484 if (MPFR_LIKELY(rnd_mode
== MPFR_RNDN
))
486 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(a
));
487 if (rnd_mode
== MPFR_RNDZ
)
489 else /* rnd_mode = AWAY */
499 /* General case: 2 <= d < p */
500 MPFR_UNSIGNED_MINUS_MODULO(sh
, p
);
501 cp
= MPFR_TMP_LIMBS_ALLOC (n
);
503 /* Shift c in temporary allocated place */
504 dm
= d
% GMP_NUMB_BITS
;
505 m
= d
/ GMP_NUMB_BITS
;
506 if (MPFR_UNLIKELY(dm
== 0))
508 /* dm = 0 and m > 0: Just copy */
510 MPN_COPY(cp
, MPFR_MANT(c
)+m
, n
-m
);
513 else if (MPFR_LIKELY(m
== 0))
515 /* dm >=2 and m == 0: just shift */
516 MPFR_ASSERTD(dm
>= 2);
517 mpn_rshift(cp
, MPFR_MANT(c
), n
, dm
);
521 /* dm > 0 and m > 0: shift and zero */
522 mpn_rshift(cp
, MPFR_MANT(c
)+m
, n
-m
, dm
);
526 DEBUG( mpfr_print_mant_binary("Before", MPFR_MANT(c
), p
) );
527 DEBUG( mpfr_print_mant_binary("B= ", MPFR_MANT(b
), p
) );
528 DEBUG( mpfr_print_mant_binary("After ", cp
, p
) );
530 /* Compute bcp=Cp and bcp1=C'p+1 */
533 /* Try to compute them from C' rather than C (FIXME: Faster?) */
534 bcp
= (cp
[0] & (MPFR_LIMB_ONE
<<(sh
-1))) ;
535 if (MPFR_LIKELY( cp
[0] & MPFR_LIMB_MASK(sh
-1) ))
539 /* We can't compute C'p+1 from C'. Compute it from C */
540 /* Start from bit x=p-d+sh in mantissa C
541 (+sh since we have already looked sh bits in C'!) */
542 mpfr_prec_t x
= p
-d
+sh
-1;
543 if (MPFR_LIKELY(x
>p
))
544 /* We are already looked at all the bits of c, so C'p+1 = 0*/
548 mp_limb_t
*tp
= MPFR_MANT(c
);
549 mp_size_t kx
= n
-1 - (x
/ GMP_NUMB_BITS
);
550 mpfr_prec_t sx
= GMP_NUMB_BITS
-1-(x
%GMP_NUMB_BITS
);
551 DEBUG (printf ("(First) x=%lu Kx=%ld Sx=%lu\n",
552 (unsigned long) x
, (long) kx
,
553 (unsigned long) sx
));
554 /* Looks at the last bits of limb kx (if sx=0 does nothing)*/
555 if (tp
[kx
] & MPFR_LIMB_MASK(sx
))
560 /*If sx==0, tp[kx] hasn't been checked*/
563 } while (kx
>=0 && tp
[kx
]==0);
571 /* Compute Cp and C'p+1 from C with sh=0 */
572 mp_limb_t
*tp
= MPFR_MANT(c
);
573 /* Start from bit x=p-d in mantissa C */
575 mp_size_t kx
= n
-1 - (x
/ GMP_NUMB_BITS
);
576 mpfr_prec_t sx
= GMP_NUMB_BITS
-1-(x
%GMP_NUMB_BITS
);
577 MPFR_ASSERTD(p
>= d
);
578 bcp
= (tp
[kx
] & (MPFR_LIMB_ONE
<<sx
));
579 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
580 if (tp
[kx
] & MPFR_LIMB_MASK(sx
))
584 /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
587 } while (kx
>=0 && tp
[kx
]==0);
591 DEBUG( printf("sh=%lu Cp=%d C'p+1=%d\n", sh
, bcp
!=0, bcp1
!=0) );
593 /* Check if we can lose a bit, and if so compute Cp+1 and C'p+2 */
595 if (MPFR_UNLIKELY((bp
[n
-1]-cp
[n
-1]) <= MPFR_LIMB_HIGHBIT
))
597 /* We can lose a bit so we precompute Cp+1 and C'p+2 */
598 /* Test for trivial case: since C'p+1=0, Cp+1=0 and C'p+2 =0 */
599 if (MPFR_LIKELY(bcp1
== 0))
606 /* We can lose a bit:
607 compute Cp+1 and C'p+2 from mantissa C */
608 mp_limb_t
*tp
= MPFR_MANT(c
);
609 /* Start from bit x=(p+1)-d in mantissa C */
610 mpfr_prec_t x
= p
+1-d
;
611 mp_size_t kx
= n
-1 - (x
/GMP_NUMB_BITS
);
612 mpfr_prec_t sx
= GMP_NUMB_BITS
-1-(x
%GMP_NUMB_BITS
);
614 DEBUG (printf ("(pre) x=%lu Kx=%ld Sx=%lu\n",
615 (unsigned long) x
, (long) kx
,
616 (unsigned long) sx
));
617 bbcp
= (tp
[kx
] & (MPFR_LIMB_ONE
<<sx
)) ;
618 /* Looks at the last bits of limb kx (If sx=0, does nothing)*/
619 /* If Cp+1=0, since C'p+1!=0, C'p+2=1 ! */
620 if (MPFR_LIKELY(bbcp
==0 || (tp
[kx
]&MPFR_LIMB_MASK(sx
))))
624 /*kx += (sx==0);*/ /*If sx==0, tp[kx] hasn't been checked*/
627 } while (kx
>=0 && tp
[kx
]==0);
629 DEBUG (printf ("(Pre) Scan done for %ld\n", (long) kx
));
631 } /*End of Bcp1 != 0*/
632 DEBUG( printf("(Pre) Cp+1=%d C'p+2=%d\n", bbcp
!=0, bbcp1
!=0) );
633 } /* End of "can lose a bit" */
635 /* Clean shifted C' */
636 mask
= ~MPFR_LIMB_MASK (sh
);
639 /* Substract the mantissa c from b in a */
641 mpn_sub_n (ap
, bp
, cp
, n
);
642 DEBUG( mpfr_print_mant_binary("Sub= ", ap
, p
) );
644 /* Normalize: we lose at max one bit*/
645 if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap
[n
-1]) == 0))
647 /* High bit is not set and we have to fix it! */
648 /* Ap >= 010000xxx001 */
649 mpn_lshift(ap
, ap
, n
, 1);
650 /* Ap >= 100000xxx010 */
651 if (MPFR_UNLIKELY(bcp
!=0)) /* Check if Cp = -1 */
652 /* Since Cp == -1, we have to substract one more */
654 mpn_sub_1(ap
, ap
, n
, MPFR_LIMB_ONE
<<sh
);
655 MPFR_ASSERTD(MPFR_LIMB_MSB(ap
[n
-1]) != 0);
657 /* Ap >= 10000xxx001 */
658 /* Final exponent -1 since we have shifted the mantissa */
660 /* Update bcp and bcp1 */
661 MPFR_ASSERTN(bbcp
!= (mp_limb_t
) -1);
662 MPFR_ASSERTN(bbcp1
!= (mp_limb_t
) -1);
665 /* We dont't have anymore a valid Cp+1!
666 But since Ap >= 100000xxx001, the final sub can't unnormalize!*/
668 MPFR_ASSERTD( !(ap
[0] & ~mask
) );
671 if (MPFR_LIKELY(rnd_mode
== MPFR_RNDN
))
673 if (MPFR_LIKELY(bcp
==0))
675 else if ((bcp1
) || ((ap
[0] & (MPFR_LIMB_ONE
<<sh
)) != 0))
681 /* Update rounding mode */
682 MPFR_UPDATE_RND_MODE(rnd_mode
, MPFR_IS_NEG(a
));
683 if (rnd_mode
== MPFR_RNDZ
&& (MPFR_LIKELY(bcp
|| bcp1
)))
687 MPFR_RET_NEVER_GO_HERE ();
689 /* Sub one ulp to the result */
691 mpn_sub_1 (ap
, ap
, n
, MPFR_LIMB_ONE
<< sh
);
692 /* Result should be smaller than exact value: inexact=-1 */
694 /* Check normalisation */
695 if (MPFR_UNLIKELY(MPFR_LIMB_MSB(ap
[n
-1]) == 0))
697 /* ap was a power of 2, and we lose a bit */
698 /* Now it is 0111111111111111111[00000 */
699 mpn_lshift(ap
, ap
, n
, 1);
701 /* And the lost bit x depends on Cp+1, and Cp */
702 /* Compute Cp+1 if it isn't already compute (ie d==1) */
703 /* FIXME: Is this case possible? */
704 if (MPFR_UNLIKELY(d
== 1))
706 DEBUG( printf("(SubOneUlp)Cp=%d, Cp+1=%d C'p+1=%d\n", bcp
!=0,bbcp
!=0,bcp1
!=0));
707 /* Compute the last bit (Since we have shifted the mantissa)
708 we need one more bit!*/
709 MPFR_ASSERTN(bbcp
!= (mp_limb_t
) -1);
710 if ( (rnd_mode
== MPFR_RNDZ
&& bcp
==0)
711 || (rnd_mode
==MPFR_RNDN
&& bbcp
==0)
712 || (bcp
&& bcp1
==0) ) /*Exact result*/
714 ap
[0] |= MPFR_LIMB_ONE
<<sh
;
715 if (rnd_mode
== MPFR_RNDN
)
717 DEBUG( printf("(SubOneUlp) Last bit set\n") );
719 /* Result could be exact if C'p+1 = 0 and rnd == Zero
720 since we have had one more bit to the result */
721 /* Fixme: rnd_mode == MPFR_RNDZ needed ? */
722 if (bcp1
==0 && rnd_mode
==MPFR_RNDZ
)
724 DEBUG( printf("(SubOneUlp) Exact result\n") );
732 /* Check if the result is an exact power of 2: 100000000000
733 in which cases, we could have to do sub_one_ulp due to some nasty reasons:
734 If Result is a Power of 2:
736 | If Cp=-1 and C'p+1 = 0, SubOneUlp and the result is EXACT.
737 If Cp=-1 and C'p+1 =-1, SubOneUlp and the result is above.
740 If Cp= 0 and Cp+1 =-1 and C'p+2=-1, SubOneUlp and the result is above
741 If cp=-1 and C'p+1 = 0, SubOneUlp and the result is exact.
743 X bit should always be set if SubOneUlp*/
744 if (MPFR_UNLIKELY(ap
[n
-1] == MPFR_LIMB_HIGHBIT
))
749 } while (k
>=0 && ap
[k
]==0);
750 if (MPFR_UNLIKELY(k
<0))
752 /* It is a power of 2! */
753 /* Compute Cp+1 if it isn't already compute (ie d==1) */
754 /* FIXME: Is this case possible? */
757 DEBUG( printf("(Truncate) Cp=%d, Cp+1=%d C'p+1=%d C'p+2=%d\n", \
758 bcp
!=0, bbcp
!=0, bcp1
!=0, bbcp1
!=0) );
759 MPFR_ASSERTN(bbcp
!= (mp_limb_t
) -1);
760 MPFR_ASSERTN((rnd_mode
!= MPFR_RNDN
) || (bcp
!= 0) || (bbcp
== 0) || (bbcp1
!= (mp_limb_t
) -1));
761 if (((rnd_mode
!= MPFR_RNDZ
) && bcp
)
763 ((rnd_mode
== MPFR_RNDN
) && (bcp
== 0) && (bbcp
) && (bbcp1
)))
765 DEBUG( printf("(Truncate) Do sub\n") );
766 mpn_sub_1 (ap
, ap
, n
, MPFR_LIMB_ONE
<< sh
);
767 mpn_lshift(ap
, ap
, n
, 1);
768 ap
[0] |= MPFR_LIMB_ONE
<<sh
;
770 /* FIXME: Explain why it works (or why not)... */
771 inexact
= (bcp1
== 0) ? 0 : (rnd_mode
==MPFR_RNDN
) ? -1 : 1;
777 /* Calcul of Inexact flag.*/
778 inexact
= MPFR_LIKELY(bcp
|| bcp1
) ? 1 : 0;
782 /* FIXME: Is this test really useful?
783 If d==0 : Exact case. This is never called.
784 if 1 < d < p : bx=MPFR_EXP(b) or MPFR_EXP(b)-1 > MPFR_EXP(c) > emin
785 if d == 1 : bx=MPFR_EXP(b). If we could lose any bits, the exact
786 normalisation is called.
787 if d >= p : bx=MPFR_EXP(b) >= MPFR_EXP(c) + p > emin
788 After SubOneUlp, we could have one bit less.
789 if 1 < d < p : bx >= MPFR_EXP(b)-2 >= MPFR_EXP(c) > emin
790 if d == 1 : bx >= MPFR_EXP(b)-1 = MPFR_EXP(c) > emin.
791 if d >= p : bx >= MPFR_EXP(b)-1 > emin since p>=2.
793 MPFR_ASSERTD( bx
>= __gmpfr_emin
);
795 if (MPFR_UNLIKELY(bx < __gmpfr_emin))
797 DEBUG( printf("(Final Underflow)\n") );
798 if (rnd_mode == MPFR_RNDN &&
799 (bx < __gmpfr_emin - 1 ||
800 (inexact >= 0 && mpfr_powerof2_raw (a))))
801 rnd_mode = MPFR_RNDZ;
802 MPFR_TMP_FREE(marker);
803 return mpfr_underflow (a, rnd_mode, MPFR_SIGN(a));
806 MPFR_SET_EXP (a
, bx
);
808 MPFR_TMP_FREE(marker
);
809 MPFR_RET (inexact
* MPFR_INT_SIGN (a
));