--- /dev/null
+/* mpfr_strtofr -- set a floating-point number from a string
+
+Copyright 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include <string.h> /* For strlen */
+#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) */
+ 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.
+ 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, GMP_RNDN);
+ mpfr_log2 (x, x, GMP_RNDN);
+ mpfr_ui_div (x, 1, x, GMP_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, GMP_RNDN);
+ mpfr_sub (x, x, y, GMP_RNDN);
+ mpfr_ui_div (x, 1, x, GMP_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) * sizeof(char) + 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[1];
+ /* the exponent digits are kept in ASCII */
+ mp_exp_t read_exp = strtol (str + 1, endptr, 10);
+ mp_exp_t sum = 0;
+ if (endptr[0] != str+1)
+ str = endptr[0];
+ MPFR_ASSERTN (read_exp == (long) read_exp);
+ MPFR_SADD_OVERFLOW (sum, read_exp, pstr->exp_base,
+ mp_exp_t, mp_exp_unsigned_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[1];
+ pstr->exp_bin = (mp_exp_t) strtol (str + 1, endptr, 10);
+ if (endptr[0] != str+1)
+ str = endptr[0];
+ }
+
+ /* Remove 0's at the beginning and end of mant_s[0..prec_s-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, mp_rnd_t rnd)
+{
+ mp_prec_t prec;
+ mp_exp_t exp;
+ mp_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 y 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 = (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; /* 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.
+ 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 == (mp_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 get 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) * BITS_PER_MP_LIMB + 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 (BITS_PER_MP_LIMB - count) bits
+ 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;
+ }
+
+ /* compute base^(exp_s-pr) on n limbs */
+ if (IS_POW2 (pstr->base))
+ {
+ /* Base: 2, 4, 8, 16, 32 */
+ int pow2;
+ mp_exp_t tmp;
+
+ count_leading_zeros (pow2, (mp_limb_t) pstr->base);
+ pow2 = BITS_PER_MP_LIMB - 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, -(mp_exp_t) pstr_size,
+ mp_exp_t, mp_exp_unsigned_t,
+ MPFR_EXP_MIN, MPFR_EXP_MAX,
+ goto overflow, goto underflow);
+ /* On some FreeBsd/Alpha, LONG_MIN/1 produces an exception
+ so we check for this before doing the division */
+ if (tmp > 0 && pow2 != 1 && MPFR_EXP_MAX/pow2 <= tmp)
+ goto overflow;
+ else if (tmp < 0 && pow2 != 1 && MPFR_EXP_MIN/pow2 >= tmp)
+ goto underflow;
+ tmp *= pow2;
+ MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin,
+ mp_exp_t, mp_exp_unsigned_t,
+ MPFR_EXP_MIN, MPFR_EXP_MAX,
+ goto overflow, goto underflow);
+ MPFR_SADD_OVERFLOW (exp, exp, tmp,
+ mp_exp_t, mp_exp_unsigned_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 > (mp_exp_t) pstr_size)
+ {
+ mp_limb_t *z;
+ mp_exp_t exp_z;
+
+ result = (mp_limb_t*) MPFR_TMP_ALLOC ((2*ysize+1)*BYTES_PER_MP_LIMB);
+
+ /* 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 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);
+
+ /* 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,
+ mp_exp_t, mp_exp_unsigned_t,
+ MPFR_EXP_MIN, MPFR_EXP_MAX,
+ goto overflow, goto underflow);
+ MPFR_SADD_OVERFLOW (exp, exp, exp_z,
+ mp_exp_t, mp_exp_unsigned_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 < (mp_exp_t) pstr_size)
+ {
+ mp_limb_t *z;
+ mp_exp_t exp_z;
+
+ 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 */
+ MPN_ZERO (y, ysize);
+
+ /* pstr_size - pstr->exp_base can overflow */
+ MPFR_SADD_OVERFLOW (exp_z, (mp_exp_t) pstr_size, -pstr->exp_base,
+ mp_exp_t, mp_exp_unsigned_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);
+ 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,
+ mp_exp_t, mp_exp_unsigned_t,
+ MPFR_EXP_MIN, MPFR_EXP_MAX,
+ goto underflow, goto overflow);
+ MPFR_SADD_OVERFLOW (exp, exp, -exp_z,
+ mp_exp_t, mp_exp_unsigned_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_s-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 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,
+ (pstr->negative) ? -1 : 1,
+ ysize_bits - err - 1,
+ GMP_RNDN, rnd, MPFR_PREC(x)))
+ 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, 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,
+ mp_exp_t, mp_exp_unsigned_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 == GMP_RNDN)
+ rnd = GMP_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,
+ mp_rnd_t rnd)
+{
+ int res;
+ struct parsed_string pstr;
+
+ MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 36));
+
+ /* If an error occured, it must return 0 */
+ MPFR_SET_ZERO (x);
+ MPFR_SET_POS (x);
+
+ /* Though bases up to MPFR_MAX_BASE are supported, we require a lower
+ limit: 36. For such values <= 36, parsing is case-insensitive. */
+ MPFR_ASSERTN (MPFR_MAX_BASE >= 36);
+ 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 == GMP_RNDN)
+ rnd = GMP_RNDZ;
+ res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1);
+ }
+#endif
+
+ if (end != NULL)
+ *end = (char *) string;
+ return res;
+}
projects / msp430-gcc.git / blobdiff