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