From d44793e4c78cf1a88d00046d32f8ff403dea655b Mon Sep 17 00:00:00 2001 From: vlefevre Date: Mon, 19 Jan 2004 03:01:10 +0000 Subject: _mpfr_ceil -> mpfr_ceil_double + check for overflow. git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@2632 280ebfd0-de03-0410-8827-d642c229c3f4 --- get_str.c | 57 +++++++++++++++++++++++++++++++-------------------------- 1 file changed, 31 insertions(+), 26 deletions(-) (limited to 'get_str.c') diff --git a/get_str.c b/get_str.c index 229350526..13beb3672 100644 --- a/get_str.c +++ b/get_str.c @@ -24,20 +24,21 @@ MA 02111-1307, USA. */ #include #include #include +#include #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" #include "mpfr.h" #include "mpfr-impl.h" -static double _mpfr_ceil _MPFR_PROTO ((double)); +static double mpfr_ceil_double _MPFR_PROTO ((double)); static int mpfr_get_str_aux _MPFR_PROTO ((char *, mp_exp_t *, mp_limb_t *, mp_size_t, mp_exp_t, long, int, size_t, mp_rnd_t)); static mp_exp_t mpfr_get_str_compute_g _MPFR_PROTO ((int, mp_exp_t)); static char num_to_text[] = "0123456789abcdefghijklmnopqrstuvwxyz"; -/* for 2 <= b <= 36, log_b2[b-2] + log_b2_low[b-2] is a 76-bit upper +/* for 2 <= b <= 36, log_b2[b-2] + log_b2_low[b-2] is a 76-bit upper approximation of log(2)/log(b), with log_b2[b-2] having 23 significative bits only. These approximations were computed with the following program. @@ -177,14 +178,18 @@ static const double log_b2_low[35] = { } static double -_mpfr_ceil (double x) +mpfr_ceil_double (double x) { - double y = (double) (long int) x; + double y; + /* Note: this function should be rewritten to avoid the possible + overflow. */ + MPFR_ASSERTN(x >= (double) LONG_MIN && x <= (double) LONG_MAX); + y = (double) (long int) x; return ((y < x) ? y + 1.0 : y); } /* Input: an approximation r*2^f of an real Y, with |r*2^f-Y| <= 2^(e+f). - Returns if possible in the string s the mantissa corresponding to + Returns if possible in the string s the mantissa corresponding to the integer nearest to Y, within the direction rnd, and returns the the exponent in exp. n is the number of limbs of r. @@ -193,7 +198,7 @@ _mpfr_ceil (double x) b is the wanted base (2 <= b <= 36). m is the number of wanted digits in the mantissa. rnd is the rounding mode. - It is assumed that b^(m-1) <= Y < b^(m+1), thus the returned value + It is assumed that b^(m-1) <= Y < b^(m+1), thus the returned value satisfies b^(m-1) <= rnd(Y) < b^(m+1). Return value: @@ -284,18 +289,18 @@ mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n, /* if size_s1 = m + 2, necessarily we have b^(m+1) as result, and the result will not change */ - /* so we have to double-round only when size_s1 = m + 1 and - (i) the result is inexact + /* so we have to double-round only when size_s1 = m + 1 and + (i) the result is inexact (ii) or the last digit is non-zero */ if ((size_s1 == m + 1) && ((dir != 0) || (str1[size_s1 - 1] != 0))) { /* rounding mode */ - rnd1 = rnd; + rnd1 = rnd; /* round to nearest case */ - if (rnd == GMP_RNDN) + if (rnd == GMP_RNDN) { - if (2 * str1[size_s1 - 1] == b) + if (2 * str1[size_s1 - 1] == b) { if (dir == -1) rnd1 = GMP_RNDU; @@ -336,7 +341,7 @@ mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n, } str1[i]++; } - dir = 1; + dir = 1; } /* round toward zero (truncate) */ else @@ -354,7 +359,7 @@ mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n, cannot_round: dir = 2; } - + TMP_FREE(marker); return (dir); @@ -366,12 +371,12 @@ mpfr_get_str_compute_g (int beta, mp_exp_t e) { double g0, g1; mp_exp_t g; - + g0 = (double) e * log_b2[beta - 2]; g1 = (double) e * log_b2_low[beta - 2]; - g = (mp_exp_t) _mpfr_ceil (g0); + g = (mp_exp_t) mpfr_ceil_double (g0); g0 -= (double) g; - return g + (mp_exp_t) _mpfr_ceil (g0 + g1); + return g + (mp_exp_t) mpfr_ceil_double (g0 + g1); } /* prints the mantissa of x in the string s, and writes the corresponding @@ -393,7 +398,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd int exact; /* exact result */ mp_exp_t exp, g; mp_exp_t prec, log_2prec; /* precision of the computation */ - long err; + long err; mp_limb_t *a; mp_exp_t exp_a; mp_limb_t *result; @@ -425,7 +430,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd If EXP(x)=-1, one bit is lost... */ { int k, lost; - + count_leading_zeros (k, (mp_limb_t) b); k = BITS_PER_MP_LIMB - k - 1; /* b = 2^k */ lost = (-MPFR_EXP(x)) % k; @@ -433,7 +438,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd lost += k; m += lost; } - m = (size_t) _mpfr_ceil (__mp_bases[b].chars_per_bit_exactly * (double) m); + m = (size_t) mpfr_ceil_double (__mp_bases[b].chars_per_bit_exactly * (double) m); if (m < 2) m = 2; } @@ -507,7 +512,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd count_leading_zeros (pow2, (mp_limb_t) b); pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */ - + /* set MPFR_EXP(x) = f*pow2 + r, 1 <= r <= pow2 */ f = (MPFR_GET_EXP (x) - 1) / pow2; r = MPFR_GET_EXP (x) - f * pow2; @@ -526,7 +531,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd nb = n * BITS_PER_MP_LIMB - prec; /* round xp to the precision prec, and put it into x1 put the carry into x1[n] */ - if ((x1[n] = mpfr_round_raw (x1, xp, MPFR_PREC(x), + if ((x1[n] = mpfr_round_raw (x1, xp, MPFR_PREC(x), MPFR_IS_STRICTNEG(x), prec, rnd, &inexp))) { @@ -541,7 +546,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd else /* 2^prec needs still m digits, but x1 may need n+1 limbs */ n ++; } - + /* it remains to shift x1 by nb limbs to the right, since mpn_get_str expects a right-normalized number */ if (nb != 0) @@ -551,7 +556,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd if (x1[n - 1] == 0) n --; } - + mpn_get_str ((unsigned char*) s, b, x1, n); for (i=0; i n) - exact = (exact && + exact = (exact && ((mpn_scan1 (xp, 0) >= (nx - n) * BITS_PER_MP_LIMB))); /* we loose one more bit in the multiplication, except when err=0 where we loose two bits */ err = (err <= 0) ? 2 : err + 1; - + /* result = a * x */ result = (mp_limb_t*) TMP_ALLOC ((n + nx1) * sizeof (mp_limb_t)); mpn_mul (result, a, n, x1, nx1); -- cgit v1.2.1