contrib/mpfr/strtofr.c

   1 /* mpfr_strtofr -- set a floating-point number from a string
   2
   3 Copyright 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 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   mp_exp_t       exp_base; /* number of digits before the point */
  40   mp_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, GMP_RNDN);
 130       mpfr_log2 (x, x, GMP_RNDN);
 131       mpfr_ui_div (x, 1, x, GMP_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, GMP_RNDN);
 137           mpfr_sub (x, x, y, GMP_RNDN);
 138           mpfr_ui_div (x, 1, x, GMP_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) * sizeof(char) + 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[1];
 385       /* the exponent digits are kept in ASCII */
 386       mp_exp_t read_exp = strtol (str + 1, endptr, 10);
 387       mp_exp_t sum = 0;
 388       if (endptr[0] != str+1)
 389         str = endptr[0];
 390       MPFR_ASSERTN (read_exp == (long) read_exp);
 391       MPFR_SADD_OVERFLOW (sum, read_exp, pstr->exp_base,
 392                           mp_exp_t, mp_exp_unsigned_t,
 393                           MPFR_EXP_MIN, MPFR_EXP_MAX,
 394                           res = 2, res = 3);
 395       /* Since exp_base was positive, read_exp + exp_base can't
 396          do a negative overflow. */
 397       MPFR_ASSERTD (res != 3);
 398       pstr->exp_base = sum;
 399     }
 400   else if ((base == 2 || base == 16)
 401            && (*str == 'p' || *str == 'P')
 402            && (!isspace((unsigned char) str[1])))
 403     {
 404       char *endptr[1];
 405       pstr->exp_bin = (mp_exp_t) strtol (str + 1, endptr, 10);
 406       if (endptr[0] != str+1)
 407         str = endptr[0];
 408     }
 409
 410   /* Remove 0's at the beginning and end of mant_s[0..prec_s-1] */
 411   mant = pstr->mantissa;
 412   for ( ; (pstr->prec > 0) && (*mant == 0) ; mant++, pstr->prec--)
 413     pstr->exp_base--;
 414   for ( ; (pstr->prec > 0) && (mant[pstr->prec - 1] == 0); pstr->prec--);
 415   pstr->mant = mant;
 416
 417   /* Check if x = 0 */
 418   if (pstr->prec == 0)
 419     {
 420       MPFR_SET_ZERO (x);
 421       if (pstr->negative)
 422         MPFR_SET_NEG(x);
 423       else
 424         MPFR_SET_POS(x);
 425       res = 0;
 426     }
 427
 428   *string = str;
 429  end:
 430   if (pstr->mantissa != NULL && res != 1)
 431     (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
 432   return res;
 433 }
 434
 435 /* Transform a parsed string to a mpfr_t according to the rounding mode
 436    and the precision of x.
 437    Returns the ternary value. */
 438 static int
 439 parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd)
 440 {
 441   mp_prec_t prec;
 442   mp_exp_t  exp;
 443   mp_exp_t  ysize_bits;
 444   mp_limb_t *y, *result;
 445   int count, exact;
 446   size_t pstr_size;
 447   mp_size_t ysize, real_ysize;
 448   int res, err;
 449   MPFR_ZIV_DECL (loop);
 450   MPFR_TMP_DECL (marker);
 451
 452   /* initialize the working precision */
 453   prec = MPFR_PREC (x) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (x));
 454
 455   /* compute y as long as rounding is not possible */
 456   MPFR_TMP_MARK(marker);
 457   MPFR_ZIV_INIT (loop, prec);
 458   for (;;)
 459     {
 460       /* Set y to the value of the ~prec most significant bits of pstr->mant
 461          (as long as we guarantee correct rounding, we don't need to get
 462          exactly prec bits). */
 463       ysize = (prec - 1) / BITS_PER_MP_LIMB + 1;
 464       /* prec bits corresponds to ysize limbs */
 465       ysize_bits = ysize * BITS_PER_MP_LIMB;
 466       /* and to ysize_bits >= prec > MPFR_PREC (x) bits */
 467       y = (mp_limb_t*) MPFR_TMP_ALLOC ((2 * ysize + 1) * sizeof (mp_limb_t));
 468       y += ysize; /* y has (ysize+1) allocated limbs */
 469
 470       /* pstr_size is the number of characters we read in pstr->mant
 471          to have at least ysize full limbs.
 472          We must have base^(pstr_size-1) >= (2^(BITS_PER_MP_LIMB))^ysize
 473          (in the worst case, the first digit is one and all others are zero).
 474          i.e., pstr_size >= 1 + ysize*BITS_PER_MP_LIMB/log2(base)
 475           Since ysize ~ prec/BITS_PER_MP_LIMB and prec < Umax/2 =>
 476           ysize*BITS_PER_MP_LIMB can not overflow.
 477          We compute pstr_size = 1 + ceil(ysize_bits * Num / Den)
 478           where Num/Den >= 1/log2(base)
 479          It is not exactly ceil(1/log2(base)) but could be one more (base 2)
 480          Quite ugly since it tries to avoid overflow:
 481          let Num = RedInvLog2Table[pstr->base-2][0]
 482          and Den = RedInvLog2Table[pstr->base-2][1],
 483          and ysize_bits = a*Den+b,
 484          then ysize_bits * Num/Den = a*Num + (b * Num)/Den,
 485          thus ceil(ysize_bits * Num/Den) = a*Num + floor(b * Num + Den - 1)/Den
 486       */
 487       {
 488         unsigned long Num = RedInvLog2Table[pstr->base-2][0];
 489         unsigned long Den = RedInvLog2Table[pstr->base-2][1];
 490         pstr_size = ((ysize_bits / Den) * Num)
 491           + (((ysize_bits % Den) * Num + Den - 1) / Den)
 492           + 1;
 493       }
 494
 495       /* since pstr_size corresponds to at least ysize_bits full bits,
 496          and ysize_bits > prec, the weight of the neglected part of
 497          pstr->mant (if any) is < ulp(y) < ulp(x) */
 498
 499       /* if the number of wanted characters is more than what we have in
 500          pstr->mant, round it down */
 501       if (pstr_size >= pstr->prec)
 502         pstr_size = pstr->prec;
 503       MPFR_ASSERTD (pstr_size == (mp_exp_t) pstr_size);
 504
 505       /* convert str into binary: note that pstr->mant is big endian,
 506          thus no offset is needed */
 507       real_ysize = mpn_set_str (y, pstr->mant, pstr_size, pstr->base);
 508       MPFR_ASSERTD (real_ysize <= ysize+1);
 509
 510       /* normalize y: warning we can get even get ysize+1 limbs! */
 511       MPFR_ASSERTD (y[real_ysize - 1] != 0); /* mpn_set_str guarantees this */
 512       count_leading_zeros (count, y[real_ysize - 1]);
 513       /* exact means that the number of limbs of the output of mpn_set_str
 514          is less or equal to ysize */
 515       exact = real_ysize <= ysize;
 516       if (exact) /* shift y to the left in that case y should be exact */
 517         {
 518           /* we have enough limbs to store {y, real_ysize} */
 519           /* shift {y, num_limb} for count bits to the left */
 520           if (count != 0)
 521             mpn_lshift (y + ysize - real_ysize, y, real_ysize, count);
 522           if (real_ysize != ysize)
 523             {
 524               if (count == 0)
 525                 MPN_COPY_DECR (y + ysize - real_ysize, y, real_ysize);
 526               MPN_ZERO (y, ysize - real_ysize);
 527             }
 528           /* for each bit shift decrease exponent of y */
 529           /* (This should not overflow) */
 530           exp = - ((ysize - real_ysize) * BITS_PER_MP_LIMB + count);
 531         }
 532       else  /* shift y to the right, by doing this we might lose some
 533                bits from the result of mpn_set_str (in addition to the
 534                characters neglected from pstr->mant) */
 535         {
 536           /* shift {y, num_limb} for (BITS_PER_MP_LIMB - count) bits
 537              to the right. FIXME: can we prove that count cannot be zero here,
 538              since mpn_rshift does not accept a shift of BITS_PER_MP_LIMB? */
 539           MPFR_ASSERTD (count != 0);
 540           exact = mpn_rshift (y, y, real_ysize, BITS_PER_MP_LIMB - count) ==
 541             MPFR_LIMB_ZERO;
 542           /* for each bit shift increase exponent of y */
 543           exp = BITS_PER_MP_LIMB - count;
 544         }
 545
 546       /* compute base^(exp_s-pr) on n limbs */
 547       if (IS_POW2 (pstr->base))
 548         {
 549           /* Base: 2, 4, 8, 16, 32 */
 550           int pow2;
 551           mp_exp_t tmp;
 552
 553           count_leading_zeros (pow2, (mp_limb_t) pstr->base);
 554           pow2 = BITS_PER_MP_LIMB - pow2 - 1; /* base = 2^pow2 */
 555           MPFR_ASSERTD (0 < pow2 && pow2 <= 5);
 556           /* exp += pow2 * (pstr->exp_base - pstr_size) + pstr->exp_bin
 557              with overflow checking
 558              and check that we can add/substract 2 to exp without overflow */
 559           MPFR_SADD_OVERFLOW (tmp, pstr->exp_base, -(mp_exp_t) pstr_size,
 560                               mp_exp_t, mp_exp_unsigned_t,
 561                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 562                               goto overflow, goto underflow);
 563           /* On some FreeBsd/Alpha, LONG_MIN/1 produces an exception
 564              so we check for this before doing the division */
 565           if (tmp > 0 && pow2 != 1 && MPFR_EXP_MAX/pow2 <= tmp)
 566             goto overflow;
 567           else if (tmp < 0 && pow2 != 1 && MPFR_EXP_MIN/pow2 >= tmp)
 568             goto underflow;
 569           tmp *= pow2;
 570           MPFR_SADD_OVERFLOW (tmp, tmp, pstr->exp_bin,
 571                               mp_exp_t, mp_exp_unsigned_t,
 572                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 573                               goto overflow, goto underflow);
 574           MPFR_SADD_OVERFLOW (exp, exp, tmp,
 575                               mp_exp_t, mp_exp_unsigned_t,
 576                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 577                               goto overflow, goto underflow);
 578           result = y;
 579           err = 0;
 580         }
 581       /* case non-power-of-two-base, and pstr->exp_base > pstr_size */
 582       else if (pstr->exp_base > (mp_exp_t) pstr_size)
 583         {
 584           mp_limb_t *z;
 585           mp_exp_t   exp_z;
 586
 587           result = (mp_limb_t*) MPFR_TMP_ALLOC ((2*ysize+1)*BYTES_PER_MP_LIMB);
 588
 589           /* z = base^(exp_base-sptr_size) using space allocated at y-ysize */
 590           z = y - ysize;
 591           /* NOTE: exp_base-pstr_size can't overflow since pstr_size > 0 */
 592           err = mpfr_mpn_exp (z, &exp_z, pstr->base,
 593                               pstr->exp_base - pstr_size, ysize);
 594           if (err == -2)
 595             goto overflow;
 596           exact = exact && (err == -1);
 597
 598           /* If exact is non zero, then z equals exactly the value of the
 599              pstr_size most significant digits from pstr->mant, i.e., the
 600              only difference can come from the neglected pstr->prec-pstr_size
 601              least significant digits of pstr->mant.
 602              If exact is zero, then z is rounded toward zero with respect
 603              to that value. */
 604
 605           /* multiply(y = 0.mant_s[0]...mant_s[pr-1])_base by base^(exp_s-g) */
 606           mpn_mul_n (result, y, z, ysize);
 607
 608           /* compute the error on the product */
 609           if (err == -1)
 610             err = 0;
 611           err ++;
 612
 613           /* compute the exponent of y */
 614           /* exp += exp_z + ysize_bits with overflow checking
 615              and check that we can add/substract 2 to exp without overflow */
 616           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 617                               mp_exp_t, mp_exp_unsigned_t,
 618                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 619                               goto overflow, goto underflow);
 620           MPFR_SADD_OVERFLOW (exp, exp, exp_z,
 621                               mp_exp_t, mp_exp_unsigned_t,
 622                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 623                               goto overflow, goto underflow);
 624
 625           /* normalize result */
 626           if (MPFR_LIMB_MSB (result[2 * ysize - 1]) == 0)
 627             {
 628               mp_limb_t *r = result + ysize - 1;
 629               mpn_lshift (r, r, ysize + 1, 1);
 630               /* Overflow checking not needed */
 631               exp --;
 632             }
 633
 634           /* if the low ysize limbs of {result, 2*ysize} are all zero,
 635              then the result is still "exact" (if it was before) */
 636           exact = exact && (mpn_scan1 (result, 0)
 637                             >= (unsigned long) ysize_bits);
 638           result += ysize;
 639         }
 640       /* case exp_base < pstr_size */
 641       else if (pstr->exp_base < (mp_exp_t) pstr_size)
 642         {
 643           mp_limb_t *z;
 644           mp_exp_t exp_z;
 645
 646           result = (mp_limb_t*) MPFR_TMP_ALLOC ((3*ysize+1) * BYTES_PER_MP_LIMB);
 647
 648           /* set y to y * K^ysize */
 649           y = y - ysize;  /* we have allocated ysize limbs at y - ysize */
 650           MPN_ZERO (y, ysize);
 651
 652           /* pstr_size - pstr->exp_base can overflow */
 653           MPFR_SADD_OVERFLOW (exp_z, (mp_exp_t) pstr_size, -pstr->exp_base,
 654                               mp_exp_t, mp_exp_unsigned_t,
 655                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 656                               goto underflow, goto overflow);
 657
 658           /* (z, exp_z) = base^(exp_base-pstr_size) */
 659           z = result + 2*ysize + 1;
 660           err = mpfr_mpn_exp (z, &exp_z, pstr->base, exp_z, ysize);
 661           exact = exact && (err == -1);
 662           if (err == -2)
 663             goto underflow; /* FIXME: Sure? */
 664           if (err == -1)
 665             err = 0;
 666
 667           /* compute y / z */
 668           /* result will be put into result + n, and remainder into result */
 669           mpn_tdiv_qr (result + ysize, result, (mp_size_t) 0, y,
 670                        2 * ysize, z, ysize);
 671
 672           /* exp -= exp_z + ysize_bits with overflow checking
 673              and check that we can add/substract 2 to exp without overflow */
 674           MPFR_SADD_OVERFLOW (exp_z, exp_z, ysize_bits,
 675                               mp_exp_t, mp_exp_unsigned_t,
 676                               MPFR_EXP_MIN, MPFR_EXP_MAX,
 677                               goto underflow, goto overflow);
 678           MPFR_SADD_OVERFLOW (exp, exp, -exp_z,
 679                               mp_exp_t, mp_exp_unsigned_t,
 680                               MPFR_EXP_MIN+2, MPFR_EXP_MAX-2,
 681                               goto overflow, goto underflow);
 682           err += 2;
 683           /* if the remainder of the division is zero, then the result is
 684              still "exact" if it was before */
 685           exact = exact && (mpn_popcount (result, ysize) == 0);
 686
 687           /* normalize result */
 688           if (result[2 * ysize] == MPFR_LIMB_ONE)
 689             {
 690               mp_limb_t *r = result + ysize;
 691               exact = exact && ((*r & MPFR_LIMB_ONE) == 0);
 692               mpn_rshift (r, r, ysize + 1, 1);
 693               /* Overflow Checking not needed */
 694               exp ++;
 695             }
 696           result += ysize;
 697         }
 698       /* case exp_base = pstr_size: no multiplication or division needed */
 699       else
 700         {
 701           /* base^(exp_s-pr) = 1             nothing to compute */
 702           result = y;
 703           err = 0;
 704         }
 705
 706       /* If result is exact, we still have to consider the neglected part
 707          of the input string. For a directed rounding, in that case we could
 708          still correctly round, since the neglected part is less than
 709          one ulp, but that would make the code more complex, and give a
 710          speedup for rare cases only. */
 711       exact = exact && (pstr_size == pstr->prec);
 712
 713       /* at this point, result is an approximation rounded toward zero
 714          of the pstr_size most significant digits of pstr->mant, with
 715          equality in case exact is non-zero. */
 716
 717       /* test if rounding is possible, and if so exit the loop */
 718       if (exact || mpfr_can_round_raw (result, ysize,
 719                                        (pstr->negative) ? -1 : 1,
 720                                        ysize_bits - err - 1,
 721                                        GMP_RNDN, rnd, MPFR_PREC(x)))
 722         break;
 723
 724       /* update the prec for next loop */
 725       MPFR_ZIV_NEXT (loop, prec);
 726     } /* loop */
 727   MPFR_ZIV_FREE (loop);
 728
 729   /* round y */
 730   if (mpfr_round_raw (MPFR_MANT (x), result,
 731                       ysize_bits,
 732                       pstr->negative, MPFR_PREC(x), rnd, &res ))
 733     {
 734       /* overflow when rounding y */
 735       MPFR_MANT (x)[MPFR_LIMB_SIZE (x) - 1] = MPFR_LIMB_HIGHBIT;
 736       /* Overflow Checking not needed */
 737       exp ++;
 738     }
 739
 740   if (res == 0) /* fix ternary value */
 741     {
 742       exact = exact && (pstr_size == pstr->prec);
 743       if (!exact)
 744         res = (pstr->negative) ? 1 : -1;
 745     }
 746
 747   /* Set sign of x before exp since check_range needs a valid sign */
 748   (pstr->negative) ? MPFR_SET_NEG (x) : MPFR_SET_POS (x);
 749
 750   /* DO NOT USE MPFR_SET_EXP. The exp may be out of range! */
 751   MPFR_SADD_OVERFLOW (exp, exp, ysize_bits,
 752                       mp_exp_t, mp_exp_unsigned_t,
 753                       MPFR_EXP_MIN, MPFR_EXP_MAX,
 754                       goto overflow, goto underflow);
 755   MPFR_EXP (x) = exp;
 756   res = mpfr_check_range (x, res, rnd);
 757   goto end;
 758
 759  underflow:
 760   /* This is called when there is a huge overflow
 761      (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 762   if (rnd == GMP_RNDN)
 763     rnd = GMP_RNDZ;
 764   res = mpfr_underflow (x, rnd, (pstr->negative) ? -1 : 1);
 765   goto end;
 766
 767  overflow:
 768   res = mpfr_overflow (x, rnd, (pstr->negative) ? -1 : 1);
 769
 770  end:
 771   MPFR_TMP_FREE (marker);
 772   return res;
 773 }
 774
 775 static void
 776 free_parsed_string (struct parsed_string *pstr)
 777 {
 778   (*__gmp_free_func) (pstr->mantissa, pstr->alloc);
 779 }
 780
 781 int
 782 mpfr_strtofr (mpfr_t x, const char *string, char **end, int base,
 783               mp_rnd_t rnd)
 784 {
 785   int res;
 786   struct parsed_string pstr;
 787
 788   MPFR_ASSERTN (base == 0 || (base >= 2 && base <= 36));
 789
 790   /* If an error occured, it must return 0 */
 791   MPFR_SET_ZERO (x);
 792   MPFR_SET_POS (x);
 793
 794   /* Though bases up to MPFR_MAX_BASE are supported, we require a lower
 795      limit: 36. For such values <= 36, parsing is case-insensitive. */
 796   MPFR_ASSERTN (MPFR_MAX_BASE >= 36);
 797   res = parse_string (x, &pstr, &string, base);
 798   /* If res == 0, then it was exact (NAN or INF),
 799      so it is also the ternary value */
 800   if (MPFR_UNLIKELY (res == -1))  /* invalid data */
 801     res = 0;  /* x is set to 0, which is exact, thus ternary value is 0 */
 802   else if (res == 1)
 803     {
 804       res = parsed_string_to_mpfr (x, &pstr, rnd);
 805       free_parsed_string (&pstr);
 806     }
 807   else if (res == 2)
 808     res = mpfr_overflow (x, rnd, (pstr.negative) ? -1 : 1);
 809   MPFR_ASSERTD (res != 3);
 810 #if 0
 811   else if (res == 3)
 812     {
 813       /* This is called when there is a huge overflow
 814          (Real expo < MPFR_EXP_MIN << __gmpfr_emin */
 815       if (rnd == GMP_RNDN)
 816         rnd = GMP_RNDZ;
 817       res = mpfr_underflow (x, rnd, (pstr.negative) ? -1 : 1);
 818     }
 819 #endif
 820
 821   if (end != NULL)
 822     *end = (char *) string;
 823   return res;
 824 }