1 /* mpfr_round_raw_generic -- Generic rounding function
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 # error "ERROR: flag must be defined (0 / 1)"
27 # error "ERROR: use_enexp must be defined (0 / 1)"
29 #ifndef mpfr_round_raw_generic
30 # error "ERROR: mpfr_round_raw_generic must be defined"
34 * If flag = 0, puts in y the value of xp (with precision xprec and
35 * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and
36 * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity
37 * (i.e. *xp != 0). In that case, the return value is a possible carry
38 * (0 or 1) that may happen during the rounding, in which case the result
41 * If inexp != NULL, put in *inexp the inexact flag of the rounding (0, 1, -1).
42 * In case of even rounding when rnd = MPFR_RNDN, put MPFR_EVEN_INEX (2) or
43 * -MPFR_EVEN_INEX (-2) in *inexp.
45 * If flag = 1, just returns whether one should add 1 or not for rounding.
47 * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal
48 * to 1. In this case, the even rounding is done away from 0, which is
49 * a natural generalization. Indeed, a number with 1-bit precision can
50 * be seen as a denormalized number with more precision.
54 mpfr_round_raw_generic(
58 const mp_limb_t
*xp
, mpfr_prec_t xprec
,
59 int neg
, mpfr_prec_t yprec
, mpfr_rnd_t rnd_mode
66 mp_limb_t himask
, lomask
, sb
;
76 MPFR_ASSERTD(inexp
!= ((int*) 0));
77 MPFR_ASSERTD(neg
== 0 || neg
== 1);
79 if (flag
&& !use_inexp
&&
80 (xprec
<= yprec
|| MPFR_IS_LIKE_RNDZ (rnd_mode
, neg
)))
83 xsize
= MPFR_PREC2LIMBS (xprec
);
84 nw
= yprec
/ GMP_NUMB_BITS
;
85 rw
= yprec
& (GMP_NUMB_BITS
- 1);
87 if (MPFR_UNLIKELY(xprec
<= yprec
))
88 { /* No rounding is necessary. */
89 /* if yp=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */
92 MPFR_ASSERTD(nw
>= 1);
93 MPFR_ASSERTD(nw
>= xsize
);
97 MPN_COPY_DECR(yp
+ (nw
- xsize
), xp
, xsize
);
98 MPN_ZERO(yp
, nw
- xsize
);
103 if (use_inexp
|| !MPFR_IS_LIKE_RNDZ(rnd_mode
, neg
))
105 mp_size_t k
= xsize
- nw
- 1;
110 lomask
= MPFR_LIMB_MASK (GMP_NUMB_BITS
- rw
);
115 lomask
= ~(mp_limb_t
) 0;
116 himask
= ~(mp_limb_t
) 0;
118 MPFR_ASSERTD(k
>= 0);
119 sb
= xp
[k
] & lomask
; /* First non-significant bits */
120 /* Rounding to nearest ? */
121 if (MPFR_LIKELY( rnd_mode
== MPFR_RNDN
) )
123 /* Rounding to nearest */
124 mp_limb_t rbmask
= MPFR_LIMB_ONE
<< (GMP_NUMB_BITS
- 1 - rw
);
125 if (sb
& rbmask
) /* rounding bit */
126 sb
&= ~rbmask
; /* it is 1, clear it */
129 /* Rounding bit is 0, behave like rounding to 0 */
132 while (MPFR_UNLIKELY(sb
== 0) && k
> 0)
134 /* rounding to nearest, with rounding bit = 1 */
135 if (MPFR_UNLIKELY(sb
== 0)) /* Even rounding. */
137 /* sb == 0 && rnd_mode == MPFR_RNDN */
138 sb
= xp
[xsize
- nw
] & (himask
^ (himask
<< 1));
142 *inexp
= 2*MPFR_EVEN_INEX
*neg
-MPFR_EVEN_INEX
;
143 /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX;*/
144 /* Since neg = 0 or 1 and sb=0*/
146 return 0 /*sb != 0 && rnd_mode != MPFR_RNDZ */;
148 MPN_COPY_INCR(yp
, xp
+ xsize
- nw
, nw
);
155 /* sb != 0 && rnd_mode == MPFR_RNDN */
157 *inexp
= MPFR_EVEN_INEX
-2*MPFR_EVEN_INEX
*neg
;
158 /*((neg!=0)^(sb!=0))? MPFR_EVEN_INEX : -MPFR_EVEN_INEX; */
159 /*Since neg= 0 or 1 and sb != 0 */
160 goto rnd_RNDN_add_one_ulp
;
163 else /* sb != 0 && rnd_mode == MPFR_RNDN*/
166 /* *inexp = (neg == 0) ? 1 : -1; but since neg = 0 or 1 */
168 rnd_RNDN_add_one_ulp
:
170 return 1; /*sb != 0 && rnd_mode != MPFR_RNDZ;*/
172 carry
= mpn_add_1 (yp
, xp
+ xsize
- nw
, nw
,
174 MPFR_LIMB_ONE
<< (GMP_NUMB_BITS
- rw
)
181 /* Rounding to Zero ? */
182 else if (MPFR_IS_LIKE_RNDZ(rnd_mode
, neg
))
184 /* rnd_mode == MPFR_RNDZ */
186 while (MPFR_UNLIKELY(sb
== 0) && k
> 0)
189 /* rnd_mode == MPFR_RNDZ and neg = 0 or 1 */
190 /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
191 *inexp
= MPFR_UNLIKELY(sb
== 0) ? 0 : (2*neg
-1);
193 return 0; /*sb != 0 && rnd_mode != MPFR_RNDZ;*/
195 MPN_COPY_INCR(yp
, xp
+ xsize
- nw
, nw
);
202 /* rnd_mode = Away */
203 while (MPFR_UNLIKELY(sb
== 0) && k
> 0)
205 if (MPFR_UNLIKELY(sb
== 0))
207 /* sb = 0 && rnd_mode != MPFR_RNDZ */
209 /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
214 MPN_COPY_INCR(yp
, xp
+ xsize
- nw
, nw
);
221 /* sb != 0 && rnd_mode != MPFR_RNDZ */
223 /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
228 carry
= mpn_add_1(yp
, xp
+ xsize
- nw
, nw
,
229 rw
? MPFR_LIMB_ONE
<< (GMP_NUMB_BITS
- rw
)
239 /* Roundind mode = Zero / No inexact flag */
241 return 0 /*sb != 0 && rnd_mode != MPFR_RNDZ*/;
246 himask
= ~MPFR_LIMB_MASK (GMP_NUMB_BITS
- rw
);
249 himask
= ~(mp_limb_t
) 0;
250 MPN_COPY_INCR(yp
, xp
+ xsize
- nw
, nw
);
259 #undef mpfr_round_raw_generic