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