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-2015 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 #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
, mpfr_prec_t prec
, mpfr_rnd_t rnd_mode
)
54 MPFR_TMP_DECL(marker
);
56 MPFR_ASSERTN(prec
>= MPFR_PREC_MIN
&& prec
<= MPFR_PREC_MAX
);
58 nw
= MPFR_PREC2LIMBS (prec
); /* needed allocated limbs */
60 /* check if x has enough allocated space for the significand */
61 /* Get the number of limbs from the precision.
62 (Compatible with all allocation methods) */
63 ow
= MPFR_LIMB_SIZE (x
);
66 /* FIXME: Variable can't be created using custom allocation,
67 MPFR_DECL_INIT or GROUP_ALLOC: How to detect? */
68 ow
= MPFR_GET_ALLOC_SIZE(x
);
71 /* Realloc significand */
72 mpfr_limb_ptr tmpx
= (mpfr_limb_ptr
) (*__gmp_reallocate_func
)
73 (MPFR_GET_REAL_PTR(x
), MPFR_MALLOC_SIZE(ow
), MPFR_MALLOC_SIZE(nw
));
74 MPFR_SET_MANT_PTR(x
, tmpx
); /* mant ptr must be set
76 MPFR_SET_ALLOC_SIZE(x
, nw
); /* new number of allocated limbs */
80 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x
) ))
82 MPFR_PREC(x
) = prec
; /* Special value: need to set prec */
85 MPFR_ASSERTD(MPFR_IS_INF(x
) || MPFR_IS_ZERO(x
));
86 return 0; /* infinity and zero are exact */
89 /* x is a non-zero real number */
91 MPFR_TMP_MARK(marker
);
92 tmp
= MPFR_TMP_LIMBS_ALLOC (nw
);
94 carry
= mpfr_round_raw (tmp
, xp
, MPFR_PREC(x
), MPFR_IS_NEG(x
),
95 prec
, rnd_mode
, &inexact
);
98 if (MPFR_UNLIKELY(carry
))
100 mpfr_exp_t exp
= MPFR_EXP (x
);
102 if (MPFR_UNLIKELY(exp
== __gmpfr_emax
))
103 (void) mpfr_overflow(x
, rnd_mode
, MPFR_SIGN(x
));
106 MPFR_ASSERTD (exp
< __gmpfr_emax
);
107 MPFR_SET_EXP (x
, exp
+ 1);
108 xp
[nw
- 1] = MPFR_LIMB_HIGHBIT
;
110 MPN_ZERO(xp
, nw
- 1);
114 MPN_COPY(xp
, tmp
, nw
);
116 MPFR_TMP_FREE(marker
);
120 /* assumption: GMP_NUMB_BITS is a power of 2 */
122 /* assuming b is an approximation to x in direction rnd1 with error at
123 most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly
124 x to precision prec with direction rnd2, and 0 otherwise.
130 mpfr_can_round (mpfr_srcptr b
, mpfr_exp_t err
, mpfr_rnd_t rnd1
,
131 mpfr_rnd_t rnd2
, mpfr_prec_t prec
)
133 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b
)))
134 return 0; /* We cannot round if Zero, Nan or Inf */
136 return mpfr_can_round_raw (MPFR_MANT(b
), MPFR_LIMB_SIZE(b
),
137 MPFR_SIGN(b
), err
, rnd1
, rnd2
, prec
);
141 mpfr_can_round_raw (const mp_limb_t
*bp
, mp_size_t bn
, int neg
, mpfr_exp_t err0
,
142 mpfr_rnd_t rnd1
, mpfr_rnd_t rnd2
, mpfr_prec_t prec
)
149 MPFR_TMP_DECL(marker
);
151 if (MPFR_UNLIKELY(err0
< 0 || (mpfr_uexp_t
) err0
<= prec
))
152 return 0; /* can't round */
153 else if (MPFR_UNLIKELY (prec
> (mpfr_prec_t
) bn
* GMP_NUMB_BITS
))
154 { /* then ulp(b) < precision < error */
155 return rnd2
== MPFR_RNDN
&& (mpfr_uexp_t
) err0
- 2 >= prec
;
156 /* can round only in rounding to the nearest and err0 >= prec + 2 */
159 MPFR_ASSERT_SIGN(neg
);
160 neg
= MPFR_IS_NEG_SIGN(neg
);
162 /* if the error is smaller than ulp(b), then anyway it will propagate
164 err
= ((mpfr_uexp_t
) err0
> (mpfr_prec_t
) bn
* GMP_NUMB_BITS
) ?
165 (mpfr_prec_t
) bn
* GMP_NUMB_BITS
: (mpfr_prec_t
) err0
;
167 /* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */
168 k
= (err
- 1) / GMP_NUMB_BITS
;
169 MPFR_UNSIGNED_MINUS_MODULO(s
, err
);
170 /* the error corresponds to bit s in limb k, the most significant limb
173 k1
= (prec
- 1) / GMP_NUMB_BITS
;
174 MPFR_UNSIGNED_MINUS_MODULO(s1
, prec
);
175 /* the last significant bit is bit s1 in limb k1 */
177 /* don't need to consider the k1 most significant limbs */
180 prec
-= (mpfr_prec_t
) k1
* GMP_NUMB_BITS
;
182 /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
183 change bp[bn-1] >> s1, then we can round */
184 MPFR_TMP_MARK(marker
);
186 k
++; /* since we work with k+1 everywhere */
187 tmp
= MPFR_TMP_LIMBS_ALLOC (tn
);
189 MPN_COPY (tmp
, bp
, bn
- k
);
191 MPFR_ASSERTD (k
> 0);
193 /* Transform RNDD and RNDU to Zero / Away */
194 MPFR_ASSERTD((neg
== 0) || (neg
==1));
195 if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd1
, neg
))
202 cc
= (bp
[bn
- 1] >> s1
) & 1;
203 /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
205 cc
^= mpfr_round_raw2 (bp
, bn
, neg
, rnd2
, prec
);
206 /* cc is the new value of bit s1 in bp[bn-1] */
207 /* now round b + 2^(MPFR_EXP(b)-err) */
208 cc2
= mpn_add_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
211 /* Round to nearest */
212 /* first round b+2^(MPFR_EXP(b)-err) */
213 cc
= mpn_add_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
214 cc
= (tmp
[bn
- 1] >> s1
) & 1; /* gives 0 when cc=1 */
215 cc
^= mpfr_round_raw2 (tmp
, bn
, neg
, rnd2
, prec
);
216 /* now round b-2^(MPFR_EXP(b)-err) */
217 cc2
= mpn_sub_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
221 cc
= (bp
[bn
- 1] >> s1
) & 1;
222 cc
^= mpfr_round_raw2 (bp
, bn
, neg
, rnd2
, prec
);
223 /* now round b +/- 2^(MPFR_EXP(b)-err) */
224 cc2
= mpn_sub_1 (tmp
+ bn
- k
, bp
+ bn
- k
, k
, MPFR_LIMB_ONE
<< s
);
228 /* if cc2 is 1, then a carry or borrow propagates to the next limb */
231 MPFR_TMP_FREE(marker
);
235 cc2
= (tmp
[bn
- 1] >> s1
) & 1;
236 cc2
^= mpfr_round_raw2 (tmp
, bn
, neg
, rnd2
, prec
);
238 MPFR_TMP_FREE(marker
);