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