source/libs/mpfr/mpfr-src/src/strtofr.c

   1 /* mpfr_strtofr -- set a floating-point number from a string
   2
   3 Copyright 2004-2016 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramba projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 3 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #include <stdlib.h> /* For strtol */
  24 #include <ctype.h>  /* For isspace */
  25
  26 #define MPFR_NEED_LONGLONG_H
  27 #include "mpfr-impl.h"
  28
  29 #define MPFR_MAX_BASE 62
  30
  31 struct parsed_string {
  32   int            negative; /* non-zero iff the number is negative */
  33   int            base;     /* base of the string */
  34   unsigned char *mantissa; /* raw significand (without any point) */
  35   unsigned char *mant;     /* stripped significand (without starting and
  36                               ending zeroes). This points inside the area
  37                               allocated for the mantissa field. */
  38   size_t         prec;     /* length of mant (zero for +/-0) */
  39   size_t         alloc;    /* allocation size of mantissa */
  40   mpfr_exp_t     exp_base; /* number of digits before the point */
  41   mpfr_exp_t     exp_bin;  /* exponent in case base=2 or 16, and the pxxx
  42                               format is used (i.e., exponent is given in
  43                               base 10) */
  44 };
  45
  46 /* This table has been generated by the following program.
  47    For 2 <= b <= MPFR_MAX_BASE,
  48    RedInvLog2Table[b-2][0] / RedInvLog2Table[b-2][1]
  49    is an upper approximation of log(2)/log(b).
  50 */
  51 static const unsigned long RedInvLog2Table[MPFR_MAX_BASE-1][2] = {
  52   {1UL, 1UL},
  53   {53UL, 84UL},
  54   {1UL, 2UL},
  55   {4004UL, 9297UL},
  56   {53UL, 137UL},
  57   {2393UL, 6718UL},
  58   {1UL, 3UL},
  59   {665UL, 2108UL},
  60   {4004UL, 13301UL},
  61   {949UL, 3283UL},
  62   {53UL, 190UL},
  63   {5231UL, 19357UL},
  64   {2393UL, 9111UL},
  65   {247UL, 965UL},
  66   {1UL, 4UL},
  67   {4036UL, 16497UL},
  68   {665UL, 2773UL},
  69   {5187UL, 22034UL},
  70   {4004UL, 17305UL},
  71   {51UL, 224UL},
  72   {949UL, 4232UL},
  73   {3077UL, 13919UL},
  74   {53UL, 243UL},
  75   {73UL, 339UL},
  76   {5231UL, 24588UL},
  77   {665UL, 3162UL},
  78   {2393UL, 11504UL},
  79   {4943UL, 24013UL},
  80   {247UL, 1212UL},
  81   {3515UL, 17414UL},
  82   {1UL, 5UL},
  83   {4415UL, 22271UL},
  84   {4036UL, 20533UL},
  85   {263UL, 1349UL},
  86   {665UL, 3438UL},
  87   {1079UL, 5621UL},
  88   {5187UL, 27221UL},
  89   {2288UL, 12093UL},
  90   {4004UL, 21309UL},
  91   {179UL, 959UL},
  92   {51UL, 275UL},
  93   {495UL, 2686UL},
  94   {949UL, 5181UL},
  95   {3621UL, 19886UL},
  96   {3077UL, 16996UL},
  97   {229UL, 1272UL},
  98   {53UL, 296UL},
  99   {109UL, 612UL},
 100   {73UL, 412UL},
 101   {1505UL, 8537UL},
 102   {5231UL, 29819UL},
 103   {283UL, 1621UL},
 104   {665UL, 3827UL},
 105   {32UL, 185UL},
 106   {2393UL, 13897UL},
 107   {1879UL, 10960UL},
 108   {4943UL, 28956UL},
 109   {409UL, 2406UL},
 110   {247UL, 1459UL},
 111   {231UL, 1370UL},
 112   {3515UL, 20929UL} };
 113 #if 0
 114 #define N 8
 115 int main ()
 116 {
 117   unsigned long tab[N];
 118   int i, n, base;
 119   mpfr_t x, y;
 120   mpq_t q1, q2;
 121   int overflow = 0, base_overflow;
 122
 123   mpfr_init2 (x, 200);
 124   mpfr_init2 (y, 200);
 125   mpq_init (q1);
 126   mpq_init (q2);
 127
 128   for (base = 2 ; base < 63 ; base ++)
 129     {
 130       mpfr_set_ui (x, base, MPFR_RNDN);
 131       mpfr_log2 (x, x, MPFR_RNDN);
 132       mpfr_ui_div (x, 1, x, MPFR_RNDN);
 133       printf ("Base: %d x=%e ", base, mpfr_get_d1 (x));
 134       for (i = 0 ; i < N ; i++)
 135         {
 136           mpfr_floor (y, x);
 137           tab[i] = mpfr_get_ui (y, MPFR_RNDN);
 138           mpfr_sub (x, x, y, MPFR_RNDN);
 139           mpfr_ui_div (x, 1, x, MPFR_RNDN);
 140         }
 141       for (i = N-1 ; i >= 0 ; i--)
 142         if (tab[i] != 0)
 143           break;
 144       mpq_set_ui (q1, tab[i], 1);
 145       for (i = i-1 ; i >= 0 ; i--)
 146           {
 147             mpq_inv (q1, q1);
 148             mpq_set_ui (q2, tab[i], 1);
 149             mpq_add (q1, q1, q2);
 150           }
 151       printf("Approx: ", base);
 152       mpq_out_str (stdout, 10, q1);
 153       printf (" = %e\n", mpq_get_d (q1) );
 154       fprintf (stderr, "{");
 155       mpz_out_str (stderr, 10, mpq_numref (q1));
 156       fprintf (stderr, "UL, ");
 157       mpz_out_str (stderr, 10, mpq_denref (q1));
 158       fprintf (stderr, "UL},\n");
 159       if (mpz_cmp_ui (mpq_numref (q1), 1<<16-1) >= 0
 160           || mpz_cmp_ui (mpq_denref (q1), 1<<16-1) >= 0)
 161         overflow = 1, base_overflow = base;
 162     }
 163
 164   mpq_clear (q2);
 165   mpq_clear (q1);
 166   mpfr_clear (y);
 167   mpfr_clear (x);
 168   if (overflow )
 169     printf ("OVERFLOW for base =%d!\n", base_overflow);
 170 }
 171 #endif
 172
 173
 174 /* Compatible with any locale, but one still assumes that 'a', 'b', 'c',
 175    ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like
 176    in any ASCII-based character set). */
 177 static int
 178 digit_value_in_base (int c, int base)
 179 {
 180   int digit;
 181
 182   MPFR_ASSERTD (base > 0 && base <= MPFR_MAX_BASE);
 183
 184   if (c >= '0' && c <= '9')
 185     digit = c - '0';
 186   else if (c >= 'a' && c <= 'z')
 187     digit = (base >= 37) ? c - 'a' + 36 : c - 'a' + 10;
 188   else if (c >= 'A' && c <= 'Z')
 189     digit = c - 'A' + 10;
 190   else
 191     return -1;
 192
 193   return MPFR_LIKELY (digit < base) ? digit : -1;
 194 }
 195
 196 /* Compatible with any locale, but one still assumes that 'a', 'b', 'c',
 197    ..., 'z', and 'A', 'B', 'C', ..., 'Z' are consecutive values (like
 198    in any ASCII-based character set). */
 199 /* TODO: support EBCDIC. */
 200 static int
 201 fast_casecmp (const char *s1, const char *s2)
 202 {
 203   unsigned char c1, c2;
 204
 205   do
 206     {
 207       c2 = *(const unsigned char *) s2++;
 208       if (c2 == '\0')
 209         return 0;
 210       c1 = *(const unsigned char *) s1++;
 211       if (c1 >= 'A' && c1 <= 'Z')
 212         c1 = c1 - 'A' + 'a';
 213     }
 214   while (c1 == c2);
 215   return 1;
 216 }
 217
 218 /* Parse a string and fill pstr.
 219    Return the advanced ptr too.
 220    It returns:
 221       -1 if invalid string,
 222       0 if special string (like nan),
 223       1 if the string is ok.
 224       2 if overflows
 225    So it doesn't return the ternary value
 226    BUT if it returns 0 (NAN or INF), the ternary value is also '0'
 227    (ie NAN and INF are exact) */
 228 static int
 229 parse_string (mpfr_t x, struct parsed_string *pstr,
 230               const char **string, int base)
 231 {
 232   const char *str = *string;
 233   unsigned char *mant;
 234   int point;
 235   int res = -1;  /* Invalid input return value */
 236   const char *prefix_str;
 237   int decimal_point;
 238
 239   decimal_point = (unsigned char) MPFR_DECIMAL_POINT;
 240
 241   /* Init variable */
 242   pstr->mantissa = NULL;
 243
 244   /* Optional leading whitespace */
 245   while (isspace((unsigned char) *str)) str++;
 246
 247   /* An optional sign `+' or `-' */
 248   pstr->negative = (*str == '-');
 249   if (*str == '-' || *str == '+')
 250     str++;
 251
 252   /* Can be case-insensitive NAN */
 253   if (fast_casecmp (str, "@nan@") == 0)
 254     {
 255       str += 5;
 256       goto set_nan;
 257     }
 258   if (base <= 16 && fast_casecmp (str, "nan") == 0)
 259     {
 260       str += 3;
 261     set_nan:
 262       /* Check for "(dummychars)" */
 263       if (*str == '(')
 264         {
 265           const char *s;
 266           for (s = str+1 ; *s != ')' ; s++)
 267             if (!(*s >= 'A' && *s <= 'Z')
 268                 && !(*s >= 'a' && *s <= 'z')
 269                 && !(*s >= '0' && *s <= '9')
 270                 && *s != '_')
 271               break;
 272           if (*s == ')')
 273             str = s+1;
 274         }
 275       *string = str;
 276       MPFR_SET_NAN(x);
 277       /* MPFR_RET_NAN not used as the return value isn't a ternary value */
 278       __gmpfr_flags |= MPFR_FLAGS_NAN;
 279       return 0;
 280     }
 281
 282   /* Can be case-insensitive INF */
 283   if (fast_casecmp (str, "@inf@") == 0)
 284     {
 285       str += 5;
 286       goto set_inf;
 287     }
 288   if (base <= 16 && fast_casecmp (str, "infinity") == 0)
 289     {
 290       str += 8;
 291       goto set_inf;
 292     }
 293   if (base <= 16 && fast_casecmp (str, "inf") == 0)
 294     {
 295       str += 3;
 296     set_inf:
 297       *string = str;
 298       MPFR_SET_INF (x);
 299       (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
 300       return 0;
 301     }
 302
 303   /* If base=0 or 16, it may include '0x' prefix */
 304   prefix_str = NULL;
 305   if ((base == 0 || base == 16) && str[0]=='0'
 306       && (str[1]=='x' || str[1] == 'X'))
 307     {
 308       prefix_str = str;
 309       base = 16;
 310       str += 2;
 311     }
 312   /* If base=0 or 2, it may include '0b' prefix */
 313   if ((base == 0 || base == 2) && str[0]=='0'
 314       && (str[1]=='b' || str[1] == 'B'))
 315     {
 316       prefix_str = str;
 317       base = 2;
 318       str += 2;
 319     }
 320   /* Else if base=0, we assume decimal base */
 321   if (base == 0)
 322     base = 10;
 323   pstr->base = base;
 324
 325   /* Alloc mantissa */
 326   pstr->alloc = (size_t) strlen (str) + 1;
 327   pstr->mantissa = (unsigned char*) (*__gmp_allocate_func) (pstr->alloc);
 328
 329   /* Read mantissa digits */
 330  parse_begin:
 331   mant = pstr->mantissa;
 332   point = 0;
 333   pstr->exp_base = 0;
 334   pstr->exp_bin  = 0;
 335
 336   for (;;) /* Loop until an invalid character is read */
 337     {
 338       int c = (unsigned char) *str++;
 339       /* The cast to unsigned char is needed because of digit_value_in_base;
 340          decimal_point uses this convention too. */
 341       if (c == '.' || c == decimal_point)
 342         {
 343           if (MPFR_UNLIKELY(point)) /* Second '.': stop parsing */
 344             break;
 345           point = 1;
 346           continue;
 347         }
 348       c = digit_value_in_base (c, base);
 349       if (c == -1)
 350         break;
 351       MPFR_ASSERTN (c >= 0); /* c is representable in an unsigned char */
 352       *mant++ = (unsigned char) c;
 353       if (!point)
 354         pstr->exp_base ++;
 355     }
 356   str--; /* The last read character was invalid */
 357
 358   /* Update the # of char in the mantissa */
 359   pstr->prec = mant - pstr->mantissa;
 360   /* Check if there are no characters in the mantissa (Invalid argument) */
 361   if (pstr->prec == 0)
 362     {
 363       /* Check if there was a prefix (in such a case, we have to read
 364          again the mantissa without skipping the prefix)
 365          The allocated mantissa is still big enough since we will
 366          read only 0, and we alloc one more char than needed.
 367          FIXME: Not really friendly. Maybe cleaner code? */
 368       if (prefix_str != NULL)
 369         {
 370           str = prefix_str;
 371           prefix_str = NULL;
 372           goto parse_begin;
 373         }
 374       goto end;
 375     }
 376
 377   /* Valid entry */
 378   res = 1;
 379   MPFR_ASSERTD (pstr->exp_base >= 0);
 380
 381   /* an optional exponent (e or E, p or P, @) */
 382   if ( (*str == '@' || (base <= 10 && (*str == 'e' || *str == 'E')))
 383        && (!isspace((unsigned char) str[1])) )
 384     {
 385       char *endptr;
 386       /* the exponent digits are kept in ASCII */
 387       mpfr_exp_t sum;
 388       long read_exp = strtol (str + 1, &endptr, 10);
 389       if (endptr != str+1)
 390         str = endptr;
 391       sum =
 392         read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) :
 393         read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) :
 394         (mpfr_exp_t) read_exp;
 395       MPFR_SADD_OVERFLOW (sum, sum, pstr->exp_base,
 396                           mpfr_exp_t, mpfr_uexp_t,
 397                           MPFR_EXP_MIN, MPFR_EXP_MAX,
 398                           res = 2, res = 3);
 399       /* Since exp_base was positive, read_exp + exp_base can't
 400          do a negative overflow. */
 401       MPFR_ASSERTD (res != 3);
 402       pstr->exp_base = sum;
 403     }
 404   else if ((base == 2 || base == 16)
 405            && (*str == 'p' || *str == 'P')
 406            && (!isspace((unsigned char) str[1])))
 407     {
 408       char *endptr;
 409       long read_exp = strtol (str + 1, &endptr, 10);
 410       if (endptr != str+1)
 411         str = endptr;
 412       pstr->exp_bin =
 413         read_exp < MPFR_EXP_MIN ? (str = endptr, MPFR_EXP_MIN) :
 414         read_exp > MPFR_EXP_MAX ? (str = endptr, MPFR_EXP_MAX) :
 415         (mpfr_exp_t) read_exp;
 416     }
 417
 418   /* Remove 0's at the beginning and end of mantissa[0..prec-1] */
 419   mant = pstr->mantissa;
 420   for ( ; (pstr->prec > 0) && (*mant == 0) ; mant++, pstr->prec--)
 421     pstr->exp_base--;
 422   for ( ; (pstr->prec > 0) && (mant[pstr->prec - 1] == 0); pstr->prec--);
 423   pstr->mant = mant;
 424
 425   /* Check if x = 0 */
 426   if (pstr->prec == 0)
 427     {
 428       MPFR_SET_ZERO (x);
 429       if (pstr->negative)
 430         MPFR_SET_NEG(x);
 431       else
 432         MPFR_SET_POS(x);
 433       res = 0;
 434     }
 435
 436   *string = str;
 437  end:
 438   if (pstr->mantissa != NULL && res != 1)
 439     (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
 440   return res;
 441 }
 442
 443 /* Transform a parsed string to a mpfr_t according to the rounding mode
 444    and the precision of x.
 445    Returns the ternary value. */
 446 static int
 447 parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mpfr_rnd_t rnd)
 448 {
 449   mpfr_prec_t prec;
 450   mpfr_exp_t exp;
 451   mpfr_exp_t ysize_bits;
 452   mp_limb_t *y, *result;
 453   int count, exact;
 454   size_t pstr_size;
 455   mp_size_t ysize, real_ysize;
 456   int res, err;
 457   MPFR_ZIV_DECL (loop);
 458   MPFR_TMP_DECL (marker);
 459
 460   /* initialize the working precision */
 461   prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x));
 462
 463   /* compute the value y of the leading characters as long as rounding is not
 464      possible */
 465   MPFR_TMP_MARK(marker);
 466   MPFR_ZIV_INIT (loop, prec);
 467   for (;;)
 468     {
 469       /* Set y to the value of the ~prec most significant bits of pstr->mant
 470          (as long as we guarantee correct rounding, we don't need to get
 471          exactly prec bits). */
 472       ysize = MPFR_PREC2LIMBS (prec);
 473       /* prec bits corresponds to ysize limbs */
 474       ysize_bits = ysize * GMP_NUMB_BITS;
 475       /* and to ysize_bits >= prec > MPFR_PREC (x) bits */
 476       /* we need to allocate one more limb to work around bug
 477          https://gmplib.org/list-archives/gmp-bugs/2013-December/003267.html */
 478       y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 2);
 479       y += ysize; /* y has (ysize+2) allocated limbs */
 480
 481       /* pstr_size is the number of characters we read in pstr->mant
 482          to have at least ysize full limbs.
 483          We must have base^(pstr_size-1) >= (2^(GMP_NUMB_BITS))^ysize
 484          (in the worst case, the first digit is one and all others are zero).
 485          i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base)
 486           Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 =>
 487           ysize*GMP_NUMB_BITS can not overflow.
 488          We compute pstr_size = 1 + ceil(ysize_bits * Num / Den)
 489           where Num/Den >= 1/log2(base)
 490          It is not exactly ceil(1/log2(base)) but could be one more (base 2)
 491          Quite ugly since it tries to avoid overflow:
 492          let Num = RedInvLog2Table[pstr->base-2][0]
 493          and Den = RedInvLog2Table[pstr->base-2][1],
 494          and ysize_bits = a*Den+b,
 495          then ysize_bits * Num/Den = a*Num + (b * Num)/Den,
 496          thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den
 497       */
 498       {
 499         unsigned long Num = RedInvLog2Table[pstr->base-2][0];
 500         unsigned long Den = RedInvLog2Table[pstr->base-2][1];
 501         pstr_size = ((ysize_bits / Den) * Num)
 502           + (((ysize_bits % Den) * Num + Den - 1) / Den)
 503           + 1;
 504       }
 505
 506       /* since pstr_size corresponds to at least ysize_bits full bits,
 507          and ysize_bits > prec, the weight of the neglected part of
 508          pstr->mant (if any) is < ulp(y) < ulp(x) */
 509
 510       /* if the number of wanted characters is more than what we have in
 511          pstr->mant, round it down */
 512       if (pstr_size >= pstr->prec)
 513         pstr_size = pstr->prec;
 514       MPFR_ASSERTD (pstr_size == (mpfr_exp_t) pstr_size);
 515
 516       /* convert str into binary: note that pstr->mant is big endian,
 517          thus no offset is needed */
 518       real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base);
 519       MPFR_ASSERTD (real_ysize <= ysize+1);
 520
 521       /* normalize y: warning we can even get ysize+1 limbs! */
 522       MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */
 523       count_leading_zeros (count, y[real_ysize - 1]);
 524       /* exact means that the number of limbs of the output of mpn_set_str
 525          is less or equal to ysize */
 526       exact = real_ysize <= ysize;
 527       if (exact) /* shift y to the left in that case y should be exact */
 528         {
 529           /* we have enough limbs to store {y, real_ysize} */
 530           /* shift {y, num_limb} for count bits to the left */
 531           if (count != 0)
 532             mpn_lshift (y + ysize - real_ysize, y, real_ysize, count);
 533           if (real_ysize != ysize)
 534             {
 535               if (count == 0)
 536                 MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize);
 537               MPN_ZERO (y, ysize - real_ysize);
 538             }
 539           /* for each bit shift decrease exponent of y */
 540           /* (This should not overflow) */
 541           exp = - ((ysize - real_ysize) * GMP_NUMB_BITS + count);
 542         }
 543       else  /* shift y to the right, by doing this we might lose some
 544                bits from the result of mpn_set_str (in addition to the
 545                characters neglected from pstr->mant) */
 546         {
 547           /* shift {y, num_limb} for (GMP_NUMB_BITS - count) bits
 548              to the right. FIXME: can we prove that count cannot be zero here,
 549              since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */
 550           MPFR_ASSERTD (count != 0);
 551           exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) ==
 552             MPFR_LIMB_ZERO;
 553           /* for each bit shift increase exponent of y */
 554           exp = GMP_NUMB_BITS - count;
 555         }
 556
 557       /* compute base^(exp_base - pstr_size) on n limbs */
 558       if (IS_POW2 (pstr->base))
 559         {
 560           /* Base: 2, 4, 8, 16, 32 */
 561           int pow2;
 562           mpfr_exp_t tmp;
 563
 564           count_leading_zeros (pow2, (mp_limb_t) pstr->base);
 565           pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */
 566           MPFR_ASSERTD (0 < pow2 && pow2 <= 5);
 567           /* exp += pow2 * (pstr->exp_base - pstr_size) + pstr->exp_bin
 568              with overflow checking
 569              and check that we can add/subtract 2 to exp without overflow */
 570           MPFR_SADD_OVERFLOW (tmp, pstr->exp_base, -(mpfr_exp_t) pstr_size,
 571                               mpfr_exp_t, mpfr_uexp_t,
 572                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 573                               goto overflow, goto underflow);
 574           /* On some FreeBsd/Alpha, LONG_MIN/1 produced an exception
 575              so we used to check for this before doing the division.
 576              Since this bug is closed now (Nov 26, 2009), we remove
 577              that check (http://www.freebsd.org/cgi/query-pr.cgi?pr=72024) */
 578           if (tmp > 0 && MPFR_EXP_MAX / pow2 <= tmp)
 579             goto overflow;
 580           else if (tmp < 0 && MPFR_EXP_MIN / pow2 >= tmp)
 581             goto underflow;
 582           tmp *= pow2;
 583           MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin,
 584                               mpfr_exp_t, mpfr_uexp_t,
 585                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 586                               goto overflow, goto underflow);
 587           MPFR_SADD_OVERFLOW (exp, exp, tmp,
 588                               mpfr_exp_t, mpfr_uexp_t,
 589                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 590                               goto overflow, goto underflow);
 591           result = y;
 592           err = 0;
 593         }
 594       /* case non-power-of-two-base, and pstr->exp_base > pstr_size */
 595       else if (pstr->exp_base > (mpfr_exp_t) pstr_size)
 596         {
 597           mp_limb_t *z;
 598           mpfr_exp_t exp_z;
 599
 600           result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1);
 601
 602           /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */
 603           z = y - ysize;
 604           /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */
 605           err = mpfr_mpn_exp (z, &exp_z, pstr->base,
 606                               pstr->exp_base - pstr_size, ysize);
 607           if (err == -2)
 608             goto overflow;
 609           exact = exact && (err == -1);
 610
 611           /* If exact is non zero, then z equals exactly the value of the
 612              pstr_size most significant digits from pstr->mant, i.e., the
 613              only difference can come from the neglected pstr->prec-pstr_size
 614              least significant digits of pstr->mant.
 615              If exact is zero, then z is rounded toward zero with respect
 616              to that value. */
 617
 618           /* multiply(y = 0.mant[0]...mant[pr-1])_base by base^(exp-g):
 619              since both y and z are rounded toward zero, so is "result" */
 620           mpn_mul_n (result, y, z, ysize);
 621
 622           /* compute the error on the product */
 623           if (err == -1)
 624             err = 0;
 625           err ++;
 626
 627           /* compute the exponent of y */
 628           /* exp += exp_z + ysize_bits with overflow checking
 629              and check that we can add/subtract 2 to exp without overflow */
 630           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 631                               mpfr_exp_t, mpfr_uexp_t,
 632                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 633                               goto overflow, goto underflow);
 634           MPFR_SADD_OVERFLOW (exp, exp, exp_z,
 635                               mpfr_exp_t, mpfr_uexp_t,
 636                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 637                               goto overflow, goto underflow);
 638
 639           /* normalize result */
 640           if (MPFR_LIMB_MSB (result[2 * ysize - 1]) == 0)
 641             {
 642               mp_limb_t *r = result + ysize - 1;
 643               mpn_lshift (r, r, ysize + 1, 1);
 644               /* Overflow checking not needed */
 645               exp --;
 646             }
 647
 648           /* if the low ysize limbs of {result, 2*ysize} are all zero,
 649              then the result is still "exact" (if it was before) */
 650           exact = exact && (mpn_scan1 (result, 0)
 651                             >= (unsigned long) ysize_bits);
 652           result += ysize;
 653         }
 654       /* case exp_base < pstr_size */
 655       else if (pstr->exp_base < (mpfr_exp_t) pstr_size)
 656         {
 657           mp_limb_t *z;
 658           mpfr_exp_t exp_z;
 659
 660           result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1);
 661
 662           /* set y to y * K^ysize */
 663           y = y - ysize;  /* we have allocated ysize limbs at y - ysize */
 664           MPN_ZERO (y, ysize);
 665
 666           /* pstr_size - pstr->exp_base can overflow */
 667           MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base,
 668                               mpfr_exp_t, mpfr_uexp_t,
 669                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 670                               goto underflow, goto overflow);
 671
 672           /* (z, exp_z) = base^(exp_base-pstr_size) */
 673           z = result + 2*ysize + 1;
 674           err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize);
 675           /* Since we want y/z rounded toward zero, we must get an upper
 676              bound of z. If err >= 0, the error on z is bounded by 2^err. */
 677           if (err >= 0)
 678             {
 679               mp_limb_t cy;
 680               unsigned long h = err / GMP_NUMB_BITS;
 681               unsigned long l = err - h * GMP_NUMB_BITS;
 682
 683               if (h >= ysize) /* not enough precision in z */
 684                 goto next_loop;
 685               cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l);
 686               if (cy != 0) /* the code below requires z on ysize limbs */
 687                 goto next_loop;
 688             }
 689           exact = exact && (err == -1);
 690           if (err == -2)
 691             goto underflow; /* FIXME: Sure? */
 692           if (err == -1)
 693             err = 0;
 694
 695           /* compute y / z */
 696           /* result will be put into result + n, and remainder into result */
 697           mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y,
 698                        2 * ysize, z, ysize);
 699
 700           /* exp -= exp_z + ysize_bits with overflow checking
 701              and check that we can add/subtract 2 to exp without overflow */
 702           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 703                               mpfr_exp_t, mpfr_uexp_t,
 704                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 705                               goto underflow, goto overflow);
 706           MPFR_SADD_OVERFLOW (exp, exp, -exp_z,
 707                               mpfr_exp_t, mpfr_uexp_t,
 708                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 709                               goto overflow, goto underflow);
 710           err += 2;
 711           /* if the remainder of the division is zero, then the result is
 712              still "exact" if it was before */
 713           exact = exact && (mpn_popcount (result, ysize) == 0);
 714
 715           /* normalize result */
 716           if (result[2 * ysize] == MPFR_LIMB_ONE)
 717             {
 718               mp_limb_t *r = result + ysize;
 719
 720               exact = exact && ((*r & MPFR_LIMB_ONE) == 0);
 721               mpn_rshift (r, r, ysize + 1, 1);
 722               /* Overflow Checking not needed */
 723               exp ++;
 724             }
 725           result += ysize;
 726         }
 727       /* case exp_base = pstr_size: no multiplication or division needed */
 728       else
 729         {
 730           /* base^(exp-pr) = 1             nothing to compute */
 731           result = y;
 732           err = 0;
 733         }
 734
 735       /* If result is exact, we still have to consider the neglected part
 736          of the input string. For a directed rounding, in that case we could
 737          still correctly round, since the neglected part is less than
 738          one ulp, but that would make the code more complex, and give a
 739          speedup for rare cases only. */
 740       exact = exact && (pstr_size == pstr->prec);
 741
 742       /* at this point, result is an approximation rounded toward zero
 743          of the pstr_size most significant digits of pstr->mant, with
 744          equality in case exact is non-zero. */
 745
 746       /* test if rounding is possible, and if so exit the loop */
 747       if (exact || mpfr_can_round_raw (result, ysize,
 748                                        (pstr->negative) ? -1 : 1,
 749                                        ysize_bits - err - 1,
 750                                        MPFR_RNDN, rnd, MPFR_PREC(x)))
 751         break;
 752
 753     next_loop:
 754       /* update the prec for next loop */
 755       MPFR_ZIV_NEXT (loop, prec);
 756     } /* loop */
 757   MPFR_ZIV_FREE (loop);
 758
 759   /* round y */
 760   if (mpfr_round_raw (MPFR_MANT (x), result,
 761                       ysize_bits,
 762                       pstr->negative, MPFR_PREC(x), rnd, &res ))
 763     {
 764       /* overflow when rounding y */
 765       MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT;
 766       /* Overflow Checking not needed */
 767       exp ++;
 768     }
 769
 770   if (res == 0) /* fix ternary value */
 771     {
 772       exact = exact && (pstr_size == pstr->prec);
 773       if (!exact)
 774         res = (pstr->negative) ? 1 : -1;
 775     }
 776
 777   /* Set sign of x before exp since check_range needs a valid sign */
 778   (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
 779
 780   /* DO NOT USE MPFR_SET_EXP. The exp may be out of range! */
 781   MPFR_SADD_OVERFLOW (exp, exp, ysize_bits,
 782                       mpfr_exp_t, mpfr_uexp_t,
 783                       MPFR_EXP_MIN, MPFR_EXP_MAX,
 784                       goto overflow, goto underflow);
 785   MPFR_EXP (x) = exp;
 786   res = mpfr_check_range (x, res, rnd);
 787   goto end;
 788
 789  underflow:
 790   /* This is called when there is a huge overflow
 791      (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 792   if (rnd == MPFR_RNDN)
 793     rnd = MPFR_RNDZ;
 794   res = mpfr_underflow (x, rnd, (pstr->negative) ? -1 : 1);
 795   goto end;
 796
 797  overflow:
 798   res = mpfr_overflow (x, rnd, (pstr->negative) ? -1 : 1);
 799
 800  end:
 801   MPFR_TMP_FREE (marker);
 802   return res;
 803 }
 804
 805 static void
 806 free_parsed_string (struct parsed_string *pstr)
 807 {
 808   (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
 809 }
 810
 811 int
 812 mpfr_strtofr (mpfr_t x, const char *string, char **end, int base,
 813               mpfr_rnd_t rnd)
 814 {
 815   int res;
 816   struct parsed_string pstr;
 817
 818   /* For base <= 36, parsing is case-insensitive. */
 819   MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 62));
 820
 821   /* If an error occured, it must return 0 */
 822   MPFR_SET_ZERO (x);
 823   MPFR_SET_POS (x);
 824
 825   MPFR_ASSERTN (MPFR_MAX_BASE >= 62);
 826   res = parse_string (x, &pstr, &string, base);
 827   /* If res == 0, then it was exact (NAN or INF),
 828      so it is also the ternary value */
 829   if (MPFR_UNLIKELY (res == -1))  /* invalid data */
 830     res = 0;  /* x is set to 0, which is exact, thus ternary value is 0 */
 831   else if (res == 1)
 832     {
 833       res = parsed_string_to_mpfr (x, &pstr, rnd);
 834       free_parsed_string (&pstr);
 835     }
 836   else if (res == 2)
 837     res = mpfr_overflow (x, rnd, (pstr.negative) ? -1 : 1);
 838   MPFR_ASSERTD (res != 3);
 839 #if 0
 840   else if (res == 3)
 841     {
 842       /* This is called when there is a huge overflow
 843          (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 844       if (rnd == MPFR_RNDN)
 845         rnd = MPFR_RNDZ;
 846       res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1);
 847     }
 848 #endif
 849
 850   if (end != NULL)
 851     *end = (char *) string;
 852   return res;
 853 }