diff options
Diffstat (limited to 'src/strtofr.c')
| -rw-r--r-- | src/strtofr.c | 851 |
1 files changed, 851 insertions, 0 deletions
diff --git a/src/strtofr.c b/src/strtofr.c new file mode 100644 index 0000000..5248c10 --- /dev/null +++ b/src/strtofr.c @@ -0,0 +1,851 @@ +/* mpfr_strtofr -- set a floating-point number from a string + +Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. +Contributed by the AriC and Caramel projects, INRIA. + +This file is part of the GNU MPFR Library. + +The GNU MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +The GNU MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see +http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., +51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ + +#include <stdlib.h> /* For strtol */ +#include <ctype.h> /* For isspace */ + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + +#define MPFR_MAX_BASE 62 + +struct parsed_string { + 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). This points inside the area + allocated for the mantissa field. */ + size_t prec; /* length of mant (zero for +/-0) */ + size_t alloc; /* allocation size of mantissa */ + mpfr_exp_t exp_base; /* number of digits before the point */ + mpfr_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. + 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}, + {1UL, 2UL}, + {4004UL, 9297UL}, + {53UL, 137UL}, + {2393UL, 6718UL}, + {1UL, 3UL}, + {665UL, 2108UL}, + {4004UL, 13301UL}, + {949UL, 3283UL}, + {53UL, 190UL}, + {5231UL, 19357UL}, + {2393UL, 9111UL}, + {247UL, 965UL}, + {1UL, 4UL}, + {4036UL, 16497UL}, + {665UL, 2773UL}, + {5187UL, 22034UL}, + {4004UL, 17305UL}, + {51UL, 224UL}, + {949UL, 4232UL}, + {3077UL, 13919UL}, + {53UL, 243UL}, + {73UL, 339UL}, + {5231UL, 24588UL}, + {665UL, 3162UL}, + {2393UL, 11504UL}, + {4943UL, 24013UL}, + {247UL, 1212UL}, + {3515UL, 17414UL}, + {1UL, 5UL}, + {4415UL, 22271UL}, + {4036UL, 20533UL}, + {263UL, 1349UL}, + {665UL, 3438UL}, + {1079UL, 5621UL}, + {5187UL, 27221UL}, + {2288UL, 12093UL}, + {4004UL, 21309UL}, + {179UL, 959UL}, + {51UL, 275UL}, + {495UL, 2686UL}, + {949UL, 5181UL}, + {3621UL, 19886UL}, + {3077UL, 16996UL}, + {229UL, 1272UL}, + {53UL, 296UL}, + {109UL, 612UL}, + {73UL, 412UL}, + {1505UL, 8537UL}, + {5231UL, 29819UL}, + {283UL, 1621UL}, + {665UL, 3827UL}, + {32UL, 185UL}, + {2393UL, 13897UL}, + {1879UL, 10960UL}, + {4943UL, 28956UL}, + {409UL, 2406UL}, + {247UL, 1459UL}, + {231UL, 1370UL}, + {3515UL, 20929UL} }; +#if 0 +#define N 8 +int main () +{ + unsigned long tab[N]; + int i, n, base; + mpfr_t x, y; + mpq_t q1, q2; + int overflow = 0, base_overflow; + + mpfr_init2 (x, 200); + mpfr_init2 (y, 200); + mpq_init (q1); + mpq_init (q2); + + for (base = 2 ; base < 63 ; base ++) + { + mpfr_set_ui (x, base, MPFR_RNDN); + mpfr_log2 (x, x, MPFR_RNDN); + mpfr_ui_div (x, 1, x, MPFR_RNDN); + printf ("Base: %d x=%e ", base, mpfr_get_d1 (x)); + for (i = 0 ; i < N ; i++) + { + mpfr_floor (y, x); + tab[i] = mpfr_get_ui (y, MPFR_RNDN); + mpfr_sub (x, x, y, MPFR_RNDN); + mpfr_ui_div (x, 1, x, MPFR_RNDN); + } + for (i = N-1 ; i >= 0 ; i--) + if (tab[i] != 0) + break; + mpq_set_ui (q1, tab[i], 1); + for (i = i-1 ; i >= 0 ; i--) + { + mpq_inv (q1, q1); + mpq_set_ui (q2, tab[i], 1); + mpq_add (q1, q1, q2); + } + printf("Approx: ", base); + mpq_out_str (stdout, 10, q1); + printf (" = %e\n", mpq_get_d (q1) ); + fprintf (stderr, "{"); + mpz_out_str (stderr, 10, mpq_numref (q1)); + fprintf (stderr, "UL, "); + mpz_out_str (stderr, 10, mpq_denref (q1)); + fprintf (stderr, "UL},\n"); + if (mpz_cmp_ui (mpq_numref (q1), 1<<16-1) >= 0 + || mpz_cmp_ui (mpq_denref (q1), 1<<16-1) >= 0) + overflow = 1, base_overflow = base; + } + + mpq_clear (q2); + mpq_clear (q1); + mpfr_clear (y); + mpfr_clear (x); + if (overflow ) + printf ("OVERFLOW for base =%d!\n", base_overflow); +} +#endif + + +/* Compatible with any locale, but one still assumes that 'a', 'b', 'c', + ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like + in any ASCII-based character set). */ +static int +digit_value_in_base (int c, int base) +{ + int digit; + + MPFR_ASSERTD (base > 0 && base <= MPFR_MAX_BASE); + + if (c >= '0' && c <= '9') + digit = c - '0'; + else if (c >= 'a' && c <= 'z') + digit = (base >= 37) ? c - 'a' + 36 : c - 'a' + 10; + else if (c >= 'A' && c <= 'Z') + digit = c - 'A' + 10; + else + return -1; + + return MPFR_LIKELY (digit < base) ? digit : -1; +} + +/* Compatible with any locale, but one still assumes that 'a', 'b', 'c', + ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like + in any ASCII-based character set). */ +/* TODO: support EBCDIC. */ +static int +fast_casecmp (const char *s1, const char *s2) +{ + unsigned char c1, c2; + + do + { + c2 = *(const unsigned char *) s2++; + if (c2 == '\0') + return 0; + c1 = *(const unsigned char *) s1++; + if (c1 >= 'A' && c1 <= 'Z') + c1 = c1 - 'A' + 'a'; + } + while (c1 == c2); + return 1; +} + +/* Parse a string and fill pstr. + Return the advanced ptr too. + It returns: + -1 if invalid string, + 0 if special string (like nan), + 1 if the string is ok. + 2 if overflows + So it doesn't return the ternary value + BUT if it returns 0 (NAN or INF), the ternary value is also '0' + (ie NAN and INF are exact) */ +static int +parse_string (mpfr_t x, struct parsed_string *pstr, + const char **string, int base) +{ + const char *str = *string; + unsigned char *mant; + int point; + int res = -1; /* Invalid input return value */ + const char *prefix_str; + int decimal_point; + + decimal_point = (unsigned char) MPFR_DECIMAL_POINT; + + /* Init variable */ + pstr->mantissa = NULL; + + /* Optional leading whitespace */ + while (isspace((unsigned char) *str)) str++; + + /* An optional sign `+' or `-' */ + pstr->negative = (*str == '-'); + if (*str == '-' || *str == '+') + str++; + + /* Can be case-insensitive NAN */ + if (fast_casecmp (str, "@nan@") == 0) + { + str += 5; + goto set_nan; + } + if (base <= 16 && fast_casecmp (str, "nan") == 0) + { + str += 3; + set_nan: + /* Check for "(dummychars)" */ + if (*str == '(') + { + const char *s; + for (s = str+1 ; *s != ')' ; s++) + if (!(*s >= 'A' && *s <= 'Z') + && !(*s >= 'a' && *s <= 'z') + && !(*s >= '0' && *s <= '9') + && *s != '_') + break; + if (*s == ')') + str = s+1; + } + *string = str; + MPFR_SET_NAN(x); + /* MPFR_RET_NAN not used as the return value isn't a ternary value */ + __gmpfr_flags |= MPFR_FLAGS_NAN; + return 0; + } + + /* Can be case-insensitive INF */ + if (fast_casecmp (str, "@inf@") == 0) + { + str += 5; + goto set_inf; + } + if (base <= 16 && fast_casecmp (str, "infinity") == 0) + { + str += 8; + goto set_inf; + } + if (base <= 16 && fast_casecmp (str, "inf") == 0) + { + str += 3; + set_inf: + *string = str; + MPFR_SET_INF (x); + (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x); + return 0; + } + + /* If base=0 or 16, it may include '0x' prefix */ + prefix_str = NULL; + if ((base == 0 || base == 16) && str[0]=='0' + && (str[1]=='x' || str[1] == 'X')) + { + prefix_str = str; + base = 16; + str += 2; + } + /* If base=0 or 2, it may include '0b' prefix */ + if ((base == 0 || base == 2) && str[0]=='0' + && (str[1]=='b' || str[1] == 'B')) + { + prefix_str = str; + base = 2; + str += 2; + } + /* Else if base=0, we assume decimal base */ + if (base == 0) + base = 10; + pstr->base = base; + + /* Alloc mantissa */ + pstr->alloc = (size_t) strlen (str) + 1; + pstr->mantissa = (unsigned char*) (*__gmp_allocate_func) (pstr->alloc); + + /* Read mantissa digits */ + parse_begin: + mant = pstr->mantissa; + point = 0; + pstr->exp_base = 0; + pstr->exp_bin = 0; + + for (;;) /* Loop until an invalid character is read */ + { + int c = (unsigned char) *str++; + /* The cast to unsigned char is needed because of digit_value_in_base; + decimal_point uses this convention too. */ + if (c == '.' || c == decimal_point) + { + if (MPFR_UNLIKELY(point)) /* Second '.': stop parsing */ + break; + point = 1; + continue; + } + c = digit_value_in_base (c, base); + if (c == -1) + break; + MPFR_ASSERTN (c >= 0); /* c is representable in an unsigned char */ + *mant++ = (unsigned char) c; + if (!point) + pstr->exp_base ++; + } + str--; /* The last read character was invalid */ + + /* Update the # of char in the mantissa */ + pstr->prec = mant - pstr->mantissa; + /* Check if there are no characters in the mantissa (Invalid argument) */ + if (pstr->prec == 0) + { + /* Check if there was a prefix (in such a case, we have to read + again the mantissa without skipping the prefix) + The allocated mantissa is still big enough since we will + read only 0, and we alloc one more char than needed. + FIXME: Not really friendly. Maybe cleaner code? */ + if (prefix_str != NULL) + { + str = prefix_str; + prefix_str = NULL; + goto parse_begin; + } + goto end; + } + + /* Valid entry */ + res = 1; + MPFR_ASSERTD (pstr->exp_base >= 0); + + /* an optional exponent (e or E, p or P, @) */ + if ( (*str == '@' || (base <= 10 && (*str == 'e' || *str == 'E'))) + && (!isspace((unsigned char) str[1])) ) + { + char *endptr; + /* the exponent digits are kept in ASCII */ + mpfr_exp_t sum; + long read_exp = strtol (str + 1, &endptr, 10); + if (endptr != str+1) + str = endptr; + sum = + read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) : + read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) : + (mpfr_exp_t) read_exp; + MPFR_SADD_OVERFLOW (sum, sum, pstr->exp_base, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + res = 2, res = 3); + /* Since exp_base was positive, read_exp + exp_base can't + do a negative overflow. */ + MPFR_ASSERTD (res != 3); + pstr->exp_base = sum; + } + else if ((base == 2 || base == 16) + && (*str == 'p' || *str == 'P') + && (!isspace((unsigned char) str[1]))) + { + char *endptr; + long read_exp = strtol (str + 1, &endptr, 10); + if (endptr != str+1) + str = endptr; + pstr->exp_bin = + read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) : + read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) : + (mpfr_exp_t) read_exp; + } + + /* Remove 0's at the beginning and end of mantissa[0..prec-1] */ + mant = pstr->mantissa; + for ( ; (pstr->prec > 0) && (*mant == 0) ; mant++, pstr->prec--) + pstr->exp_base--; + for ( ; (pstr->prec > 0) && (mant[pstr->prec - 1] == 0); pstr->prec--); + pstr->mant = mant; + + /* Check if x = 0 */ + if (pstr->prec == 0) + { + MPFR_SET_ZERO (x); + if (pstr->negative) + MPFR_SET_NEG(x); + else + MPFR_SET_POS(x); + res = 0; + } + + *string = str; + end: + if (pstr->mantissa != NULL && res != 1) + (*__gmp_free_func) (pstr->mantissa, pstr->alloc); + return res; +} + +/* Transform a parsed string to a mpfr_t according to the rounding mode + and the precision of x. + Returns the ternary value. */ +static int +parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd) +{ + mpfr_prec_t prec; + mpfr_exp_t exp; + mpfr_exp_t ysize_bits; + mp_limb_t *y, *result; + int count, exact; + size_t pstr_size; + mp_size_t ysize, real_ysize; + int res, err; + MPFR_ZIV_DECL (loop); + MPFR_TMP_DECL (marker); + + /* initialize the working precision */ + prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)); + + /* compute the value y of the leading characters as long as rounding is not + possible */ + MPFR_TMP_MARK(marker); + MPFR_ZIV_INIT (loop, prec); + for (;;) + { + /* 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 = MPFR_PREC2LIMBS (prec); + /* prec bits corresponds to ysize limbs */ + ysize_bits = ysize * GMP_NUMB_BITS; + /* and to ysize_bits >= prec > MPFR_PREC (x) bits */ + y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1); + 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^(GMP_NUMB_BITS))^ysize + (in the worst case, the first digit is one and all others are zero). + i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base) + Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 => + ysize*GMP_NUMB_BITS can not overflow. + 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: + 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 == (mpfr_exp_t) pstr_size); + + /* 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); + + /* normalize y: warning we can even get ysize+1 limbs! */ + MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ + count_leading_zeros (count, y[real_ysize - 1]); + /* 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 + ysize - real_ysize, y, real_ysize, count); + if (real_ysize != 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) * GMP_NUMB_BITS + count); + } + 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 (GMP_NUMB_BITS - count) bits + to the right. FIXME: can we prove that count cannot be zero here, + since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */ + MPFR_ASSERTD (count != 0); + exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) == + MPFR_LIMB_ZERO; + /* for each bit shift increase exponent of y */ + exp = GMP_NUMB_BITS - count; + } + + /* compute base^(exp_base - pstr_size) on n limbs */ + if (IS_POW2 (pstr->base)) + { + /* Base: 2, 4, 8, 16, 32 */ + int pow2; + mpfr_exp_t tmp; + + count_leading_zeros (pow2, (mp_limb_t) pstr->base); + pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */ + MPFR_ASSERTD (0 < pow2 && pow2 <= 5); + /* exp += pow2 * (pstr->exp_base - pstr_size) + pstr->exp_bin + with overflow checking + and check that we can add/substract 2 to exp without overflow */ + MPFR_SADD_OVERFLOW (tmp, pstr->exp_base, -(mpfr_exp_t) pstr_size, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto overflow, goto underflow); + /* On some FreeBsd/Alpha, LONG_MIN/1 produced an exception + so we used to check for this before doing the division. + Since this bug is closed now (Nov 26, 2009), we remove + that check (http://www.freebsd.org/cgi/query-pr.cgi?pr=72024) */ + if (tmp > 0 && MPFR_EXP_MAX / pow2 <= tmp) + goto overflow; + else if (tmp < 0 && MPFR_EXP_MIN / pow2 >= tmp) + goto underflow; + tmp *= pow2; + MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto overflow, goto underflow); + MPFR_SADD_OVERFLOW (exp, exp, tmp, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, + goto overflow, goto underflow); + result = y; + err = 0; + } + /* case non-power-of-two-base, and pstr->exp_base > pstr_size */ + else if (pstr->exp_base > (mpfr_exp_t) pstr_size) + { + mp_limb_t *z; + mpfr_exp_t exp_z; + + result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1); + + /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */ + z = y - ysize; + /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */ + err = mpfr_mpn_exp (z, &exp_z, pstr->base, + pstr->exp_base - pstr_size, ysize); + if (err == -2) + 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 toward zero with respect + to that value. */ + + /* multiply(y = 0.mant[0]...mant[pr-1])_base by base^(exp-g): + since both y and z are rounded toward zero, so is "result" */ + mpn_mul_n (result, y, z, ysize); + + /* compute the error on the product */ + if (err == -1) + err = 0; + err ++; + + /* compute the exponent of y */ + /* exp += exp_z + ysize_bits with overflow checking + and check that we can add/substract 2 to exp without overflow */ + MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto overflow, goto underflow); + MPFR_SADD_OVERFLOW (exp, exp, exp_z, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN+2, MPFR_EXP_MAX-2, + goto overflow, goto underflow); + + /* normalize result */ + if (MPFR_LIMB_MSB (result[2 * ysize - 1]) == 0) + { + mp_limb_t *r = result + ysize - 1; + mpn_lshift (r, r, ysize + 1, 1); + /* Overflow checking not needed */ + 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; + } + /* case exp_base < pstr_size */ + else if (pstr->exp_base < (mpfr_exp_t) pstr_size) + { + mp_limb_t *z; + mpfr_exp_t exp_z; + + result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1); + + /* set y to y * K^ysize */ + y = y - ysize; /* we have allocated ysize limbs at y - ysize */ + MPN_ZERO (y, ysize); + + /* pstr_size - pstr->exp_base can overflow */ + MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto underflow, goto overflow); + + /* (z, exp_z) = base^(exp_base-pstr_size) */ + z = result + 2*ysize + 1; + err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize); + /* Since we want y/z rounded toward zero, we must get an upper + bound of z. If err >= 0, the error on z is bounded by 2^err. */ + if (err >= 0) + { + mp_limb_t cy; + unsigned long h = err / GMP_NUMB_BITS; + unsigned long l = err - h * GMP_NUMB_BITS; + + if (h >= ysize) /* not enough precision in z */ + goto next_loop; + cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l); + if (cy != 0) /* the code below requires z on ysize limbs */ + goto next_loop; + } + exact = exact && (err == -1); + if (err == -2) + goto underflow; /* FIXME: Sure? */ + if (err == -1) + err = 0; + + /* 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); + + /* exp -= exp_z + ysize_bits with overflow checking + and check that we can add/substract 2 to exp without overflow */ + MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto underflow, goto overflow); + MPFR_SADD_OVERFLOW (exp, exp, -exp_z, + mpfr_exp_t, mpfr_uexp_t, + 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 */ + if (result[2 * ysize] == MPFR_LIMB_ONE) + { + mp_limb_t *r = result + ysize; + + exact = exact && ((*r & MPFR_LIMB_ONE) == 0); + mpn_rshift (r, r, ysize + 1, 1); + /* Overflow Checking not needed */ + exp ++; + } + result += ysize; + } + /* case exp_base = pstr_size: no multiplication or division needed */ + else + { + /* base^(exp-pr) = 1 nothing to compute */ + result = y; + err = 0; + } + + /* 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); + + /* at this point, result is an approximation rounded toward 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, + (pstr->negative) ? -1 : 1, + ysize_bits - err - 1, + MPFR_RNDN, rnd, MPFR_PREC(x))) + break; + + next_loop: + /* update the prec for next loop */ + MPFR_ZIV_NEXT (loop, prec); + } /* loop */ + MPFR_ZIV_FREE (loop); + + /* round y */ + if (mpfr_round_raw (MPFR_MANT (x), result, + ysize_bits, + pstr->negative, MPFR_PREC(x), rnd, &res )) + { + /* overflow when rounding y */ + MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT; + /* Overflow Checking not needed */ + 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); + + /* DO NOT USE MPFR_SET_EXP. The exp may be out of range! */ + MPFR_SADD_OVERFLOW (exp, exp, ysize_bits, + mpfr_exp_t, mpfr_uexp_t, + MPFR_EXP_MIN, MPFR_EXP_MAX, + goto overflow, goto underflow); + MPFR_EXP (x) = exp; + res = mpfr_check_range (x, res, rnd); + goto end; + + underflow: + /* This is called when there is a huge overflow + (Real expo < MPFR_EXP_MIN << __gmpfr_emin */ + if (rnd == MPFR_RNDN) + rnd = MPFR_RNDZ; + res = mpfr_underflow (x, rnd, (pstr->negative) ? -1 : 1); + goto end; + + overflow: + res = mpfr_overflow (x, rnd, (pstr->negative) ? -1 : 1); + + end: + MPFR_TMP_FREE (marker); + return res; +} + +static void +free_parsed_string (struct parsed_string *pstr) +{ + (*__gmp_free_func) (pstr->mantissa, pstr->alloc); +} + +int +mpfr_strtofr (mpfr_t x, const char *string, char **end, int base, + mpfr_rnd_t rnd) +{ + int res; + struct parsed_string pstr; + + /* For base <= 36, parsing is case-insensitive. */ + MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 62)); + + /* If an error occured, it must return 0 */ + MPFR_SET_ZERO (x); + MPFR_SET_POS (x); + + MPFR_ASSERTN (MPFR_MAX_BASE >= 62); + res = parse_string (x, &pstr, &string, base); + /* If res == 0, then it was exact (NAN or INF), + so it is also the ternary value */ + if (MPFR_UNLIKELY (res == -1)) /* invalid data */ + res = 0; /* x is set to 0, which is exact, thus ternary value is 0 */ + else if (res == 1) + { + res = parsed_string_to_mpfr (x, &pstr, rnd); + free_parsed_string (&pstr); + } + else if (res == 2) + res = mpfr_overflow (x, rnd, (pstr.negative) ? -1 : 1); + MPFR_ASSERTD (res != 3); +#if 0 + else if (res == 3) + { + /* This is called when there is a huge overflow + (Real expo < MPFR_EXP_MIN << __gmpfr_emin */ + if (rnd == MPFR_RNDN) + rnd = MPFR_RNDZ; + res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1); + } +#endif + + if (end != NULL) + *end = (char *) string; + return res; +} |