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