diff options
Diffstat (limited to 'Build/source/libs/mpfr/mpfr-src/src/strtofr.c')
-rw-r--r-- | Build/source/libs/mpfr/mpfr-src/src/strtofr.c | 978 |
1 files changed, 0 insertions, 978 deletions
diff --git a/Build/source/libs/mpfr/mpfr-src/src/strtofr.c b/Build/source/libs/mpfr/mpfr-src/src/strtofr.c deleted file mode 100644 index c6aeb739133..00000000000 --- a/Build/source/libs/mpfr/mpfr-src/src/strtofr.c +++ /dev/null @@ -1,978 +0,0 @@ -/* mpfr_strtofr -- set a floating-point number from a string - -Copyright 2004-2020 Free Software Foundation, Inc. -Contributed by the AriC and Caramba 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 -https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., -51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ - -#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, + exponent - except in case of binary exponent (exp_bin) */ - mpfr_exp_t exp_bin; /* binary exponent of the pxxx format for - base = 2 or 16 */ -}; - -/* 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 to log(2)/log(b), no larger than 1. - Note: these numbers must fit on 16 bits, thus unsigned int is OK. -*/ -static const unsigned int RedInvLog2Table[MPFR_MAX_BASE-1][2] = { - {1, 1}, - {53, 84}, - {1, 2}, - {4004, 9297}, - {53, 137}, - {2393, 6718}, - {1, 3}, - {665, 2108}, - {4004, 13301}, - {949, 3283}, - {53, 190}, - {5231, 19357}, - {2393, 9111}, - {247, 965}, - {1, 4}, - {4036, 16497}, - {665, 2773}, - {5187, 22034}, - {4004, 17305}, - {51, 224}, - {949, 4232}, - {3077, 13919}, - {53, 243}, - {73, 339}, - {5231, 24588}, - {665, 3162}, - {2393, 11504}, - {4943, 24013}, - {247, 1212}, - {3515, 17414}, - {1, 5}, - {4415, 22271}, - {4036, 20533}, - {263, 1349}, - {665, 3438}, - {1079, 5621}, - {5187, 27221}, - {2288, 12093}, - {4004, 21309}, - {179, 959}, - {51, 275}, - {495, 2686}, - {949, 5181}, - {3621, 19886}, - {3077, 16996}, - {229, 1272}, - {53, 296}, - {109, 612}, - {73, 412}, - {1505, 8537}, - {5231, 29819}, - {283, 1621}, - {665, 3827}, - {32, 185}, - {2393, 13897}, - {1879, 10960}, - {4943, 28956}, - {409, 2406}, - {247, 1459}, - {231, 1370}, - {3515, 20929} }; -#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*) mpfr_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); - - /* FIXME: In the code below (both cases), if the exponent from the - string is large, it will be replaced by MPFR_EXP_MIN or MPFR_EXP_MAX, - i.e. it will have a different value. This may not change the result - in most cases, but there is no guarantee on very long strings when - mpfr_exp_t is a 32-bit type, as the exponent could be brought back - to the current exponent range. */ - - /* 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) - mpfr_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 precx, prec, ysize_bits, pstr_size; - mpfr_exp_t exp; - mp_limb_t *result; - int count, exact; - mp_size_t ysize, real_ysize, diff_ysize; - int res, err; - const int extra_limbs = GMP_NUMB_BITS >= 12 ? 1 : 2; /* see below */ - MPFR_ZIV_DECL (loop); - MPFR_TMP_DECL (marker); - - /* initialize the working precision */ - precx = MPFR_GET_PREC (x); - prec = precx + MPFR_INT_CEIL_LOG2 (precx); - - /* Compute the value y of the leading characters as long as rounding is not - possible. - Note: We have some integer overflow checking using MPFR_EXP_MIN and - MPFR_EXP_MAX in this loop. Thanks to the large margin between these - extremal values of the mpfr_exp_t type and the valid minimum/maximum - exponents, such integer overflows would correspond to real underflow - or overflow on the result (possibly except in huge precisions, which - are disregarded here; anyway, in practice, such issues could occur - only with 32-bit precision and exponent types). Such checks could be - extended to real early underflow/overflow checking, in order to avoid - useless computations in such cases; in such a case, be careful that - the approximation errors need to be taken into account. */ - MPFR_TMP_MARK(marker); - MPFR_ZIV_INIT (loop, prec); - for (;;) - { - mp_limb_t *y0, *y; - - /* y will be regarded as a number with precision prec. */ - ysize = MPFR_PREC2LIMBS (prec); - /* prec bits corresponds to ysize limbs */ - ysize_bits = (mpfr_prec_t) ysize * GMP_NUMB_BITS; - MPFR_ASSERTD (ysize_bits >= prec); - /* and to ysize_bits >= prec > precx bits. */ - /* We need to allocate one more limb as specified by mpn_set_str - (a limb may be written in rp[rn]). Note that the manual of GMP - up to 5.1.3 was incorrect on this point. - See the following discussion: - https://gmplib.org/list-archives/gmp-bugs/2013-December/003267.html */ - y0 = MPFR_TMP_LIMBS_ALLOC (2 * ysize + extra_limbs + 1); - y = y0 + ysize; /* y has (ysize + extra_limbs + 1) allocated limbs */ - - /* pstr_size is the number of bytes we want to read from pstr->mant - to fill at least ysize full limbs with mpn_set_str. - 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 1/log2(base) <= Num/Den <= 1 - 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 - - Note: denoting m = pstr_size and n = ysize_bits, assuming we have - m = 1 + ceil(n/log2(b)), i.e., b^(m-1) >= 2^n > b^(m-2), then - b^(m-1)/2^n < b, and since we consider m characters of the input, - the corresponding part is less than b^m < b^2*2^n. - This implies that if b^2 < 2^GMP_NUMB_BITS, which for b <= 62 holds - for GMP_NUMB_BITS >= 12, we have real_ysize <= ysize+1 below - (this also implies that for GMP_NUMB_BITS >= 13, the number of bits - of y[real_ysize-1] below is less than GMP_NUMB_BITS, thus - count < GMP_NUMB_BITS). - Warning: for GMP_NUMB_BITS=8, we can have real_ysize = ysize + 2! - Hence the allocation above for ysize + extra_limbs limbs. - */ - { - unsigned int Num = RedInvLog2Table[pstr->base-2][0]; - unsigned int Den = RedInvLog2Table[pstr->base-2][1]; - MPFR_ASSERTD (Num <= Den && Den <= 65535); /* thus no overflow */ - pstr_size = (ysize_bits / Den) * Num - + ((unsigned long) (ysize_bits % Den) * Num + Den - 1) / Den - + 1; - } - - /* Since pstr_size corresponds to at least ysize_bits 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 bytes is more than what is available - in pstr->mant, i.e. pstr->prec, reduce it to pstr->prec. */ - if (pstr_size > pstr->prec) - pstr_size = pstr->prec; - - /* Convert str (potentially truncated to pstr_size) 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); - - /* See above for the explanation of the following assertion. */ - MPFR_ASSERTD (real_ysize <= ysize + extra_limbs); - - /* The Boolean "exact" will attempt to track exactness of the result: - If it is true, then this means that the result is exact, allowing - termination, even though the rounding test may not succeed. - Conversely, if the result is exact, then "exact" will not - necessarily be true at the end of the Ziv loop, but we will need - to make sure that at some point, "exact" will be true in order to - guarantee termination. FIXME: check that. */ - /* First, consider the part of the input string that has been ignored. - Note that the trailing zeros have been removed in parse_string, so - that if something has been ignored, it must be non-zero. */ - exact = pstr_size == pstr->prec; - - /* Normalize y and set the initial value of its exponent exp, which - is 0 when y is not shifted. - Since pstr->mant was normalized, mpn_set_str guarantees that - the most significant limb is non-zero. */ - MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */ - count_leading_zeros (count, y[real_ysize - 1]); - diff_ysize = ysize - real_ysize; - MPFR_LOG_MSG (("diff_ysize = %ld\n", (long) diff_ysize)); - if (diff_ysize >= 0) - { - /* We have enough limbs to store {y, real_ysize} exactly - in {y, ysize}, so that we can do a left shift, without - losing any information ("exact" will not change). */ - if (count != 0) - mpn_lshift (y + diff_ysize, y, real_ysize, count); - if (diff_ysize > 0) - { - if (count == 0) - mpn_copyd (y + diff_ysize, y, real_ysize); - MPN_ZERO (y, diff_ysize); - } - /* exp = negation of the total shift count, avoiding overflows. */ - exp = - ((mpfr_exp_t) diff_ysize * GMP_NUMB_BITS + count); - } - else - { - /* Shift {y, real_ysize} for (GMP_NUMB_BITS - count) bits to the - right, and put the ysize most significant limbs into {y, ysize}. - We have either real_ysize = ysize + 1 or real_ysize = ysize + 2 - (only possible with extra_limbs == 2). */ - MPFR_ASSERTD (diff_ysize == -1 || - (extra_limbs == 2 && diff_ysize == -2)); - if (count != 0) - { - /* Before doing the shift, consider the limb that will entirely - be lost if real_ysize = ysize + 2. */ - exact = exact && (diff_ysize == -1 || y[0] == MPFR_LIMB_ZERO); - /* mpn_rshift allows overlap, provided destination <= source */ - /* FIXME: The bits lost due to mpn_rshift are not taken - into account in the error analysis below! */ - if (mpn_rshift (y, y - (diff_ysize + 1), real_ysize, - GMP_NUMB_BITS - count) != MPFR_LIMB_ZERO) - exact = 0; /* some non-zero bits have been shifted out */ - } - else - { - /* the case real_ysize = ysize + 2 with count = 0 cannot happen - even with GMP_NUMB_BITS = 8 since 62^2 < 256^2/2 */ - MPFR_ASSERTD (diff_ysize == -1); - exact = exact && y[0] == MPFR_LIMB_ZERO; - /* copy {y+real_ysize-ysize, ysize} to {y, ysize} */ - mpn_copyi (y, y + 1, real_ysize - 1); - } - /* exp = shift count */ - /* TODO: add some explanations about what exp means exactly. */ - exp = GMP_NUMB_BITS * (- diff_ysize) - 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; - - MPFR_LOG_MSG (("case 1 (base = power of 2)\n", 0)); - - 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/subtract 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; - - MPFR_LOG_MSG (("case 2 (exp_base > pstr_size)\n", 0)); - - result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1); - - /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */ - z = y0; - /* 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/subtract 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) >= 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; - - MPFR_LOG_MSG (("case 3 (exp_base < pstr_size)\n", 0)); - - result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1); - - /* y0 = y * K^ysize */ - MPN_ZERO (y0, 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^(pstr_size - exp_base) */ - z = result + 2*ysize + 1; - err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize); - - /* Now {z, ysize} * 2^(exp_z_out - ysize_bits) is an approximation - to base^exp_z_in (denoted b^e below), rounded toward zero, with: - * if err = -1, the result is exact; - * if err = -2, an overflow occurred in the computation of exp_z; - * otherwise the error is bounded by 2^err ulps. - Thus the exact value of b^e is between z and z + 2^err, where - z is {z, ysize} properly scaled by a power of 2. Then the error - will be: - y/b^e - trunc(y/z) = eps1 + eps2 - with - eps1 = y/b^e - y/z <= 0 - eps2 = y/z - trunc(y/z) >= 0 - thus the errors will (partly) compensate, giving a bound - max(|eps1|,|eps2|). - In addition, there is a 3rd error eps3 since y might be the - conversion of only a part of the character string, and/or y - might be truncated by the mpn_rshift call above: - eps3 = exact_y/b^e - y/b^e >= 0. - */ - if (err == -2) - goto underflow; /* FIXME: Sure? */ - else if (err == -1) - err = 0; /* see the note below */ - else - exact = 0; - - /* exp -= exp_z + ysize_bits with overflow checking - and check that we can add/subtract 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); - - /* Compute the integer division y/z rounded toward zero. - The quotient will be put at result + ysize (size: ysize + 1), - and the remainder at result (size: ysize). - Both the dividend {y, 2*ysize} and the divisor {z, ysize} are - normalized, i.e., the most significant bit of their most - significant limb is 1. */ - MPFR_ASSERTD (MPFR_LIMB_MSB (y0[2 * ysize - 1]) != 0); - MPFR_ASSERTD (MPFR_LIMB_MSB (z[ysize - 1]) != 0); - mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y0, - 2 * ysize, z, ysize); - - /* The truncation error of the mpn_tdiv_qr call (eps2 above) is at - most 1 ulp. Idem for the error eps3, which has the same sign, - thus eps2 + eps3 <= 2 ulps. - FIXME: For eps3, this is not obvious and should be explained. - For the error eps1 coming from the approximation to b^e, - we have (still up to a power-of-2 normalization): - y/z - y/b^e = y * (b^e-z) / (z * b^e) <= y * 2^err / (z * b^e). - We have to convert that error in terms of ulp(trunc(y/z)). - We first have ulp(trunc(y/z)) = ulp(y/z). - - FIXME: There must be some discussion about the exponents, - because up to a power of 2, 1/2 <= |y/z| < 1 and - 1 <= |y/z| < 2 are equivalent and give no information. - Moreover 1/2 <= b^e < 1 has not been explained and may - hide mistakes since one may have 1/2 <= z < 1 < b^e. - - Since both y and z are normalized, the quotient - {result+ysize, ysize+1} has exactly ysize limbs, plus maybe one - bit (this corresponds to the MPFR_ASSERTD below): - * if the quotient has exactly ysize limbs, then 1/2 <= |y/z| < 1 - (up to a power of 2) and since 1/2 <= b^e < 1, the error is at - most 2^(err+1) ulps; - * if the quotient has one extra bit, then 1 <= |y/z| < 2 - (up to a power of 2) and since 1/2 <= b^e < 1, the error is at - most 2^(err+2) ulps; but since we will shift the result right - below by one bit, the final error will be at most 2^(err+1) ulps - too. - - Thus the error is: - * at most 2^(err+1) ulps for eps1 - * at most 2 ulps for eps2 + eps3, which is of opposite sign - and we can bound the error by 2^(err+1) ulps in all cases. - - Note: If eps1 was 0, the error would be bounded by 2 ulps, - thus replacing err = -1 by err = 0 above was the right thing - to do, since 2^(0+1) = 2. - */ - MPFR_ASSERTD (result[2 * ysize] <= 1); - - err += 1; /* see above for the explanation of the +1 term */ - - /* 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 - { - MPFR_LOG_MSG (("case 4 (exp_base = pstr_size)\n", 0)); - - /* base^(exp-pr) = 1 nothing to compute */ - result = y; - err = 0; - } - - MPFR_LOG_MSG (("exact = %d, err = %d, precx = %Pu\n", - exact, err, precx)); - - /* 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. - Note: we also need to be able to determine the correct ternary value, - thus we use the precx + (rnd == MPFR_RNDN) trick. - For example if result = xxx...xxx111...111 and rnd = RNDN, - then we know the correct rounding is xxx...xx(x+1), but we cannot know - the correct ternary value. */ - if (exact || mpfr_round_p (result, ysize, ysize_bits - err - 1, - precx + (rnd == MPFR_RNDN))) - break; - - /* 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, precx, rnd, &res)) - { - /* overflow when rounding y */ - MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT; - /* Overflow Checking not needed */ - exp ++; - } - - /* Note: if exact <> 0, then the approximation {result, ysize} is exact, - thus no double-rounding can occur: - (a) either the ternary value res is non-zero, and it is the correct - ternary value that we should return - (b) or the ternary value res is zero, and we should return 0. */ - - /* 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) -{ - mpfr_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 occurred, it must return 0. */ - MPFR_SET_ZERO (x); - MPFR_SET_POS (x); - - MPFR_STAT_STATIC_ASSERT (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; -} |