contrib/mpfr/src/strtofr.c

   1 /* mpfr_strtofr -- set a floating-point number from a string
   2
   3 Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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       y = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1);
 477       y += ysize; /* y has (ysize+1) allocated limbs */
 478
 479       /* pstr_size is the number of characters we read in pstr->mant
 480          to have at least ysize full limbs.
 481          We must have base^(pstr_size-1) >= (2^(GMP_NUMB_BITS))^ysize
 482          (in the worst case, the first digit is one and all others are zero).
 483          i.e., pstr_size >= 1 + ysize*GMP_NUMB_BITS/log2(base)
 484           Since ysize ~ prec/GMP_NUMB_BITS and prec < Umax/2 =>
 485           ysize*GMP_NUMB_BITS can not overflow.
 486          We compute pstr_size = 1 + ceil(ysize_bits * Num / Den)
 487           where Num/Den >= 1/log2(base)
 488          It is not exactly ceil(1/log2(base)) but could be one more (base 2)
 489          Quite ugly since it tries to avoid overflow:
 490          let Num = RedInvLog2Table[pstr->base-2][0]
 491          and Den = RedInvLog2Table[pstr->base-2][1],
 492          and ysize_bits = a*Den+b,
 493          then ysize_bits * Num/Den = a*Num + (b * Num)/Den,
 494          thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den
 495       */
 496       {
 497         unsigned long Num = RedInvLog2Table[pstr->base-2][0];
 498         unsigned long Den = RedInvLog2Table[pstr->base-2][1];
 499         pstr_size = ((ysize_bits / Den) * Num)
 500           + (((ysize_bits % Den) * Num + Den - 1) / Den)
 501           + 1;
 502       }
 503
 504       /* since pstr_size corresponds to at least ysize_bits full bits,
 505          and ysize_bits > prec, the weight of the neglected part of
 506          pstr->mant (if any) is < ulp(y) < ulp(x) */
 507
 508       /* if the number of wanted characters is more than what we have in
 509          pstr->mant, round it down */
 510       if (pstr_size >= pstr->prec)
 511         pstr_size = pstr->prec;
 512       MPFR_ASSERTD (pstr_size == (mpfr_exp_t) pstr_size);
 513
 514       /* convert str into binary: note that pstr->mant is big endian,
 515          thus no offset is needed */
 516       real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base);
 517       MPFR_ASSERTD (real_ysize <= ysize+1);
 518
 519       /* normalize y: warning we can even get ysize+1 limbs! */
 520       MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */
 521       count_leading_zeros (count, y[real_ysize - 1]);
 522       /* exact means that the number of limbs of the output of mpn_set_str
 523          is less or equal to ysize */
 524       exact = real_ysize <= ysize;
 525       if (exact) /* shift y to the left in that case y should be exact */
 526         {
 527           /* we have enough limbs to store {y, real_ysize} */
 528           /* shift {y, num_limb} for count bits to the left */
 529           if (count != 0)
 530             mpn_lshift (y + ysize - real_ysize, y, real_ysize, count);
 531           if (real_ysize != ysize)
 532             {
 533               if (count == 0)
 534                 MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize);
 535               MPN_ZERO (y, ysize - real_ysize);
 536             }
 537           /* for each bit shift decrease exponent of y */
 538           /* (This should not overflow) */
 539           exp = - ((ysize - real_ysize) * GMP_NUMB_BITS + count);
 540         }
 541       else  /* shift y to the right, by doing this we might lose some
 542                bits from the result of mpn_set_str (in addition to the
 543                characters neglected from pstr->mant) */
 544         {
 545           /* shift {y, num_limb} for (GMP_NUMB_BITS - count) bits
 546              to the right. FIXME: can we prove that count cannot be zero here,
 547              since mpn_rshift does not accept a shift of GMP_NUMB_BITS? */
 548           MPFR_ASSERTD (count != 0);
 549           exact = mpn_rshift (y, y, real_ysize, GMP_NUMB_BITS - count) ==
 550             MPFR_LIMB_ZERO;
 551           /* for each bit shift increase exponent of y */
 552           exp = GMP_NUMB_BITS - count;
 553         }
 554
 555       /* compute base^(exp_base - pstr_size) on n limbs */
 556       if (IS_POW2 (pstr->base))
 557         {
 558           /* Base: 2, 4, 8, 16, 32 */
 559           int pow2;
 560           mpfr_exp_t tmp;
 561
 562           count_leading_zeros (pow2, (mp_limb_t) pstr->base);
 563           pow2 = GMP_NUMB_BITS - pow2 - 1; /* base = 2^pow2 */
 564           MPFR_ASSERTD (0 < pow2 && pow2 <= 5);
 565           /* exp += pow2 * (pstr->exp_base - pstr_size) + pstr->exp_bin
 566              with overflow checking
 567              and check that we can add/substract 2 to exp without overflow */
 568           MPFR_SADD_OVERFLOW (tmp, pstr->exp_base, -(mpfr_exp_t) pstr_size,
 569                               mpfr_exp_t, mpfr_uexp_t,
 570                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 571                               goto overflow, goto underflow);
 572           /* On some FreeBsd/Alpha, LONG_MIN/1 produced an exception
 573              so we used to check for this before doing the division.
 574              Since this bug is closed now (Nov 26, 2009), we remove
 575              that check (http://www.freebsd.org/cgi/query-pr.cgi?pr=72024) */
 576           if (tmp > 0 && MPFR_EXP_MAX / pow2 <= tmp)
 577             goto overflow;
 578           else if (tmp < 0 && MPFR_EXP_MIN / pow2 >= tmp)
 579             goto underflow;
 580           tmp *= pow2;
 581           MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin,
 582                               mpfr_exp_t, mpfr_uexp_t,
 583                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 584                               goto overflow, goto underflow);
 585           MPFR_SADD_OVERFLOW (exp, exp, tmp,
 586                               mpfr_exp_t, mpfr_uexp_t,
 587                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 588                               goto overflow, goto underflow);
 589           result = y;
 590           err = 0;
 591         }
 592       /* case non-power-of-two-base, and pstr->exp_base > pstr_size */
 593       else if (pstr->exp_base > (mpfr_exp_t) pstr_size)
 594         {
 595           mp_limb_t *z;
 596           mpfr_exp_t exp_z;
 597
 598           result = MPFR_TMP_LIMBS_ALLOC (2 * ysize + 1);
 599
 600           /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */
 601           z = y - ysize;
 602           /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */
 603           err = mpfr_mpn_exp (z, &exp_z, pstr->base,
 604                               pstr->exp_base - pstr_size, ysize);
 605           if (err == -2)
 606             goto overflow;
 607           exact = exact && (err == -1);
 608
 609           /* If exact is non zero, then z equals exactly the value of the
 610              pstr_size most significant digits from pstr->mant, i.e., the
 611              only difference can come from the neglected pstr->prec-pstr_size
 612              least significant digits of pstr->mant.
 613              If exact is zero, then z is rounded toward zero with respect
 614              to that value. */
 615
 616           /* multiply(y = 0.mant[0]...mant[pr-1])_base by base^(exp-g):
 617              since both y and z are rounded toward zero, so is "result" */
 618           mpn_mul_n (result, y, z, ysize);
 619
 620           /* compute the error on the product */
 621           if (err == -1)
 622             err = 0;
 623           err ++;
 624
 625           /* compute the exponent of y */
 626           /* exp += exp_z + ysize_bits with overflow checking
 627              and check that we can add/substract 2 to exp without overflow */
 628           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 629                               mpfr_exp_t, mpfr_uexp_t,
 630                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 631                               goto overflow, goto underflow);
 632           MPFR_SADD_OVERFLOW (exp, exp, exp_z,
 633                               mpfr_exp_t, mpfr_uexp_t,
 634                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 635                               goto overflow, goto underflow);
 636
 637           /* normalize result */
 638           if (MPFR_LIMB_MSB (result[2 * ysize - 1]) == 0)
 639             {
 640               mp_limb_t *r = result + ysize - 1;
 641               mpn_lshift (r, r, ysize + 1, 1);
 642               /* Overflow checking not needed */
 643               exp --;
 644             }
 645
 646           /* if the low ysize limbs of {result, 2*ysize} are all zero,
 647              then the result is still "exact" (if it was before) */
 648           exact = exact && (mpn_scan1 (result, 0)
 649                             >= (unsigned long) ysize_bits);
 650           result += ysize;
 651         }
 652       /* case exp_base < pstr_size */
 653       else if (pstr->exp_base < (mpfr_exp_t) pstr_size)
 654         {
 655           mp_limb_t *z;
 656           mpfr_exp_t exp_z;
 657
 658           result = MPFR_TMP_LIMBS_ALLOC (3 * ysize + 1);
 659
 660           /* set y to y * K^ysize */
 661           y = y - ysize;  /* we have allocated ysize limbs at y - ysize */
 662           MPN_ZERO (y, ysize);
 663
 664           /* pstr_size - pstr->exp_base can overflow */
 665           MPFR_SADD_OVERFLOW (exp_z, (mpfr_exp_t) pstr_size, -pstr->exp_base,
 666                               mpfr_exp_t, mpfr_uexp_t,
 667                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 668                               goto underflow, goto overflow);
 669
 670           /* (z, exp_z) = base^(exp_base-pstr_size) */
 671           z = result + 2*ysize + 1;
 672           err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize);
 673           /* Since we want y/z rounded toward zero, we must get an upper
 674              bound of z. If err >= 0, the error on z is bounded by 2^err. */
 675           if (err >= 0)
 676             {
 677               mp_limb_t cy;
 678               unsigned long h = err / GMP_NUMB_BITS;
 679               unsigned long l = err - h * GMP_NUMB_BITS;
 680
 681               if (h >= ysize) /* not enough precision in z */
 682                 goto next_loop;
 683               cy = mpn_add_1 (z, z, ysize - h, MPFR_LIMB_ONE << l);
 684               if (cy != 0) /* the code below requires z on ysize limbs */
 685                 goto next_loop;
 686             }
 687           exact = exact && (err == -1);
 688           if (err == -2)
 689             goto underflow; /* FIXME: Sure? */
 690           if (err == -1)
 691             err = 0;
 692
 693           /* compute y / z */
 694           /* result will be put into result + n, and remainder into result */
 695           mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y,
 696                        2 * ysize, z, ysize);
 697
 698           /* exp -= exp_z + ysize_bits with overflow checking
 699              and check that we can add/substract 2 to exp without overflow */
 700           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 701                               mpfr_exp_t, mpfr_uexp_t,
 702                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 703                               goto underflow, goto overflow);
 704           MPFR_SADD_OVERFLOW (exp, exp, -exp_z,
 705                               mpfr_exp_t, mpfr_uexp_t,
 706                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 707                               goto overflow, goto underflow);
 708           err += 2;
 709           /* if the remainder of the division is zero, then the result is
 710              still "exact" if it was before */
 711           exact = exact && (mpn_popcount (result, ysize) == 0);
 712
 713           /* normalize result */
 714           if (result[2 * ysize] == MPFR_LIMB_ONE)
 715             {
 716               mp_limb_t *r = result + ysize;
 717
 718               exact = exact && ((*r & MPFR_LIMB_ONE) == 0);
 719               mpn_rshift (r, r, ysize + 1, 1);
 720               /* Overflow Checking not needed */
 721               exp ++;
 722             }
 723           result += ysize;
 724         }
 725       /* case exp_base = pstr_size: no multiplication or division needed */
 726       else
 727         {
 728           /* base^(exp-pr) = 1             nothing to compute */
 729           result = y;
 730           err = 0;
 731         }
 732
 733       /* If result is exact, we still have to consider the neglected part
 734          of the input string. For a directed rounding, in that case we could
 735          still correctly round, since the neglected part is less than
 736          one ulp, but that would make the code more complex, and give a
 737          speedup for rare cases only. */
 738       exact = exact && (pstr_size == pstr->prec);
 739
 740       /* at this point, result is an approximation rounded toward zero
 741          of the pstr_size most significant digits of pstr->mant, with
 742          equality in case exact is non-zero. */
 743
 744       /* test if rounding is possible, and if so exit the loop */
 745       if (exact || mpfr_can_round_raw (result, ysize,
 746                                        (pstr->negative) ? -1 : 1,
 747                                        ysize_bits - err - 1,
 748                                        MPFR_RNDN, rnd, MPFR_PREC(x)))
 749         break;
 750
 751     next_loop:
 752       /* update the prec for next loop */
 753       MPFR_ZIV_NEXT (loop, prec);
 754     } /* loop */
 755   MPFR_ZIV_FREE (loop);
 756
 757   /* round y */
 758   if (mpfr_round_raw (MPFR_MANT (x), result,
 759                       ysize_bits,
 760                       pstr->negative, MPFR_PREC(x), rnd, &res ))
 761     {
 762       /* overflow when rounding y */
 763       MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT;
 764       /* Overflow Checking not needed */
 765       exp ++;
 766     }
 767
 768   if (res == 0) /* fix ternary value */
 769     {
 770       exact = exact && (pstr_size == pstr->prec);
 771       if (!exact)
 772         res = (pstr->negative) ? 1 : -1;
 773     }
 774
 775   /* Set sign of x before exp since check_range needs a valid sign */
 776   (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
 777
 778   /* DO NOT USE MPFR_SET_EXP. The exp may be out of range! */
 779   MPFR_SADD_OVERFLOW (exp, exp, ysize_bits,
 780                       mpfr_exp_t, mpfr_uexp_t,
 781                       MPFR_EXP_MIN, MPFR_EXP_MAX,
 782                       goto overflow, goto underflow);
 783   MPFR_EXP (x) = exp;
 784   res = mpfr_check_range (x, res, rnd);
 785   goto end;
 786
 787  underflow:
 788   /* This is called when there is a huge overflow
 789      (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 790   if (rnd == MPFR_RNDN)
 791     rnd = MPFR_RNDZ;
 792   res = mpfr_underflow (x, rnd, (pstr->negative) ? -1 : 1);
 793   goto end;
 794
 795  overflow:
 796   res = mpfr_overflow (x, rnd, (pstr->negative) ? -1 : 1);
 797
 798  end:
 799   MPFR_TMP_FREE (marker);
 800   return res;
 801 }
 802
 803 static void
 804 free_parsed_string (struct parsed_string *pstr)
 805 {
 806   (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
 807 }
 808
 809 int
 810 mpfr_strtofr (mpfr_t x, const char *string, char **end, int base,
 811               mpfr_rnd_t rnd)
 812 {
 813   int res;
 814   struct parsed_string pstr;
 815
 816   /* For base <= 36, parsing is case-insensitive. */
 817   MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 62));
 818
 819   /* If an error occured, it must return 0 */
 820   MPFR_SET_ZERO (x);
 821   MPFR_SET_POS (x);
 822
 823   MPFR_ASSERTN (MPFR_MAX_BASE >= 62);
 824   res = parse_string (x, &pstr, &string, base);
 825   /* If res == 0, then it was exact (NAN or INF),
 826      so it is also the ternary value */
 827   if (MPFR_UNLIKELY (res == -1))  /* invalid data */
 828     res = 0;  /* x is set to 0, which is exact, thus ternary value is 0 */
 829   else if (res == 1)
 830     {
 831       res = parsed_string_to_mpfr (x, &pstr, rnd);
 832       free_parsed_string (&pstr);
 833     }
 834   else if (res == 2)
 835     res = mpfr_overflow (x, rnd, (pstr.negative) ? -1 : 1);
 836   MPFR_ASSERTD (res != 3);
 837 #if 0
 838   else if (res == 3)
 839     {
 840       /* This is called when there is a huge overflow
 841          (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 842       if (rnd == MPFR_RNDN)
 843         rnd = MPFR_RNDZ;
 844       res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1);
 845     }
 846 #endif
 847
 848   if (end != NULL)
 849     *end = (char *) string;
 850   return res;
 851 }