14a238c70SJohn Marino /* mpfr_strtofr -- set a floating-point number from a string
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
54a238c70SJohn Marino
64a238c70SJohn Marino This file is part of the GNU MPFR Library.
74a238c70SJohn Marino
84a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
94a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
104a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
114a238c70SJohn Marino option) any later version.
124a238c70SJohn Marino
134a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
144a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
154a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
164a238c70SJohn Marino License for more details.
174a238c70SJohn Marino
184a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
194a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
204a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
214a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
224a238c70SJohn Marino
234a238c70SJohn Marino #include <stdlib.h> /* For strtol */
244a238c70SJohn Marino #include <ctype.h> /* For isspace */
254a238c70SJohn Marino
264a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H
274a238c70SJohn Marino #include "mpfr-impl.h"
284a238c70SJohn Marino
294a238c70SJohn Marino #define MPFR_MAX_BASE 62
304a238c70SJohn Marino
314a238c70SJohn Marino struct parsed_string {
324a238c70SJohn Marino int negative; /* non-zero iff the number is negative */
334a238c70SJohn Marino int base; /* base of the string */
344a238c70SJohn Marino unsigned char *mantissa; /* raw significand (without any point) */
354a238c70SJohn Marino unsigned char *mant; /* stripped significand (without starting and
36*ab6d115fSJohn Marino ending zeroes). This points inside the area
37*ab6d115fSJohn Marino allocated for the mantissa field. */
384a238c70SJohn Marino size_t prec; /* length of mant (zero for +/-0) */
394a238c70SJohn Marino size_t alloc; /* allocation size of mantissa */
404a238c70SJohn Marino mpfr_exp_t exp_base; /* number of digits before the point */
414a238c70SJohn Marino mpfr_exp_t exp_bin; /* exponent in case base=2 or 16, and the pxxx
424a238c70SJohn Marino format is used (i.e., exponent is given in
434a238c70SJohn Marino base 10) */
444a238c70SJohn Marino };
454a238c70SJohn Marino
464a238c70SJohn Marino /* This table has been generated by the following program.
474a238c70SJohn Marino For 2 <= b <= MPFR_MAX_BASE,
484a238c70SJohn Marino RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1]
494a238c70SJohn Marino is an upper approximation of log(2)/log(b).
504a238c70SJohn Marino */
514a238c70SJohn Marino static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = {
524a238c70SJohn Marino {1UL, 1UL},
534a238c70SJohn Marino {53UL, 84UL},
544a238c70SJohn Marino {1UL, 2UL},
554a238c70SJohn Marino {4004UL, 9297UL},
564a238c70SJohn Marino {53UL, 137UL},
574a238c70SJohn Marino {2393UL, 6718UL},
584a238c70SJohn Marino {1UL, 3UL},
594a238c70SJohn Marino {665UL, 2108UL},
604a238c70SJohn Marino {4004UL, 13301UL},
614a238c70SJohn Marino {949UL, 3283UL},
624a238c70SJohn Marino {53UL, 190UL},
634a238c70SJohn Marino {5231UL, 19357UL},
644a238c70SJohn Marino {2393UL, 9111UL},
654a238c70SJohn Marino {247UL, 965UL},
664a238c70SJohn Marino {1UL, 4UL},
674a238c70SJohn Marino {4036UL, 16497UL},
684a238c70SJohn Marino {665UL, 2773UL},
694a238c70SJohn Marino {5187UL, 22034UL},
704a238c70SJohn Marino {4004UL, 17305UL},
714a238c70SJohn Marino {51UL, 224UL},
724a238c70SJohn Marino {949UL, 4232UL},
734a238c70SJohn Marino {3077UL, 13919UL},
744a238c70SJohn Marino {53UL, 243UL},
754a238c70SJohn Marino {73UL, 339UL},
764a238c70SJohn Marino {5231UL, 24588UL},
774a238c70SJohn Marino {665UL, 3162UL},
784a238c70SJohn Marino {2393UL, 11504UL},
794a238c70SJohn Marino {4943UL, 24013UL},
804a238c70SJohn Marino {247UL, 1212UL},
814a238c70SJohn Marino {3515UL, 17414UL},
824a238c70SJohn Marino {1UL, 5UL},
834a238c70SJohn Marino {4415UL, 22271UL},
844a238c70SJohn Marino {4036UL, 20533UL},
854a238c70SJohn Marino {263UL, 1349UL},
864a238c70SJohn Marino {665UL, 3438UL},
874a238c70SJohn Marino {1079UL, 5621UL},
884a238c70SJohn Marino {5187UL, 27221UL},
894a238c70SJohn Marino {2288UL, 12093UL},
904a238c70SJohn Marino {4004UL, 21309UL},
914a238c70SJohn Marino {179UL, 959UL},
924a238c70SJohn Marino {51UL, 275UL},
934a238c70SJohn Marino {495UL, 2686UL},
944a238c70SJohn Marino {949UL, 5181UL},
954a238c70SJohn Marino {3621UL, 19886UL},
964a238c70SJohn Marino {3077UL, 16996UL},
974a238c70SJohn Marino {229UL, 1272UL},
984a238c70SJohn Marino {53UL, 296UL},
994a238c70SJohn Marino {109UL, 612UL},
1004a238c70SJohn Marino {73UL, 412UL},
1014a238c70SJohn Marino {1505UL, 8537UL},
1024a238c70SJohn Marino {5231UL, 29819UL},
1034a238c70SJohn Marino {283UL, 1621UL},
1044a238c70SJohn Marino {665UL, 3827UL},
1054a238c70SJohn Marino {32UL, 185UL},
1064a238c70SJohn Marino {2393UL, 13897UL},
1074a238c70SJohn Marino {1879UL, 10960UL},
1084a238c70SJohn Marino {4943UL, 28956UL},
1094a238c70SJohn Marino {409UL, 2406UL},
1104a238c70SJohn Marino {247UL, 1459UL},
1114a238c70SJohn Marino {231UL, 1370UL},
1124a238c70SJohn Marino {3515UL, 20929UL} };
1134a238c70SJohn Marino #if 0
1144a238c70SJohn Marino #define N 8
1154a238c70SJohn Marino int main ()
1164a238c70SJohn Marino {
1174a238c70SJohn Marino unsigned long tab[N];
1184a238c70SJohn Marino int i, n, base;
1194a238c70SJohn Marino mpfr_t x, y;
1204a238c70SJohn Marino mpq_t q1, q2;
1214a238c70SJohn Marino int overflow = 0, base_overflow;
1224a238c70SJohn Marino
1234a238c70SJohn Marino mpfr_init2 (x, 200);
1244a238c70SJohn Marino mpfr_init2 (y, 200);
1254a238c70SJohn Marino mpq_init (q1);
1264a238c70SJohn Marino mpq_init (q2);
1274a238c70SJohn Marino
1284a238c70SJohn Marino for (base = 2 ; base < 63 ; base ++)
1294a238c70SJohn Marino {
1304a238c70SJohn Marino mpfr_set_ui (x, base, MPFR_RNDN);
1314a238c70SJohn Marino mpfr_log2 (x, x, MPFR_RNDN);
1324a238c70SJohn Marino mpfr_ui_div (x, 1, x, MPFR_RNDN);
1334a238c70SJohn Marino printf ("Base: %d x=%e ", base, mpfr_get_d1 (x));
1344a238c70SJohn Marino for (i = 0 ; i < N ; i++)
1354a238c70SJohn Marino {
1364a238c70SJohn Marino mpfr_floor (y, x);
1374a238c70SJohn Marino tab[i] = mpfr_get_ui (y, MPFR_RNDN);
1384a238c70SJohn Marino mpfr_sub (x, x, y, MPFR_RNDN);
1394a238c70SJohn Marino mpfr_ui_div (x, 1, x, MPFR_RNDN);
1404a238c70SJohn Marino }
1414a238c70SJohn Marino for (i = N-1 ; i >= 0 ; i--)
1424a238c70SJohn Marino if (tab[i] != 0)
1434a238c70SJohn Marino break;
1444a238c70SJohn Marino mpq_set_ui (q1, tab[i], 1);
1454a238c70SJohn Marino for (i = i-1 ; i >= 0 ; i--)
1464a238c70SJohn Marino {
1474a238c70SJohn Marino mpq_inv (q1, q1);
1484a238c70SJohn Marino mpq_set_ui (q2, tab[i], 1);
1494a238c70SJohn Marino mpq_add (q1, q1, q2);
1504a238c70SJohn Marino }
1514a238c70SJohn Marino printf("Approx: ", base);
1524a238c70SJohn Marino mpq_out_str (stdout, 10, q1);
1534a238c70SJohn Marino printf (" = %e\n", mpq_get_d (q1) );
1544a238c70SJohn Marino fprintf (stderr, "{");
1554a238c70SJohn Marino mpz_out_str (stderr, 10, mpq_numref (q1));
1564a238c70SJohn Marino fprintf (stderr, "UL, ");
1574a238c70SJohn Marino mpz_out_str (stderr, 10, mpq_denref (q1));
1584a238c70SJohn Marino fprintf (stderr, "UL},\n");
1594a238c70SJohn Marino if (mpz_cmp_ui (mpq_numref (q1), 1<<16-1) >= 0
1604a238c70SJohn Marino || mpz_cmp_ui (mpq_denref (q1), 1<<16-1) >= 0)
1614a238c70SJohn Marino overflow = 1, base_overflow = base;
1624a238c70SJohn Marino }
1634a238c70SJohn Marino
1644a238c70SJohn Marino mpq_clear (q2);
1654a238c70SJohn Marino mpq_clear (q1);
1664a238c70SJohn Marino mpfr_clear (y);
1674a238c70SJohn Marino mpfr_clear (x);
1684a238c70SJohn Marino if (overflow )
1694a238c70SJohn Marino printf ("OVERFLOW for base =%d!\n", base_overflow);
1704a238c70SJohn Marino }
1714a238c70SJohn Marino #endif
1724a238c70SJohn Marino
1734a238c70SJohn Marino
1744a238c70SJohn Marino /* Compatible with any locale, but one still assumes that 'a', 'b', 'c',
1754a238c70SJohn Marino ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like
1764a238c70SJohn Marino in any ASCII-based character set). */
1774a238c70SJohn Marino static int
digit_value_in_base(int c,int base)1784a238c70SJohn Marino digit_value_in_base (int c, int base)
1794a238c70SJohn Marino {
1804a238c70SJohn Marino int digit;
1814a238c70SJohn Marino
1824a238c70SJohn Marino MPFR_ASSERTD (base > 0 && base <= MPFR_MAX_BASE);
1834a238c70SJohn Marino
1844a238c70SJohn Marino if (c >= '0' && c <= '9')
1854a238c70SJohn Marino digit = c - '0';
1864a238c70SJohn Marino else if (c >= 'a' && c <= 'z')
1874a238c70SJohn Marino digit = (base >= 37) ? c - 'a' + 36 : c - 'a' + 10;
1884a238c70SJohn Marino else if (c >= 'A' && c <= 'Z')
1894a238c70SJohn Marino digit = c - 'A' + 10;
1904a238c70SJohn Marino else
1914a238c70SJohn Marino return -1;
1924a238c70SJohn Marino
1934a238c70SJohn Marino return MPFR_LIKELY (digit < base) ? digit : -1;
1944a238c70SJohn Marino }
1954a238c70SJohn Marino
1964a238c70SJohn Marino /* Compatible with any locale, but one still assumes that 'a', 'b', 'c',
1974a238c70SJohn Marino ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like
1984a238c70SJohn Marino in any ASCII-based character set). */
1994a238c70SJohn Marino /* TODO: support EBCDIC. */
2004a238c70SJohn Marino static int
fast_casecmp(const char * s1,const char * s2)2014a238c70SJohn Marino fast_casecmp (const char *s1, const char *s2)
2024a238c70SJohn Marino {
2034a238c70SJohn Marino unsigned char c1, c2;
2044a238c70SJohn Marino
2054a238c70SJohn Marino do
2064a238c70SJohn Marino {
2074a238c70SJohn Marino c2 = *(const unsigned char *) s2++;
2084a238c70SJohn Marino if (c2 == '\0')
2094a238c70SJohn Marino return 0;
2104a238c70SJohn Marino c1 = *(const unsigned char *) s1++;
2114a238c70SJohn Marino if (c1 >= 'A' && c1 <= 'Z')
2124a238c70SJohn Marino c1 = c1 - 'A' + 'a';
2134a238c70SJohn Marino }
2144a238c70SJohn Marino while (c1 == c2);
2154a238c70SJohn Marino return 1;
2164a238c70SJohn Marino }
2174a238c70SJohn Marino
2184a238c70SJohn Marino /* Parse a string and fill pstr.
2194a238c70SJohn Marino Return the advanced ptr too.
2204a238c70SJohn Marino It returns:
2214a238c70SJohn Marino -1 if invalid string,
2224a238c70SJohn Marino 0 if special string (like nan),
2234a238c70SJohn Marino 1 if the string is ok.
2244a238c70SJohn Marino 2 if overflows
2254a238c70SJohn Marino So it doesn't return the ternary value
2264a238c70SJohn Marino BUT if it returns 0 (NAN or INF), the ternary value is also '0'
2274a238c70SJohn Marino (ie NAN and INF are exact) */
2284a238c70SJohn Marino static int
parse_string(mpfr_t x,struct parsed_string * pstr,const char ** string,int base)2294a238c70SJohn Marino parse_string (mpfr_t x, struct parsed_string *pstr,
2304a238c70SJohn Marino const char **string, int base)
2314a238c70SJohn Marino {
2324a238c70SJohn Marino const char *str = *string;
2334a238c70SJohn Marino unsigned char *mant;
2344a238c70SJohn Marino int point;
2354a238c70SJohn Marino int res = -1; /* Invalid input return value */
2364a238c70SJohn Marino const char *prefix_str;
2374a238c70SJohn Marino int decimal_point;
2384a238c70SJohn Marino
2394a238c70SJohn Marino decimal_point = (unsigned char) MPFR_DECIMAL_POINT;
2404a238c70SJohn Marino
2414a238c70SJohn Marino /* Init variable */
2424a238c70SJohn Marino pstr->mantissa = NULL;
2434a238c70SJohn Marino
2444a238c70SJohn Marino /* Optional leading whitespace */
2454a238c70SJohn Marino while (isspace((unsigned char) *str)) str++;
2464a238c70SJohn Marino
2474a238c70SJohn Marino /* An optional sign `+' or `-' */
2484a238c70SJohn Marino pstr->negative = (*str == '-');
2494a238c70SJohn Marino if (*str == '-' || *str == '+')
2504a238c70SJohn Marino str++;
2514a238c70SJohn Marino
2524a238c70SJohn Marino /* Can be case-insensitive NAN */
2534a238c70SJohn Marino if (fast_casecmp (str, "@nan@") == 0)
2544a238c70SJohn Marino {
2554a238c70SJohn Marino str += 5;
2564a238c70SJohn Marino goto set_nan;
2574a238c70SJohn Marino }
2584a238c70SJohn Marino if (base <= 16 && fast_casecmp (str, "nan") == 0)
2594a238c70SJohn Marino {
2604a238c70SJohn Marino str += 3;
2614a238c70SJohn Marino set_nan:
2624a238c70SJohn Marino /* Check for "(dummychars)" */
2634a238c70SJohn Marino if (*str == '(')
2644a238c70SJohn Marino {
2654a238c70SJohn Marino const char *s;
2664a238c70SJohn Marino for (s = str+1 ; *s != ')' ; s++)
2674a238c70SJohn Marino if (!(*s >= 'A' && *s <= 'Z')
2684a238c70SJohn Marino && !(*s >= 'a' && *s <= 'z')
2694a238c70SJohn Marino && !(*s >= '0' && *s <= '9')
2704a238c70SJohn Marino && *s != '_')
2714a238c70SJohn Marino break;
2724a238c70SJohn Marino if (*s == ')')
2734a238c70SJohn Marino str = s+1;
2744a238c70SJohn Marino }
2754a238c70SJohn Marino *string = str;
2764a238c70SJohn Marino MPFR_SET_NAN(x);
2774a238c70SJohn Marino /* MPFR_RET_NAN not used as the return value isn't a ternary value */
2784a238c70SJohn Marino __gmpfr_flags |= MPFR_FLAGS_NAN;
2794a238c70SJohn Marino return 0;
2804a238c70SJohn Marino }
2814a238c70SJohn Marino
2824a238c70SJohn Marino /* Can be case-insensitive INF */
2834a238c70SJohn Marino if (fast_casecmp (str, "@inf@") == 0)
2844a238c70SJohn Marino {
2854a238c70SJohn Marino str += 5;
2864a238c70SJohn Marino goto set_inf;
2874a238c70SJohn Marino }
2884a238c70SJohn Marino if (base <= 16 && fast_casecmp (str, "infinity") == 0)
2894a238c70SJohn Marino {
2904a238c70SJohn Marino str += 8;
2914a238c70SJohn Marino goto set_inf;
2924a238c70SJohn Marino }
2934a238c70SJohn Marino if (base <= 16 && fast_casecmp (str, "inf") == 0)
2944a238c70SJohn Marino {
2954a238c70SJohn Marino str += 3;
2964a238c70SJohn Marino set_inf:
2974a238c70SJohn Marino *string = str;
2984a238c70SJohn Marino MPFR_SET_INF (x);
2994a238c70SJohn Marino (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
3004a238c70SJohn Marino return 0;
3014a238c70SJohn Marino }
3024a238c70SJohn Marino
3034a238c70SJohn Marino /* If base=0 or 16, it may include '0x' prefix */
3044a238c70SJohn Marino prefix_str = NULL;
3054a238c70SJohn Marino if ((base == 0 || base == 16) && str[0]=='0'
3064a238c70SJohn Marino && (str[1]=='x' || str[1] == 'X'))
3074a238c70SJohn Marino {
3084a238c70SJohn Marino prefix_str = str;
3094a238c70SJohn Marino base = 16;
3104a238c70SJohn Marino str += 2;
3114a238c70SJohn Marino }
3124a238c70SJohn Marino /* If base=0 or 2, it may include '0b' prefix */
3134a238c70SJohn Marino if ((base == 0 || base == 2) && str[0]=='0'
3144a238c70SJohn Marino && (str[1]=='b' || str[1] == 'B'))
3154a238c70SJohn Marino {
3164a238c70SJohn Marino prefix_str = str;
3174a238c70SJohn Marino base = 2;
3184a238c70SJohn Marino str += 2;
3194a238c70SJohn Marino }
3204a238c70SJohn Marino /* Else if base=0, we assume decimal base */
3214a238c70SJohn Marino if (base == 0)
3224a238c70SJohn Marino base = 10;
3234a238c70SJohn Marino pstr->base = base;
3244a238c70SJohn Marino
3254a238c70SJohn Marino /* Alloc mantissa */
3264a238c70SJohn Marino pstr->alloc = (size_t) strlen (str) + 1;
3274a238c70SJohn Marino pstr->mantissa = (unsigned char*) (*__gmp_allocate_func) (pstr->alloc);
3284a238c70SJohn Marino
3294a238c70SJohn Marino /* Read mantissa digits */
3304a238c70SJohn Marino parse_begin:
3314a238c70SJohn Marino mant = pstr->mantissa;
3324a238c70SJohn Marino point = 0;
3334a238c70SJohn Marino pstr->exp_base = 0;
3344a238c70SJohn Marino pstr->exp_bin = 0;
3354a238c70SJohn Marino
3364a238c70SJohn Marino for (;;) /* Loop until an invalid character is read */
3374a238c70SJohn Marino {
3384a238c70SJohn Marino int c = (unsigned char) *str++;
3394a238c70SJohn Marino /* The cast to unsigned char is needed because of digit_value_in_base;
3404a238c70SJohn Marino decimal_point uses this convention too. */
3414a238c70SJohn Marino if (c == '.' || c == decimal_point)
3424a238c70SJohn Marino {
3434a238c70SJohn Marino if (MPFR_UNLIKELY(point)) /* Second '.': stop parsing */
3444a238c70SJohn Marino break;
3454a238c70SJohn Marino point = 1;
3464a238c70SJohn Marino continue;
3474a238c70SJohn Marino }
3484a238c70SJohn Marino c = digit_value_in_base (c, base);
3494a238c70SJohn Marino if (c == -1)
3504a238c70SJohn Marino break;
3514a238c70SJohn Marino MPFR_ASSERTN (c >= 0); /* c is representable in an unsigned char */
3524a238c70SJohn Marino *mant++ = (unsigned char) c;
3534a238c70SJohn Marino if (!point)
3544a238c70SJohn Marino pstr->exp_base ++;
3554a238c70SJohn Marino }
3564a238c70SJohn Marino str--; /* The last read character was invalid */
3574a238c70SJohn Marino
3584a238c70SJohn Marino /* Update the # of char in the mantissa */
3594a238c70SJohn Marino pstr->prec = mant - pstr->mantissa;
3604a238c70SJohn Marino /* Check if there are no characters in the mantissa (Invalid argument) */
3614a238c70SJohn Marino if (pstr->prec == 0)
3624a238c70SJohn Marino {
3634a238c70SJohn Marino /* Check if there was a prefix (in such a case, we have to read
3644a238c70SJohn Marino again the mantissa without skipping the prefix)
3654a238c70SJohn Marino The allocated mantissa is still big enough since we will
3664a238c70SJohn Marino read only 0, and we alloc one more char than needed.
3674a238c70SJohn Marino FIXME: Not really friendly. Maybe cleaner code? */
3684a238c70SJohn Marino if (prefix_str != NULL)
3694a238c70SJohn Marino {
3704a238c70SJohn Marino str = prefix_str;
3714a238c70SJohn Marino prefix_str = NULL;
3724a238c70SJohn Marino goto parse_begin;
3734a238c70SJohn Marino }
3744a238c70SJohn Marino goto end;
3754a238c70SJohn Marino }
3764a238c70SJohn Marino
3774a238c70SJohn Marino /* Valid entry */
3784a238c70SJohn Marino res = 1;
3794a238c70SJohn Marino MPFR_ASSERTD (pstr->exp_base >= 0);
3804a238c70SJohn Marino
3814a238c70SJohn Marino /* an optional exponent (e or E, p or P, @) */
3824a238c70SJohn Marino if ( (*str == '@' || (base <= 10 && (*str == 'e' || *str == 'E')))
3834a238c70SJohn Marino && (!isspace((unsigned char) str[1])) )
3844a238c70SJohn Marino {
3854a238c70SJohn Marino char *endptr;
3864a238c70SJohn Marino /* the exponent digits are kept in ASCII */
3874a238c70SJohn Marino mpfr_exp_t sum;
3884a238c70SJohn Marino long read_exp = strtol (str + 1, &endptr, 10);
3894a238c70SJohn Marino if (endptr != str+1)
3904a238c70SJohn Marino str = endptr;
3914a238c70SJohn Marino sum =
3924a238c70SJohn Marino read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) :
3934a238c70SJohn Marino read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) :
3944a238c70SJohn Marino (mpfr_exp_t) read_exp;
3954a238c70SJohn Marino MPFR_SADD_OVERFLOW (sum, sum, pstr->exp_base,
3964a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
3974a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
3984a238c70SJohn Marino res = 2, res = 3);
3994a238c70SJohn Marino /* Since exp_base was positive, read_exp + exp_base can't
4004a238c70SJohn Marino do a negative overflow. */
4014a238c70SJohn Marino MPFR_ASSERTD (res != 3);
4024a238c70SJohn Marino pstr->exp_base = sum;
4034a238c70SJohn Marino }
4044a238c70SJohn Marino else if ((base == 2 || base == 16)
4054a238c70SJohn Marino && (*str == 'p' || *str == 'P')
4064a238c70SJohn Marino && (!isspace((unsigned char) str[1])))
4074a238c70SJohn Marino {
4084a238c70SJohn Marino char *endptr;
4094a238c70SJohn Marino long read_exp = strtol (str + 1, &endptr, 10);
4104a238c70SJohn Marino if (endptr != str+1)
4114a238c70SJohn Marino str = endptr;
4124a238c70SJohn Marino pstr->exp_bin =
4134a238c70SJohn Marino read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) :
4144a238c70SJohn Marino read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) :
4154a238c70SJohn Marino (mpfr_exp_t) read_exp;
4164a238c70SJohn Marino }
4174a238c70SJohn Marino
418*ab6d115fSJohn Marino /* Remove 0's at the beginning and end of mantissa[0..prec-1] */
4194a238c70SJohn Marino mant = pstr->mantissa;
4204a238c70SJohn Marino for ( ; (pstr->prec > 0) && (*mant == 0) ; mant++, pstr->prec--)
4214a238c70SJohn Marino pstr->exp_base--;
4224a238c70SJohn Marino for ( ; (pstr->prec > 0) && (mant[pstr->prec - 1] == 0); pstr->prec--);
4234a238c70SJohn Marino pstr->mant = mant;
4244a238c70SJohn Marino
4254a238c70SJohn Marino /* Check if x = 0 */
4264a238c70SJohn Marino if (pstr->prec == 0)
4274a238c70SJohn Marino {
4284a238c70SJohn Marino MPFR_SET_ZERO (x);
4294a238c70SJohn Marino if (pstr->negative)
4304a238c70SJohn Marino MPFR_SET_NEG(x);
4314a238c70SJohn Marino else
4324a238c70SJohn Marino MPFR_SET_POS(x);
4334a238c70SJohn Marino res = 0;
4344a238c70SJohn Marino }
4354a238c70SJohn Marino
4364a238c70SJohn Marino *string = str;
4374a238c70SJohn Marino end:
4384a238c70SJohn Marino if (pstr->mantissa != NULL && res != 1)
4394a238c70SJohn Marino (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
4404a238c70SJohn Marino return res;
4414a238c70SJohn Marino }
4424a238c70SJohn Marino
4434a238c70SJohn Marino /* Transform a parsed string to a mpfr_t according to the rounding mode
4444a238c70SJohn Marino and the precision of x.
4454a238c70SJohn Marino Returns the ternary value. */
4464a238c70SJohn Marino static int
parsed_string_to_mpfr(mpfr_t x,struct parsed_string * pstr,mpfr_rnd_t rnd)4474a238c70SJohn Marino parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
4484a238c70SJohn Marino {
4494a238c70SJohn Marino mpfr_prec_t prec;
4504a238c70SJohn Marino mpfr_exp_t exp;
4514a238c70SJohn Marino mpfr_exp_t ysize_bits;
4524a238c70SJohn Marino mp_limb_t *y, *result;
4534a238c70SJohn Marino int count, exact;
4544a238c70SJohn Marino size_t pstr_size;
4554a238c70SJohn Marino mp_size_t ysize, real_ysize;
4564a238c70SJohn Marino int res, err;
4574a238c70SJohn Marino MPFR_ZIV_DECL (loop);
4584a238c70SJohn Marino MPFR_TMP_DECL (marker);
4594a238c70SJohn Marino
4604a238c70SJohn Marino /* initialize the working precision */
4614a238c70SJohn Marino prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x));
4624a238c70SJohn Marino
463*ab6d115fSJohn Marino /* compute the value y of the leading characters as long as rounding is not
464*ab6d115fSJohn Marino possible */
4654a238c70SJohn Marino MPFR_TMP_MARK(marker);
4664a238c70SJohn Marino MPFR_ZIV_INIT (loop, prec);
4674a238c70SJohn Marino for (;;)
4684a238c70SJohn Marino {
4694a238c70SJohn Marino /* Set y to the value of the ~prec most significant bits of pstr->mant
4704a238c70SJohn Marino (as long as we guarantee correct rounding, we don't need to get
4714a238c70SJohn Marino exactly prec bits). */
472*ab6d115fSJohn Marino ysize = MPFR_PREC2LIMBS (prec);
4734a238c70SJohn Marino /* prec bits corresponds to ysize limbs */
4744a238c70SJohn Marino ysize_bits = ysize * GMP_NUMB_BITS;
4754a238c70SJohn Marino /* and to ysize_bits >= prec > MPFR_PREC (x) bits */
4764a238c70SJohn Marino y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1);
4774a238c70SJohn Marino y += ysize; /* y has (ysize+1) allocated limbs */
4784a238c70SJohn Marino
4794a238c70SJohn Marino /* pstr_size is the number of characters we read in pstr->mant
4804a238c70SJohn Marino to have at least ysize full limbs.
4814a238c70SJohn Marino We must have base^(pstr_size-1) >= (2^(GMP_NUMB_BITS))^ysize
4824a238c70SJohn Marino (in the worst case, the first digit is one and all others are zero).
4834a238c70SJohn Marino i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base)
4844a238c70SJohn Marino Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 =>
4854a238c70SJohn Marino ysize*GMP_NUMB_BITS can not overflow.
4864a238c70SJohn Marino We compute pstr_size = 1 + ceil(ysize_bits * Num / Den)
4874a238c70SJohn Marino where Num/Den >= 1/log2(base)
4884a238c70SJohn Marino It is not exactly ceil(1/log2(base)) but could be one more (base 2)
4894a238c70SJohn Marino Quite ugly since it tries to avoid overflow:
4904a238c70SJohn Marino let Num = RedInvLog2Table[pstr->base-2][0]
4914a238c70SJohn Marino and Den = RedInvLog2Table[pstr->base-2][1],
4924a238c70SJohn Marino and ysize_bits = a*Den+b,
4934a238c70SJohn Marino then ysize_bits * Num/Den = a*Num + (b * Num)/Den,
4944a238c70SJohn Marino thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den
4954a238c70SJohn Marino */
4964a238c70SJohn Marino {
4974a238c70SJohn Marino unsigned long Num = RedInvLog2Table[pstr->base-2][0];
4984a238c70SJohn Marino unsigned long Den = RedInvLog2Table[pstr->base-2][1];
4994a238c70SJohn Marino pstr_size = ((ysize_bits / Den) * Num)
5004a238c70SJohn Marino + (((ysize_bits % Den) * Num + Den - 1) / Den)
5014a238c70SJohn Marino + 1;
5024a238c70SJohn Marino }
5034a238c70SJohn Marino
5044a238c70SJohn Marino /* since pstr_size corresponds to at least ysize_bits full bits,
5054a238c70SJohn Marino and ysize_bits > prec, the weight of the neglected part of
5064a238c70SJohn Marino pstr->mant (if any) is < ulp(y) < ulp(x) */
5074a238c70SJohn Marino
5084a238c70SJohn Marino /* if the number of wanted characters is more than what we have in
5094a238c70SJohn Marino pstr->mant, round it down */
5104a238c70SJohn Marino if (pstr_size >= pstr->prec)
5114a238c70SJohn Marino pstr_size = pstr->prec;
5124a238c70SJohn Marino MPFR_ASSERTD (pstr_size == (mpfr_exp_t) pstr_size);
5134a238c70SJohn Marino
5144a238c70SJohn Marino /* convert str into binary: note that pstr->mant is big endian,
5154a238c70SJohn Marino thus no offset is needed */
5164a238c70SJohn Marino real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base);
5174a238c70SJohn Marino MPFR_ASSERTD (real_ysize <= ysize+1);
5184a238c70SJohn Marino
519*ab6d115fSJohn Marino /* normalize y: warning we can even get ysize+1 limbs! */
5204a238c70SJohn Marino MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */
5214a238c70SJohn Marino count_leading_zeros (count, y[real_ysize - 1]);
5224a238c70SJohn Marino /* exact means that the number of limbs of the output of mpn_set_str
5234a238c70SJohn Marino is less or equal to ysize */
5244a238c70SJohn Marino exact = real_ysize <= ysize;
5254a238c70SJohn Marino if (exact) /* shift y to the left in that case y should be exact */
5264a238c70SJohn Marino {
5274a238c70SJohn Marino /* we have enough limbs to store {y, real_ysize} */
5284a238c70SJohn Marino /* shift {y, num_limb} for count bits to the left */
5294a238c70SJohn Marino if (count != 0)
5304a238c70SJohn Marino mpn_lshift (y + ysize - real_ysize, y, real_ysize, count);
5314a238c70SJohn Marino if (real_ysize != ysize)
5324a238c70SJohn Marino {
5334a238c70SJohn Marino if (count == 0)
5344a238c70SJohn Marino MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize);
5354a238c70SJohn Marino MPN_ZERO (y, ysize - real_ysize);
5364a238c70SJohn Marino }
5374a238c70SJohn Marino /* for each bit shift decrease exponent of y */
5384a238c70SJohn Marino /* (This should not overflow) */
5394a238c70SJohn Marino exp = - ((ysize - real_ysize) * GMP_NUMB_BITS + count);
5404a238c70SJohn Marino }
5414a238c70SJohn Marino else /* shift y to the right, by doing this we might lose some
5424a238c70SJohn Marino bits from the result of mpn_set_str (in addition to the
5434a238c70SJohn Marino characters neglected from pstr->mant) */
5444a238c70SJohn Marino {
5454a238c70SJohn Marino /* shift {y, num_limb} for (GMP_NUMB_BITS - count) bits
5464a238c70SJohn Marino to the right. FIXME: can we prove that count cannot be zero here,
5474a238c70SJohn Marino since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */
5484a238c70SJohn Marino MPFR_ASSERTD (count != 0);
5494a238c70SJohn Marino exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) ==
5504a238c70SJohn Marino MPFR_LIMB_ZERO;
5514a238c70SJohn Marino /* for each bit shift increase exponent of y */
5524a238c70SJohn Marino exp = GMP_NUMB_BITS - count;
5534a238c70SJohn Marino }
5544a238c70SJohn Marino
555*ab6d115fSJohn Marino /* compute base^(exp_base - pstr_size) on n limbs */
5564a238c70SJohn Marino if (IS_POW2 (pstr->base))
5574a238c70SJohn Marino {
5584a238c70SJohn Marino /* Base: 2, 4, 8, 16, 32 */
5594a238c70SJohn Marino int pow2;
5604a238c70SJohn Marino mpfr_exp_t tmp;
5614a238c70SJohn Marino
5624a238c70SJohn Marino count_leading_zeros (pow2, (mp_limb_t) pstr->base);
5634a238c70SJohn Marino pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */
5644a238c70SJohn Marino MPFR_ASSERTD (0 < pow2 && pow2 <= 5);
5654a238c70SJohn Marino /* exp += pow2 * (pstr->exp_base - pstr_size) + pstr->exp_bin
5664a238c70SJohn Marino with overflow checking
5674a238c70SJohn Marino and check that we can add/substract 2 to exp without overflow */
5684a238c70SJohn Marino MPFR_SADD_OVERFLOW (tmp, pstr->exp_base, -(mpfr_exp_t) pstr_size,
5694a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
5704a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
5714a238c70SJohn Marino goto overflow, goto underflow);
5724a238c70SJohn Marino /* On some FreeBsd/Alpha, LONG_MIN/1 produced an exception
5734a238c70SJohn Marino so we used to check for this before doing the division.
5744a238c70SJohn Marino Since this bug is closed now (Nov 26, 2009), we remove
5754a238c70SJohn Marino that check (http://www.freebsd.org/cgi/query-pr.cgi?pr=72024) */
5764a238c70SJohn Marino if (tmp > 0 && MPFR_EXP_MAX / pow2 <= tmp)
5774a238c70SJohn Marino goto overflow;
5784a238c70SJohn Marino else if (tmp < 0 && MPFR_EXP_MIN / pow2 >= tmp)
5794a238c70SJohn Marino goto underflow;
5804a238c70SJohn Marino tmp *= pow2;
5814a238c70SJohn Marino MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin,
5824a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
5834a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
5844a238c70SJohn Marino goto overflow, goto underflow);
5854a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp, exp, tmp,
5864a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
5874a238c70SJohn Marino MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
5884a238c70SJohn Marino goto overflow, goto underflow);
5894a238c70SJohn Marino result = y;
5904a238c70SJohn Marino err = 0;
5914a238c70SJohn Marino }
5924a238c70SJohn Marino /* case non-power-of-two-base, and pstr->exp_base > pstr_size */
5934a238c70SJohn Marino else if (pstr->exp_base > (mpfr_exp_t) pstr_size)
5944a238c70SJohn Marino {
5954a238c70SJohn Marino mp_limb_t *z;
5964a238c70SJohn Marino mpfr_exp_t exp_z;
5974a238c70SJohn Marino
5984a238c70SJohn Marino result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1);
5994a238c70SJohn Marino
6004a238c70SJohn Marino /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */
6014a238c70SJohn Marino z = y - ysize;
6024a238c70SJohn Marino /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */
6034a238c70SJohn Marino err = mpfr_mpn_exp (z, &exp_z, pstr->base,
6044a238c70SJohn Marino pstr->exp_base - pstr_size, ysize);
6054a238c70SJohn Marino if (err == -2)
6064a238c70SJohn Marino goto overflow;
6074a238c70SJohn Marino exact = exact && (err == -1);
6084a238c70SJohn Marino
6094a238c70SJohn Marino /* If exact is non zero, then z equals exactly the value of the
6104a238c70SJohn Marino pstr_size most significant digits from pstr->mant, i.e., the
6114a238c70SJohn Marino only difference can come from the neglected pstr->prec-pstr_size
6124a238c70SJohn Marino least significant digits of pstr->mant.
6134a238c70SJohn Marino If exact is zero, then z is rounded toward zero with respect
6144a238c70SJohn Marino to that value. */
6154a238c70SJohn Marino
616*ab6d115fSJohn Marino /* multiply(y = 0.mant[0]...mant[pr-1])_base by base^(exp-g):
617*ab6d115fSJohn Marino since both y and z are rounded toward zero, so is "result" */
6184a238c70SJohn Marino mpn_mul_n (result, y, z, ysize);
6194a238c70SJohn Marino
6204a238c70SJohn Marino /* compute the error on the product */
6214a238c70SJohn Marino if (err == -1)
6224a238c70SJohn Marino err = 0;
6234a238c70SJohn Marino err ++;
6244a238c70SJohn Marino
6254a238c70SJohn Marino /* compute the exponent of y */
6264a238c70SJohn Marino /* exp += exp_z + ysize_bits with overflow checking
6274a238c70SJohn Marino and check that we can add/substract 2 to exp without overflow */
6284a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
6294a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
6304a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
6314a238c70SJohn Marino goto overflow, goto underflow);
6324a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp, exp, exp_z,
6334a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
6344a238c70SJohn Marino MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
6354a238c70SJohn Marino goto overflow, goto underflow);
6364a238c70SJohn Marino
6374a238c70SJohn Marino /* normalize result */
6384a238c70SJohn Marino if (MPFR_LIMB_MSB (result[2 * ysize - 1]) == 0)
6394a238c70SJohn Marino {
6404a238c70SJohn Marino mp_limb_t *r = result + ysize - 1;
6414a238c70SJohn Marino mpn_lshift (r, r, ysize + 1, 1);
6424a238c70SJohn Marino /* Overflow checking not needed */
6434a238c70SJohn Marino exp --;
6444a238c70SJohn Marino }
6454a238c70SJohn Marino
6464a238c70SJohn Marino /* if the low ysize limbs of {result, 2*ysize} are all zero,
6474a238c70SJohn Marino then the result is still "exact" (if it was before) */
6484a238c70SJohn Marino exact = exact && (mpn_scan1 (result, 0)
6494a238c70SJohn Marino >= (unsigned long) ysize_bits);
6504a238c70SJohn Marino result += ysize;
6514a238c70SJohn Marino }
6524a238c70SJohn Marino /* case exp_base < pstr_size */
6534a238c70SJohn Marino else if (pstr->exp_base < (mpfr_exp_t) pstr_size)
6544a238c70SJohn Marino {
6554a238c70SJohn Marino mp_limb_t *z;
6564a238c70SJohn Marino mpfr_exp_t exp_z;
6574a238c70SJohn Marino
6584a238c70SJohn Marino result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1);
6594a238c70SJohn Marino
6604a238c70SJohn Marino /* set y to y * K^ysize */
6614a238c70SJohn Marino y = y - ysize; /* we have allocated ysize limbs at y - ysize */
6624a238c70SJohn Marino MPN_ZERO (y, ysize);
6634a238c70SJohn Marino
6644a238c70SJohn Marino /* pstr_size - pstr->exp_base can overflow */
6654a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base,
6664a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
6674a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
6684a238c70SJohn Marino goto underflow, goto overflow);
6694a238c70SJohn Marino
6704a238c70SJohn Marino /* (z, exp_z) = base^(exp_base-pstr_size) */
6714a238c70SJohn Marino z = result + 2*ysize + 1;
6724a238c70SJohn Marino err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize);
673*ab6d115fSJohn Marino /* Since we want y/z rounded toward zero, we must get an upper
674*ab6d115fSJohn Marino bound of z. If err >= 0, the error on z is bounded by 2^err. */
675*ab6d115fSJohn Marino if (err >= 0)
676*ab6d115fSJohn Marino {
677*ab6d115fSJohn Marino mp_limb_t cy;
678*ab6d115fSJohn Marino unsigned long h = err / GMP_NUMB_BITS;
679*ab6d115fSJohn Marino unsigned long l = err - h * GMP_NUMB_BITS;
680*ab6d115fSJohn Marino
681*ab6d115fSJohn Marino if (h >= ysize) /* not enough precision in z */
682*ab6d115fSJohn Marino goto next_loop;
683*ab6d115fSJohn Marino cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l);
684*ab6d115fSJohn Marino if (cy != 0) /* the code below requires z on ysize limbs */
685*ab6d115fSJohn Marino goto next_loop;
686*ab6d115fSJohn Marino }
6874a238c70SJohn Marino exact = exact && (err == -1);
6884a238c70SJohn Marino if (err == -2)
6894a238c70SJohn Marino goto underflow; /* FIXME: Sure? */
6904a238c70SJohn Marino if (err == -1)
6914a238c70SJohn Marino err = 0;
6924a238c70SJohn Marino
6934a238c70SJohn Marino /* compute y / z */
6944a238c70SJohn Marino /* result will be put into result + n, and remainder into result */
6954a238c70SJohn Marino mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y,
6964a238c70SJohn Marino 2 * ysize, z, ysize);
6974a238c70SJohn Marino
6984a238c70SJohn Marino /* exp -= exp_z + ysize_bits with overflow checking
6994a238c70SJohn Marino and check that we can add/substract 2 to exp without overflow */
7004a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
7014a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
7024a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
7034a238c70SJohn Marino goto underflow, goto overflow);
7044a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp, exp, -exp_z,
7054a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
7064a238c70SJohn Marino MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
7074a238c70SJohn Marino goto overflow, goto underflow);
7084a238c70SJohn Marino err += 2;
7094a238c70SJohn Marino /* if the remainder of the division is zero, then the result is
7104a238c70SJohn Marino still "exact" if it was before */
7114a238c70SJohn Marino exact = exact && (mpn_popcount (result, ysize) == 0);
7124a238c70SJohn Marino
7134a238c70SJohn Marino /* normalize result */
7144a238c70SJohn Marino if (result[2 * ysize] == MPFR_LIMB_ONE)
7154a238c70SJohn Marino {
7164a238c70SJohn Marino mp_limb_t *r = result + ysize;
717*ab6d115fSJohn Marino
7184a238c70SJohn Marino exact = exact && ((*r & MPFR_LIMB_ONE) == 0);
7194a238c70SJohn Marino mpn_rshift (r, r, ysize + 1, 1);
7204a238c70SJohn Marino /* Overflow Checking not needed */
7214a238c70SJohn Marino exp ++;
7224a238c70SJohn Marino }
7234a238c70SJohn Marino result += ysize;
7244a238c70SJohn Marino }
7254a238c70SJohn Marino /* case exp_base = pstr_size: no multiplication or division needed */
7264a238c70SJohn Marino else
7274a238c70SJohn Marino {
728*ab6d115fSJohn Marino /* base^(exp-pr) = 1 nothing to compute */
7294a238c70SJohn Marino result = y;
7304a238c70SJohn Marino err = 0;
7314a238c70SJohn Marino }
7324a238c70SJohn Marino
7334a238c70SJohn Marino /* If result is exact, we still have to consider the neglected part
7344a238c70SJohn Marino of the input string. For a directed rounding, in that case we could
7354a238c70SJohn Marino still correctly round, since the neglected part is less than
7364a238c70SJohn Marino one ulp, but that would make the code more complex, and give a
7374a238c70SJohn Marino speedup for rare cases only. */
7384a238c70SJohn Marino exact = exact && (pstr_size == pstr->prec);
7394a238c70SJohn Marino
7404a238c70SJohn Marino /* at this point, result is an approximation rounded toward zero
7414a238c70SJohn Marino of the pstr_size most significant digits of pstr->mant, with
7424a238c70SJohn Marino equality in case exact is non-zero. */
7434a238c70SJohn Marino
7444a238c70SJohn Marino /* test if rounding is possible, and if so exit the loop */
7454a238c70SJohn Marino if (exact || mpfr_can_round_raw (result, ysize,
7464a238c70SJohn Marino (pstr->negative) ? -1 : 1,
7474a238c70SJohn Marino ysize_bits - err - 1,
7484a238c70SJohn Marino MPFR_RNDN, rnd, MPFR_PREC(x)))
7494a238c70SJohn Marino break;
7504a238c70SJohn Marino
751*ab6d115fSJohn Marino next_loop:
7524a238c70SJohn Marino /* update the prec for next loop */
7534a238c70SJohn Marino MPFR_ZIV_NEXT (loop, prec);
7544a238c70SJohn Marino } /* loop */
7554a238c70SJohn Marino MPFR_ZIV_FREE (loop);
7564a238c70SJohn Marino
7574a238c70SJohn Marino /* round y */
7584a238c70SJohn Marino if (mpfr_round_raw (MPFR_MANT (x), result,
7594a238c70SJohn Marino ysize_bits,
7604a238c70SJohn Marino pstr->negative, MPFR_PREC(x), rnd, &res ))
7614a238c70SJohn Marino {
7624a238c70SJohn Marino /* overflow when rounding y */
7634a238c70SJohn Marino MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT;
7644a238c70SJohn Marino /* Overflow Checking not needed */
7654a238c70SJohn Marino exp ++;
7664a238c70SJohn Marino }
7674a238c70SJohn Marino
7684a238c70SJohn Marino if (res == 0) /* fix ternary value */
7694a238c70SJohn Marino {
7704a238c70SJohn Marino exact = exact && (pstr_size == pstr->prec);
7714a238c70SJohn Marino if (!exact)
7724a238c70SJohn Marino res = (pstr->negative) ? 1 : -1;
7734a238c70SJohn Marino }
7744a238c70SJohn Marino
7754a238c70SJohn Marino /* Set sign of x before exp since check_range needs a valid sign */
7764a238c70SJohn Marino (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
7774a238c70SJohn Marino
7784a238c70SJohn Marino /* DO NOT USE MPFR_SET_EXP. The exp may be out of range! */
7794a238c70SJohn Marino MPFR_SADD_OVERFLOW (exp, exp, ysize_bits,
7804a238c70SJohn Marino mpfr_exp_t, mpfr_uexp_t,
7814a238c70SJohn Marino MPFR_EXP_MIN, MPFR_EXP_MAX,
7824a238c70SJohn Marino goto overflow, goto underflow);
7834a238c70SJohn Marino MPFR_EXP (x) = exp;
7844a238c70SJohn Marino res = mpfr_check_range (x, res, rnd);
7854a238c70SJohn Marino goto end;
7864a238c70SJohn Marino
7874a238c70SJohn Marino underflow:
7884a238c70SJohn Marino /* This is called when there is a huge overflow
7894a238c70SJohn Marino (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
7904a238c70SJohn Marino if (rnd == MPFR_RNDN)
7914a238c70SJohn Marino rnd = MPFR_RNDZ;
7924a238c70SJohn Marino res = mpfr_underflow (x, rnd, (pstr->negative) ? -1 : 1);
7934a238c70SJohn Marino goto end;
7944a238c70SJohn Marino
7954a238c70SJohn Marino overflow:
7964a238c70SJohn Marino res = mpfr_overflow (x, rnd, (pstr->negative) ? -1 : 1);
7974a238c70SJohn Marino
7984a238c70SJohn Marino end:
7994a238c70SJohn Marino MPFR_TMP_FREE (marker);
8004a238c70SJohn Marino return res;
8014a238c70SJohn Marino }
8024a238c70SJohn Marino
8034a238c70SJohn Marino static void
free_parsed_string(struct parsed_string * pstr)8044a238c70SJohn Marino free_parsed_string (struct parsed_string *pstr)
8054a238c70SJohn Marino {
8064a238c70SJohn Marino (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
8074a238c70SJohn Marino }
8084a238c70SJohn Marino
8094a238c70SJohn Marino int
mpfr_strtofr(mpfr_t x,const char * string,char ** end,int base,mpfr_rnd_t rnd)8104a238c70SJohn Marino mpfr_strtofr (mpfr_t x, const char *string, char **end, int base,
8114a238c70SJohn Marino mpfr_rnd_t rnd)
8124a238c70SJohn Marino {
8134a238c70SJohn Marino int res;
8144a238c70SJohn Marino struct parsed_string pstr;
8154a238c70SJohn Marino
8164a238c70SJohn Marino /* For base <= 36, parsing is case-insensitive. */
8174a238c70SJohn Marino MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 62));
8184a238c70SJohn Marino
8194a238c70SJohn Marino /* If an error occured, it must return 0 */
8204a238c70SJohn Marino MPFR_SET_ZERO (x);
8214a238c70SJohn Marino MPFR_SET_POS (x);
8224a238c70SJohn Marino
8234a238c70SJohn Marino MPFR_ASSERTN (MPFR_MAX_BASE >= 62);
8244a238c70SJohn Marino res = parse_string (x, &pstr, &string, base);
8254a238c70SJohn Marino /* If res == 0, then it was exact (NAN or INF),
8264a238c70SJohn Marino so it is also the ternary value */
8274a238c70SJohn Marino if (MPFR_UNLIKELY (res == -1)) /* invalid data */
8284a238c70SJohn Marino res = 0; /* x is set to 0, which is exact, thus ternary value is 0 */
8294a238c70SJohn Marino else if (res == 1)
8304a238c70SJohn Marino {
8314a238c70SJohn Marino res = parsed_string_to_mpfr (x, &pstr, rnd);
8324a238c70SJohn Marino free_parsed_string (&pstr);
8334a238c70SJohn Marino }
8344a238c70SJohn Marino else if (res == 2)
8354a238c70SJohn Marino res = mpfr_overflow (x, rnd, (pstr.negative) ? -1 : 1);
8364a238c70SJohn Marino MPFR_ASSERTD (res != 3);
8374a238c70SJohn Marino #if 0
8384a238c70SJohn Marino else if (res == 3)
8394a238c70SJohn Marino {
8404a238c70SJohn Marino /* This is called when there is a huge overflow
8414a238c70SJohn Marino (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
8424a238c70SJohn Marino if (rnd == MPFR_RNDN)
8434a238c70SJohn Marino rnd = MPFR_RNDZ;
8444a238c70SJohn Marino res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1);
8454a238c70SJohn Marino }
8464a238c70SJohn Marino #endif
8474a238c70SJohn Marino
8484a238c70SJohn Marino if (end != NULL)
8494a238c70SJohn Marino *end = (char *) string;
8504a238c70SJohn Marino return res;
8514a238c70SJohn Marino }
852