diff options
Diffstat (limited to 'strtofr.c')
-rw-r--r-- | strtofr.c | 154 |
1 files changed, 104 insertions, 50 deletions
@@ -31,14 +31,24 @@ MA 02110-1301, USA. */ #define MPFR_MAX_BASE 62 struct parsed_string { - int negative; - int base; - unsigned char *mantissa, *mant; - size_t prec, alloc; - mp_exp_t exp_base, exp_bin; + int negative; /* non-zero iff the number is negative */ + int base; /* base of the string */ + unsigned char *mantissa; /* raw significand (without any point) */ + unsigned char *mant; /* stripped significand (without starting and + ending zeroes) */ + size_t prec; /* length of mant (zero for +/-0) */ + size_t alloc; /* allocation size of mantissa */ + mp_exp_t exp_base; /* number of digits before the point */ + mp_exp_t exp_bin; /* exponent in case base=2 or 16, and the pxxx + format is used (i.e., exponent is given in + base 10) */ }; -/* This table has been generated by the following program */ +/* This table has been generated by the following program. + For 2 <= b <= MPFR_MAX_BASE, + RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1] + is an upper approximation of log(2)/log(b). +*/ static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = { {1UL, 1UL}, {53UL, 84UL}, @@ -441,7 +451,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_ZIV_DECL (loop); MPFR_TMP_DECL (marker); - /* determine the minimal precision for the computation */ + /* initialize the working precision */ prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); /* compute y as long as rounding is not possible */ @@ -449,60 +459,88 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_ZIV_INIT (loop, prec); for (;;) { - /* initialize y to the value of 0.mant_s[0]...mant_s[pr-1] */ + /* Set y to the value of the ~prec most significant bits of pstr->mant + (as long as we guarantee correct rounding, we don't need to get + exactly prec bits). */ ysize = (prec - 1) / BITS_PER_MP_LIMB + 1; + /* prec bits corresponds to ysize limbs */ ysize_bits = ysize * BITS_PER_MP_LIMB; + /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ y = (mp_limb_t*) MPFR_TMP_ALLOC ((2 * ysize + 1) * sizeof (mp_limb_t)); - y += ysize; - - /* required precision for str to have ~ysize limbs - We must have (2^(BITS_PER_MP_LIMB))^ysize ~= base^pstr_size - aka pstr_size = ceil (ysize*BITS_PER_MP_LIMB/log2(base)) - ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => + y += ysize; /* y has (ysize+1) allocated limbs */ + + /* pstr_size is the number of characters we read in pstr->mant + to have at least ysize full limbs. + We must have base^(pstr_size-1) >= (2^(BITS_PER_MP_LIMB))^ysize + (in the worst case, the first digit is one and all others are zero). + i.e., pstr_size >= 1 + ysize*BITS_PER_MP_LIMB/log2(base) + Since ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 => ysize*BITS_PER_MP_LIMB can not overflow. - Compute pstr_size = ysize_bits * Num / Den - Where Num/Den ~ 1/log2(base) + We compute pstr_size = 1 + ceil(ysize_bits * Num / Den) + where Num/Den >= 1/log2(base) It is not exactly ceil(1/log2(base)) but could be one more (base 2) - Quite ugly since it tries to avoid overflow */ - pstr_size = ( ysize_bits / RedInvLog2Table[pstr->base-2][1] - * RedInvLog2Table[pstr->base-2][0] ) - + (( ysize_bits % RedInvLog2Table[pstr->base-2][1]) - * RedInvLog2Table[pstr->base-2][0] - / RedInvLog2Table[pstr->base-2][1] ) - + 1; - + Quite ugly since it tries to avoid overflow: + let Num = RedInvLog2Table[pstr->base-2][0] + and Den = RedInvLog2Table[pstr->base-2][1], + and ysize_bits = a*Den+b, + then ysize_bits * Num/Den = a*Num + (b * Num)/Den, + thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den + */ + { + unsigned long Num = RedInvLog2Table[pstr->base-2][0]; + unsigned long Den = RedInvLog2Table[pstr->base-2][1]; + pstr_size = ((ysize_bits / Den) * Num) + + (((ysize_bits % Den) * Num + Den - 1) / Den) + + 1; + } + + /* since pstr_size corresponds to at least ysize_bits full bits, + and ysize_bits > prec, the weight of the neglected part of + pstr->mant (if any) is < ulp(y) < ulp(x) */ + + /* if the number of wanted characters is more than what we have in + pstr->mant, round it down */ if (pstr_size >= pstr->prec) pstr_size = pstr->prec; MPFR_ASSERTD (pstr_size == (mp_exp_t) pstr_size); - /* convert str into binary */ + /* convert str into binary: note that pstr->mant is big endian, + thus no offset is needed */ real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base); - MPFR_ASSERTD ( real_ysize <= ysize+1); + MPFR_ASSERTD (real_ysize <= ysize+1); /* normalize y: warning we can get even get ysize+1 limbs! */ - MPFR_ASSERTD (y[real_ysize - 1] != 0); + MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ count_leading_zeros (count, y[real_ysize - 1]); - exact = (real_ysize <= ysize); - if (exact != 0) /* shift y to the left in that case y should be exact */ + /* exact means that the number of limbs of the output of mpn_set_str + is less or equal to ysize */ + exact = real_ysize <= ysize; + if (exact) /* shift y to the left in that case y should be exact */ { + /* we have enough limbs to store {y, real_ysize} */ /* shift {y, num_limb} for count bits to the left */ if (count != 0) - mpn_lshift (y, y, real_ysize, count); - /* shift {y, num_limb} for (ysize-num_limb) limbs to the left */ + mpn_lshift (y + ysize - real_ysize, y, real_ysize, count); if (real_ysize != ysize) { - MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); + if (count == 0) + MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize); MPN_ZERO (y, ysize - real_ysize); } /* for each bit shift decrease exponent of y */ /* (This should not overflow) */ exp = - ((ysize - real_ysize) * BITS_PER_MP_LIMB + count); } - else /* shift y for the right */ + else /* shift y to the right, by doing this we might lose some + bits from the result of mpn_set_str (in addition to the + characters neglected from pstr->mant) */ { /* shift {y, num_limb} for (BITS_PER_MP_LIMB - count) bits - to the right */ - mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count); + to the right. FIXME: can we prove that count cannot be zero here, + since mpn_rshift does not accept a shift of BITS_PER_MP_LIMB? */ + MPFR_ASSERTD (count != 0); + exact = mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count) == + MPFR_LIMB_ZERO; /* for each bit shift increase exponent of y */ exp = BITS_PER_MP_LIMB - count; } @@ -542,7 +580,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) result = y; err = 0; } - /* case pstr->exp_base > pstr_size */ + /* case non-power-of-two-base, and pstr->exp_base > pstr_size */ else if (pstr->exp_base > (mp_exp_t) pstr_size) { mp_limb_t *z; @@ -559,6 +597,13 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) goto overflow; exact = exact && (err == -1); + /* If exact is non zero, then z equals exactly the value of the + pstr_size most significant digits from pstr->mant, i.e., the + only difference can come from the neglected pstr->prec-pstr_size + least significant digits of pstr->mant. + If exact is zero, then z is rounded towards zero with respect + to that value. */ + /* multiply(y = 0.mant_s[0]...mant_s[pr-1])_base by base^(exp_s-g) */ mpn_mul_n (result, y, z, ysize); @@ -588,6 +633,8 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) exp --; } + /* if the low ysize limbs of {result, 2*ysize} are all zero, + then the result is still "exact" (if it was before) */ exact = exact && (mpn_scan1 (result, 0) >= (unsigned long) ysize_bits); result += ysize; @@ -598,7 +645,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) mp_limb_t *z; mp_exp_t exp_z; - result = (mp_limb_t*) MPFR_TMP_ALLOC ( (3*ysize+1) * BYTES_PER_MP_LIMB); + result = (mp_limb_t*) MPFR_TMP_ALLOC ((3*ysize+1) * BYTES_PER_MP_LIMB); /* set y to y * K^ysize */ y = y - ysize; /* we have allocated ysize limbs at y - ysize */ @@ -622,7 +669,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) /* compute y / z */ /* result will be put into result + n, and remainder into result */ mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y, - 2*ysize, z, ysize); + 2 * ysize, z, ysize); /* exp -= exp_z + ysize_bits with overflow checking and check that we can add/substract 2 to exp without overflow */ @@ -635,6 +682,8 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, goto overflow, goto underflow); err += 2; + /* if the remainder of the division is zero, then the result is + still "exact" if it was before */ exact = exact && (mpn_popcount (result, ysize) == 0); /* normalize result */ @@ -648,7 +697,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) } result += ysize; } - /* case exp_base = pstr_size */ + /* case exp_base = pstr_size: no multiplication or division needed */ else { /* base^(exp_s-pr) = 1 nothing to compute */ @@ -656,18 +705,16 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) err = 0; } - if (pstr_size < pstr->prec && exact - && ((rnd == GMP_RNDN) - || ((rnd == GMP_RNDD && pstr->negative) - || (rnd == GMP_RNDU && !pstr->negative)))) - { - /* Some significant digits might have been forgotten, if so result - is not exact. */ - size_t i; + /* If result is exact, we still have to consider the neglected part + of the input string. For a directed rounding, in that case we could + still correctly round, since the neglected part is less than + one ulp, but that would make the code more complex, and give a + speedup for rare cases only. */ + exact = exact && (pstr_size == pstr->prec); - for (i = pstr_size; exact && i < pstr->prec; i++) - exact = pstr->mant[i] == 0; - } + /* at this point, result is an approximation rounded towards zero + of the pstr_size most significant digits of pstr->mant, with + equality in case exact is non-zero. */ /* test if rounding is possible, and if so exit the loop */ if (exact || mpfr_can_round_raw (result, ysize, @@ -692,6 +739,13 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd) exp ++; } + if (res == 0) /* fix ternary value */ + { + exact = exact && (pstr_size == pstr->prec); + if (!exact) + res = (pstr->negative) ? 1 : -1; + } + /* Set sign of x before exp since check_range needs a valid sign */ (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); |