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