1 /* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_prec_round,
2 mpfr_can_round, mpfr_can_round_raw -- various rounding functions
4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
22 MA 02110-1301, USA. */
24 #include "mpfr-impl.h"
26 #define mpfr_round_raw_generic mpfr_round_raw
29 #include "round_raw_generic.c"
31 #define mpfr_round_raw_generic mpfr_round_raw_2
34 #include "round_raw_generic.c"
36 /* Seems to be unused. Remove comment to implement it.
37 #define mpfr_round_raw_generic mpfr_round_raw_3
40 #include "round_raw_generic.c"
43 #define mpfr_round_raw_generic mpfr_round_raw_4
46 #include "round_raw_generic.c"
49 mpfr_prec_round (mpfr_ptr x
, mp_prec_t prec
, mp_rnd_t rnd_mode
)
54 MPFR_TMP_DECL(marker
);
56 MPFR_ASSERTN(prec
>= MPFR_PREC_MIN
&& prec
<= MPFR_PREC_MAX
);
58 nw
= 1 + (prec
- 1) / BITS_PER_MP_LIMB
; /* needed allocated limbs */
60 /* check if x has enough allocated space for the mantissa */
61 ow
= MPFR_GET_ALLOC_SIZE(x
);
64 /* Realloc mantissa */
65 mp_ptr tmpx
= (mp_ptr
) (*__gmp_reallocate_func
)
66 (MPFR_GET_REAL_PTR(x
), MPFR_MALLOC_SIZE(ow
), MPFR_MALLOC_SIZE(nw
));
67 MPFR_SET_MANT_PTR(x
, tmpx
); /* mant ptr must be set before alloc size */
68 MPFR_SET_ALLOC_SIZE(x
, nw
); /* new number of allocated limbs */
71 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x
) ))
73 MPFR_PREC(x
) = prec
; /* Special value: need to set prec */
76 MPFR_ASSERTD(MPFR_IS_INF(x
) || MPFR_IS_ZERO(x
));
77 return 0; /* infinity and zero are exact */
80 /* x is a non-zero real number */
82 MPFR_TMP_MARK(marker
);
83 tmp
= (mp_limb_t
*) MPFR_TMP_ALLOC (nw
* BYTES_PER_MP_LIMB
);
85 carry
= mpfr_round_raw (tmp
, xp
, MPFR_PREC(x
), MPFR_IS_NEG(x
),
86 prec
, rnd_mode
, &inexact
);
89 if (MPFR_UNLIKELY(carry
))
91 mp_exp_t exp
= MPFR_EXP (x
);
93 if (MPFR_UNLIKELY(exp
== __gmpfr_emax
))
94 (void) mpfr_overflow(x
, rnd_mode
, MPFR_SIGN(x
));
97 MPFR_ASSERTD (exp
< __gmpfr_emax
);
98 MPFR_SET_EXP (x
, exp
+ 1);
99 xp
[nw
- 1] = MPFR_LIMB_HIGHBIT
;
101 MPN_ZERO(xp
, nw
- 1);
105 MPN_COPY(xp
, tmp
, nw
);
107 MPFR_TMP_FREE(marker
);
111 /* assumption: BITS_PER_MP_LIMB is a power of 2 */
113 /* assuming b is an approximation to x in direction rnd1 with error at
114 most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly
115 x to precision prec with direction rnd2, and 0 otherwise.
121 mpfr_can_round (mpfr_srcptr b
, mp_exp_t err
, mp_rnd_t rnd1
,
122 mp_rnd_t rnd2
, mp_prec_t prec
)
124 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b
)))
125 return 0; /* We cannot round if Zero, Nan or Inf */
127 return mpfr_can_round_raw(MPFR_MANT(b
), MPFR_LIMB_SIZE(b
),
128 MPFR_SIGN(b
), err
, rnd1
, rnd2
, prec
);
132 mpfr_can_round_raw (const mp_limb_t
*bp
, mp_size_t bn
, int neg
, mp_exp_t err0
,
133 mp_rnd_t rnd1
, mp_rnd_t rnd2
, mp_prec_t prec
)
140 MPFR_TMP_DECL(marker
);
142 if (MPFR_UNLIKELY(err0
< 0 || (mp_exp_unsigned_t
) err0
<= prec
))
143 return 0; /* can't round */
144 else if (MPFR_UNLIKELY (prec
> (mp_prec_t
) bn
* BITS_PER_MP_LIMB
))
145 { /* then ulp(b) < precision < error */
146 return rnd2
== GMP_RNDN
&& (mp_exp_unsigned_t
) err0
- 2 >= prec
;
147 /* can round only in rounding to the nearest and err0 >= prec + 2 */
150 MPFR_ASSERT_SIGN(neg
);
151 neg
= MPFR_IS_NEG_SIGN(neg
);
153 /* if the error is smaller than ulp(b), then anyway it will propagate
155 err
= ((mp_exp_unsigned_t
) err0
> (mp_prec_t
) bn
* BITS_PER_MP_LIMB
) ?
156 (mp_prec_t
) bn
* BITS_PER_MP_LIMB
: (mp_prec_t
) err0
;
158 /* warning: if k = m*BITS_PER_MP_LIMB, consider limb m-1 and not m */
159 k
= (err
- 1) / BITS_PER_MP_LIMB
;
160 MPFR_UNSIGNED_MINUS_MODULO(s
, err
);
162 /* the error corresponds to bit s in limb k, the most significant limb
164 k1
= (prec
- 1) / BITS_PER_MP_LIMB
;
165 MPFR_UNSIGNED_MINUS_MODULO(s1
, prec
);
167 /* the last significant bit is bit s1 in limb k1 */
169 /* don't need to consider the k1 most significant limbs */
172 prec
-= (mp_prec_t
) k1
* BITS_PER_MP_LIMB
;
173 /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
174 change bp[bn-1] >> s1, then we can round */
176 MPFR_TMP_MARK(marker
);
178 k
++; /* since we work with k+1 everywhere */
179 tmp
= (mp_limb_t
*) MPFR_TMP_ALLOC(tn
* BYTES_PER_MP_LIMB
);
181 MPN_COPY (tmp
, bp
, bn
- k
);
183 MPFR_ASSERTD (k
> 0);
185 /* Transform RNDD and RNDU to Zero / Away */
186 MPFR_ASSERTD((neg
== 0) || (neg
==1));
187 if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd1
, neg
))
194 cc
= (bp
[bn
- 1] >> s1
) & 1;
195 cc
^= mpfr_round_raw2(bp
, bn
, neg
, rnd2
, prec
);
196 /* now round b +/- 2^(MPFR_EXP(b)-err) */
197 cc2
= mpn_add_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
200 /* Round to nearest */
201 /* first round b+2^(MPFR_EXP(b)-err) */
202 cc
= mpn_add_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
203 cc
= (tmp
[bn
- 1] >> s1
) & 1; /* gives 0 when cc=1 */
204 cc
^= mpfr_round_raw2 (tmp
, bn
, neg
, rnd2
, prec
);
205 /* now round b-2^(MPFR_EXP(b)-err) */
206 cc2
= mpn_sub_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
210 cc
= (bp
[bn
- 1] >> s1
) & 1;
211 cc
^= mpfr_round_raw2(bp
, bn
, neg
, rnd2
, prec
);
212 /* now round b +/- 2^(MPFR_EXP(b)-err) */
213 cc2
= mpn_sub_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
219 MPFR_TMP_FREE(marker
);
223 cc2
= (tmp
[bn
- 1] >> s1
) & 1;
224 cc2
^= mpfr_round_raw2 (tmp
, bn
, neg
, rnd2
, prec
);
226 MPFR_TMP_FREE(marker
);