1 /* Convert string representing a number to float value, using given locale. 2 Copyright (C) 1997-2012 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. 5 6 The GNU C Library is free software; you can redistribute it and/or 7 modify it under the terms of the GNU Lesser General Public 8 License as published by the Free Software Foundation; either 9 version 2.1 of the License, or (at your option) any later version. 10 11 The GNU C Library is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 Lesser General Public License for more details. 15 16 You should have received a copy of the GNU Lesser General Public 17 License along with the GNU C Library; if not, see 18 <http://www.gnu.org/licenses/>. */ 19 20 #include <config.h> 21 #include <stdarg.h> 22 #include <string.h> 23 #include <stdint.h> 24 #include <stdbool.h> 25 #include <float.h> 26 #include <math.h> 27 #define NDEBUG 1 28 #include <assert.h> 29 #ifdef HAVE_ERRNO_H 30 #include <errno.h> 31 #endif 32 33 #ifdef HAVE_FENV_H 34 #include <fenv.h> 35 #endif 36 37 #ifdef HAVE_FENV_H 38 #include "quadmath-rounding-mode.h" 39 #endif 40 #include "../printf/quadmath-printf.h" 41 #include "../printf/fpioconst.h" 42 43 #undef L_ 44 #ifdef USE_WIDE_CHAR 45 # define STRING_TYPE wchar_t 46 # define CHAR_TYPE wint_t 47 # define L_(Ch) L##Ch 48 # define ISSPACE(Ch) __iswspace_l ((Ch), loc) 49 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc) 50 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc) 51 # define TOLOWER(Ch) __towlower_l ((Ch), loc) 52 # define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr) 53 # define STRNCASECMP(S1, S2, N) \ 54 __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr) 55 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc) 56 #else 57 # define STRING_TYPE char 58 # define CHAR_TYPE char 59 # define L_(Ch) Ch 60 # define ISSPACE(Ch) isspace (Ch) 61 # define ISDIGIT(Ch) isdigit (Ch) 62 # define ISXDIGIT(Ch) isxdigit (Ch) 63 # define TOLOWER(Ch) tolower (Ch) 64 # define TOLOWER_C(Ch) \ 65 ({__typeof(Ch) __tlc = (Ch); \ 66 (__tlc >= 'A' && __tlc <= 'Z') ? __tlc - 'A' + 'a' : __tlc; }) 67 # define STRNCASECMP(S1, S2, N) \ 68 __quadmath_strncasecmp_c (S1, S2, N) 69 # ifdef HAVE_STRTOULL 70 # define STRTOULL(S, E, B) strtoull (S, E, B) 71 # else 72 # define STRTOULL(S, E, B) strtoul (S, E, B) 73 # endif 74 75 static inline int 76 __quadmath_strncasecmp_c (const char *s1, const char *s2, size_t n) 77 { 78 const unsigned char *p1 = (const unsigned char *) s1; 79 const unsigned char *p2 = (const unsigned char *) s2; 80 int result; 81 if (p1 == p2 || n == 0) 82 return 0; 83 while ((result = TOLOWER_C (*p1) - TOLOWER_C (*p2++)) == 0) 84 if (*p1++ == '\0' || --n == 0) 85 break; 86 87 return result; 88 } 89 #endif 90 91 92 /* Constants we need from float.h; select the set for the FLOAT precision. */ 93 #define MANT_DIG PASTE(FLT,_MANT_DIG) 94 #define DIG PASTE(FLT,_DIG) 95 #define MAX_EXP PASTE(FLT,_MAX_EXP) 96 #define MIN_EXP PASTE(FLT,_MIN_EXP) 97 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP) 98 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP) 99 #define MAX_VALUE PASTE(FLT,_MAX) 100 #define MIN_VALUE PASTE(FLT,_MIN) 101 102 /* Extra macros required to get FLT expanded before the pasting. */ 103 #define PASTE(a,b) PASTE1(a,b) 104 #define PASTE1(a,b) a##b 105 106 /* Function to construct a floating point number from an MP integer 107 containing the fraction bits, a base 2 exponent, and a sign flag. */ 108 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative); 109 110 /* Definitions according to limb size used. */ 111 #if BITS_PER_MP_LIMB == 32 112 # define MAX_DIG_PER_LIMB 9 113 # define MAX_FAC_PER_LIMB 1000000000UL 114 #elif BITS_PER_MP_LIMB == 64 115 # define MAX_DIG_PER_LIMB 19 116 # define MAX_FAC_PER_LIMB 10000000000000000000ULL 117 #else 118 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for" 119 #endif 120 121 #define _tens_in_limb __quadmath_tens_in_limb 122 extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] attribute_hidden; 123 124 #ifndef howmany 125 #define howmany(x,y) (((x)+((y)-1))/(y)) 126 #endif 127 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; }) 128 129 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG) 130 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG) 131 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB) 132 133 #define RETURN(val,end) \ 134 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \ 135 return val; } while (0) 136 137 /* Maximum size necessary for mpn integers to hold floating point 138 numbers. The largest number we need to hold is 10^n where 2^-n is 139 1/4 ulp of the smallest representable value (that is, n = MANT_DIG 140 - MIN_EXP + 2). Approximate using 10^3 < 2^10. */ 141 #define MPNSIZE (howmany (1 + ((MANT_DIG - MIN_EXP + 2) * 10) / 3, \ 142 BITS_PER_MP_LIMB) + 2) 143 /* Declare an mpn integer variable that big. */ 144 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size 145 /* Copy an mpn integer value. */ 146 #define MPN_ASSIGN(dst, src) \ 147 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t)) 148 149 /* Set errno and return an overflowing value with sign specified by 150 NEGATIVE. */ 151 static FLOAT 152 overflow_value (int negative) 153 { 154 #if defined HAVE_ERRNO_H && defined ERANGE 155 errno = ERANGE; 156 #endif 157 FLOAT result = (negative ? -MAX_VALUE : MAX_VALUE) * MAX_VALUE; 158 return result; 159 } 160 161 /* Set errno and return an underflowing value with sign specified by 162 NEGATIVE. */ 163 static FLOAT 164 underflow_value (int negative) 165 { 166 #if defined HAVE_ERRNO_H && defined ERANGE 167 errno = ERANGE; 168 #endif 169 FLOAT result = (negative ? -MIN_VALUE : MIN_VALUE) * MIN_VALUE; 170 return result; 171 } 172 173 /* Return a floating point number of the needed type according to the given 174 multi-precision number after possible rounding. */ 175 static FLOAT 176 round_and_return (mp_limb_t *retval, intmax_t exponent, int negative, 177 mp_limb_t round_limb, mp_size_t round_bit, int more_bits) 178 { 179 #ifdef HAVE_FENV_H 180 int mode = get_rounding_mode (); 181 #endif 182 183 if (exponent < MIN_EXP - 1) 184 { 185 mp_size_t shift; 186 bool is_tiny; 187 188 if (exponent < MIN_EXP - 1 - MANT_DIG) 189 return underflow_value (negative); 190 191 shift = MIN_EXP - 1 - exponent; 192 is_tiny = true; 193 194 more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0; 195 if (shift == MANT_DIG) 196 /* This is a special case to handle the very seldom case where 197 the mantissa will be empty after the shift. */ 198 { 199 int i; 200 201 round_limb = retval[RETURN_LIMB_SIZE - 1]; 202 round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; 203 for (i = 0; i < RETURN_LIMB_SIZE; ++i) 204 more_bits |= retval[i] != 0; 205 MPN_ZERO (retval, RETURN_LIMB_SIZE); 206 } 207 else if (shift >= BITS_PER_MP_LIMB) 208 { 209 int i; 210 211 round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB]; 212 round_bit = (shift - 1) % BITS_PER_MP_LIMB; 213 for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i) 214 more_bits |= retval[i] != 0; 215 more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) 216 != 0); 217 218 (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB], 219 RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB), 220 shift % BITS_PER_MP_LIMB); 221 MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)], 222 shift / BITS_PER_MP_LIMB); 223 } 224 else if (shift > 0) 225 { 226 #ifdef HAVE_FENV_H 227 if (TININESS_AFTER_ROUNDING && shift == 1) 228 { 229 /* Whether the result counts as tiny depends on whether, 230 after rounding to the normal precision, it still has 231 a subnormal exponent. */ 232 mp_limb_t retval_normal[RETURN_LIMB_SIZE]; 233 if (round_away (negative, 234 (retval[0] & 1) != 0, 235 (round_limb 236 & (((mp_limb_t) 1) << round_bit)) != 0, 237 (more_bits 238 || ((round_limb 239 & ((((mp_limb_t) 1) << round_bit) - 1)) 240 != 0)), 241 mode)) 242 { 243 mp_limb_t cy = mpn_add_1 (retval_normal, retval, 244 RETURN_LIMB_SIZE, 1); 245 246 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || 247 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && 248 ((retval_normal[RETURN_LIMB_SIZE - 1] 249 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) 250 != 0))) 251 is_tiny = false; 252 } 253 } 254 #endif 255 round_limb = retval[0]; 256 round_bit = shift - 1; 257 (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift); 258 } 259 /* This is a hook for the m68k long double format, where the 260 exponent bias is the same for normalized and denormalized 261 numbers. */ 262 #ifndef DENORM_EXP 263 # define DENORM_EXP (MIN_EXP - 2) 264 #endif 265 exponent = DENORM_EXP; 266 if (is_tiny 267 && ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0 268 || more_bits 269 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0)) 270 { 271 #if defined HAVE_ERRNO_H && defined ERANGE 272 errno = ERANGE; 273 #endif 274 volatile FLOAT force_underflow_exception = MIN_VALUE * MIN_VALUE; 275 (void) force_underflow_exception; 276 } 277 } 278 279 if (exponent > MAX_EXP) 280 goto overflow; 281 282 #ifdef HAVE_FENV_H 283 if (round_away (negative, 284 (retval[0] & 1) != 0, 285 (round_limb & (((mp_limb_t) 1) << round_bit)) != 0, 286 (more_bits 287 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0), 288 mode)) 289 { 290 mp_limb_t cy = mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1); 291 292 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || 293 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && 294 (retval[RETURN_LIMB_SIZE - 1] 295 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)) 296 { 297 ++exponent; 298 (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1); 299 retval[RETURN_LIMB_SIZE - 1] 300 |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB); 301 } 302 else if (exponent == DENORM_EXP 303 && (retval[RETURN_LIMB_SIZE - 1] 304 & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) 305 != 0) 306 /* The number was denormalized but now normalized. */ 307 exponent = MIN_EXP - 1; 308 } 309 #endif 310 311 if (exponent > MAX_EXP) 312 overflow: 313 return overflow_value (negative); 314 315 return MPN2FLOAT (retval, exponent, negative); 316 } 317 318 319 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits 320 into N. Return the size of the number limbs in NSIZE at the first 321 character od the string that is not part of the integer as the function 322 value. If the EXPONENT is small enough to be taken as an additional 323 factor for the resulting number (see code) multiply by it. */ 324 static const STRING_TYPE * 325 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize, 326 intmax_t *exponent 327 #ifndef USE_WIDE_CHAR 328 , const char *decimal, size_t decimal_len, const char *thousands 329 #endif 330 331 ) 332 { 333 /* Number of digits for actual limb. */ 334 int cnt = 0; 335 mp_limb_t low = 0; 336 mp_limb_t start; 337 338 *nsize = 0; 339 assert (digcnt > 0); 340 do 341 { 342 if (cnt == MAX_DIG_PER_LIMB) 343 { 344 if (*nsize == 0) 345 { 346 n[0] = low; 347 *nsize = 1; 348 } 349 else 350 { 351 mp_limb_t cy; 352 cy = mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB); 353 cy += mpn_add_1 (n, n, *nsize, low); 354 if (cy != 0) 355 { 356 assert (*nsize < MPNSIZE); 357 n[*nsize] = cy; 358 ++(*nsize); 359 } 360 } 361 cnt = 0; 362 low = 0; 363 } 364 365 /* There might be thousands separators or radix characters in 366 the string. But these all can be ignored because we know the 367 format of the number is correct and we have an exact number 368 of characters to read. */ 369 #ifdef USE_WIDE_CHAR 370 if (*str < L'0' || *str > L'9') 371 ++str; 372 #else 373 if (*str < '0' || *str > '9') 374 { 375 int inner = 0; 376 if (thousands != NULL && *str == *thousands 377 && ({ for (inner = 1; thousands[inner] != '\0'; ++inner) 378 if (thousands[inner] != str[inner]) 379 break; 380 thousands[inner] == '\0'; })) 381 str += inner; 382 else 383 str += decimal_len; 384 } 385 #endif 386 low = low * 10 + *str++ - L_('0'); 387 ++cnt; 388 } 389 while (--digcnt > 0); 390 391 if (*exponent > 0 && *exponent <= MAX_DIG_PER_LIMB - cnt) 392 { 393 low *= _tens_in_limb[*exponent]; 394 start = _tens_in_limb[cnt + *exponent]; 395 *exponent = 0; 396 } 397 else 398 start = _tens_in_limb[cnt]; 399 400 if (*nsize == 0) 401 { 402 n[0] = low; 403 *nsize = 1; 404 } 405 else 406 { 407 mp_limb_t cy; 408 cy = mpn_mul_1 (n, n, *nsize, start); 409 cy += mpn_add_1 (n, n, *nsize, low); 410 if (cy != 0) 411 { 412 assert (*nsize < MPNSIZE); 413 n[(*nsize)++] = cy; 414 } 415 } 416 417 return str; 418 } 419 420 421 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits 422 with the COUNT most significant bits of LIMB. 423 424 Implemented as a macro, so that __builtin_constant_p works even at -O0. 425 426 Tege doesn't like this macro so I have to write it here myself. :) 427 --drepper */ 428 #define mpn_lshift_1(ptr, size, count, limb) \ 429 do \ 430 { \ 431 mp_limb_t *__ptr = (ptr); \ 432 if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \ 433 { \ 434 mp_size_t i; \ 435 for (i = (size) - 1; i > 0; --i) \ 436 __ptr[i] = __ptr[i - 1]; \ 437 __ptr[0] = (limb); \ 438 } \ 439 else \ 440 { \ 441 /* We assume count > 0 && count < BITS_PER_MP_LIMB here. */ \ 442 unsigned int __count = (count); \ 443 (void) mpn_lshift (__ptr, __ptr, size, __count); \ 444 __ptr[0] |= (limb) >> (BITS_PER_MP_LIMB - __count); \ 445 } \ 446 } \ 447 while (0) 448 449 450 #define INTERNAL(x) INTERNAL1(x) 451 #define INTERNAL1(x) __##x##_internal 452 #ifndef ____STRTOF_INTERNAL 453 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF) 454 #endif 455 456 /* This file defines a function to check for correct grouping. */ 457 #include "grouping.h" 458 459 460 /* Return a floating point number with the value of the given string NPTR. 461 Set *ENDPTR to the character after the last used one. If the number is 462 smaller than the smallest representable number, set `errno' to ERANGE and 463 return 0.0. If the number is too big to be represented, set `errno' to 464 ERANGE and return HUGE_VAL with the appropriate sign. */ 465 FLOAT 466 ____STRTOF_INTERNAL (nptr, endptr, group) 467 const STRING_TYPE *nptr; 468 STRING_TYPE **endptr; 469 int group; 470 { 471 int negative; /* The sign of the number. */ 472 MPN_VAR (num); /* MP representation of the number. */ 473 intmax_t exponent; /* Exponent of the number. */ 474 475 /* Numbers starting `0X' or `0x' have to be processed with base 16. */ 476 int base = 10; 477 478 /* When we have to compute fractional digits we form a fraction with a 479 second multi-precision number (and we sometimes need a second for 480 temporary results). */ 481 MPN_VAR (den); 482 483 /* Representation for the return value. */ 484 mp_limb_t retval[RETURN_LIMB_SIZE]; 485 /* Number of bits currently in result value. */ 486 int bits; 487 488 /* Running pointer after the last character processed in the string. */ 489 const STRING_TYPE *cp, *tp; 490 /* Start of significant part of the number. */ 491 const STRING_TYPE *startp, *start_of_digits; 492 /* Points at the character following the integer and fractional digits. */ 493 const STRING_TYPE *expp; 494 /* Total number of digit and number of digits in integer part. */ 495 size_t dig_no, int_no, lead_zero; 496 /* Contains the last character read. */ 497 CHAR_TYPE c; 498 499 /* We should get wint_t from <stddef.h>, but not all GCC versions define it 500 there. So define it ourselves if it remains undefined. */ 501 #ifndef _WINT_T 502 typedef unsigned int wint_t; 503 #endif 504 /* The radix character of the current locale. */ 505 #ifdef USE_WIDE_CHAR 506 wchar_t decimal; 507 #else 508 const char *decimal; 509 size_t decimal_len; 510 #endif 511 /* The thousands character of the current locale. */ 512 #ifdef USE_WIDE_CHAR 513 wchar_t thousands = L'\0'; 514 #else 515 const char *thousands = NULL; 516 #endif 517 /* The numeric grouping specification of the current locale, 518 in the format described in <locale.h>. */ 519 const char *grouping; 520 /* Used in several places. */ 521 int cnt; 522 523 #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO 524 const struct lconv *lc = localeconv (); 525 #endif 526 527 if (__builtin_expect (group, 0)) 528 { 529 #ifdef USE_NL_LANGINFO 530 grouping = nl_langinfo (GROUPING); 531 if (*grouping <= 0 || *grouping == CHAR_MAX) 532 grouping = NULL; 533 else 534 { 535 /* Figure out the thousands separator character. */ 536 #ifdef USE_WIDE_CHAR 537 thousands = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC); 538 if (thousands == L'\0') 539 grouping = NULL; 540 #else 541 thousands = nl_langinfo (THOUSANDS_SEP); 542 if (*thousands == '\0') 543 { 544 thousands = NULL; 545 grouping = NULL; 546 } 547 #endif 548 } 549 #elif defined USE_LOCALECONV 550 grouping = lc->grouping; 551 if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX) 552 grouping = NULL; 553 else 554 { 555 /* Figure out the thousands separator character. */ 556 thousands = lc->thousands_sep; 557 if (thousands == NULL || *thousands == '\0') 558 { 559 thousands = NULL; 560 grouping = NULL; 561 } 562 } 563 #else 564 grouping = NULL; 565 #endif 566 } 567 else 568 grouping = NULL; 569 570 /* Find the locale's decimal point character. */ 571 #ifdef USE_WIDE_CHAR 572 decimal = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); 573 assert (decimal != L'\0'); 574 # define decimal_len 1 575 #else 576 #ifdef USE_NL_LANGINFO 577 decimal = nl_langinfo (DECIMAL_POINT); 578 decimal_len = strlen (decimal); 579 assert (decimal_len > 0); 580 #elif defined USE_LOCALECONV 581 decimal = lc->decimal_point; 582 if (decimal == NULL || *decimal == '\0') 583 decimal = "."; 584 decimal_len = strlen (decimal); 585 #else 586 decimal = "."; 587 decimal_len = 1; 588 #endif 589 #endif 590 591 /* Prepare number representation. */ 592 exponent = 0; 593 negative = 0; 594 bits = 0; 595 596 /* Parse string to get maximal legal prefix. We need the number of 597 characters of the integer part, the fractional part and the exponent. */ 598 cp = nptr - 1; 599 /* Ignore leading white space. */ 600 do 601 c = *++cp; 602 while (ISSPACE (c)); 603 604 /* Get sign of the result. */ 605 if (c == L_('-')) 606 { 607 negative = 1; 608 c = *++cp; 609 } 610 else if (c == L_('+')) 611 c = *++cp; 612 613 /* Return 0.0 if no legal string is found. 614 No character is used even if a sign was found. */ 615 #ifdef USE_WIDE_CHAR 616 if (c == (wint_t) decimal 617 && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9') 618 { 619 /* We accept it. This funny construct is here only to indent 620 the code correctly. */ 621 } 622 #else 623 for (cnt = 0; decimal[cnt] != '\0'; ++cnt) 624 if (cp[cnt] != decimal[cnt]) 625 break; 626 if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9') 627 { 628 /* We accept it. This funny construct is here only to indent 629 the code correctly. */ 630 } 631 #endif 632 else if (c < L_('0') || c > L_('9')) 633 { 634 /* Check for `INF' or `INFINITY'. */ 635 CHAR_TYPE lowc = TOLOWER_C (c); 636 637 if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0) 638 { 639 /* Return +/- infinity. */ 640 if (endptr != NULL) 641 *endptr = (STRING_TYPE *) 642 (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0 643 ? 8 : 3)); 644 645 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; 646 } 647 648 if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0) 649 { 650 /* Return NaN. */ 651 FLOAT retval = NAN; 652 653 cp += 3; 654 655 /* Match `(n-char-sequence-digit)'. */ 656 if (*cp == L_('(')) 657 { 658 const STRING_TYPE *startp = cp; 659 do 660 ++cp; 661 while ((*cp >= L_('0') && *cp <= L_('9')) 662 || ({ CHAR_TYPE lo = TOLOWER (*cp); 663 lo >= L_('a') && lo <= L_('z'); }) 664 || *cp == L_('_')); 665 666 if (*cp != L_(')')) 667 /* The closing brace is missing. Only match the NAN 668 part. */ 669 cp = startp; 670 else 671 { 672 /* This is a system-dependent way to specify the 673 bitmask used for the NaN. We expect it to be 674 a number which is put in the mantissa of the 675 number. */ 676 STRING_TYPE *endp; 677 unsigned long long int mant; 678 679 mant = STRTOULL (startp + 1, &endp, 0); 680 if (endp == cp) 681 SET_MANTISSA (retval, mant); 682 683 /* Consume the closing brace. */ 684 ++cp; 685 } 686 } 687 688 if (endptr != NULL) 689 *endptr = (STRING_TYPE *) cp; 690 691 return retval; 692 } 693 694 /* It is really a text we do not recognize. */ 695 RETURN (0.0, nptr); 696 } 697 698 /* First look whether we are faced with a hexadecimal number. */ 699 if (c == L_('0') && TOLOWER (cp[1]) == L_('x')) 700 { 701 /* Okay, it is a hexa-decimal number. Remember this and skip 702 the characters. BTW: hexadecimal numbers must not be 703 grouped. */ 704 base = 16; 705 cp += 2; 706 c = *cp; 707 grouping = NULL; 708 } 709 710 /* Record the start of the digits, in case we will check their grouping. */ 711 start_of_digits = startp = cp; 712 713 /* Ignore leading zeroes. This helps us to avoid useless computations. */ 714 #ifdef USE_WIDE_CHAR 715 while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands)) 716 c = *++cp; 717 #else 718 if (__builtin_expect (thousands == NULL, 1)) 719 while (c == '0') 720 c = *++cp; 721 else 722 { 723 /* We also have the multibyte thousands string. */ 724 while (1) 725 { 726 if (c != '0') 727 { 728 for (cnt = 0; thousands[cnt] != '\0'; ++cnt) 729 if (thousands[cnt] != cp[cnt]) 730 break; 731 if (thousands[cnt] != '\0') 732 break; 733 cp += cnt - 1; 734 } 735 c = *++cp; 736 } 737 } 738 #endif 739 740 /* If no other digit but a '0' is found the result is 0.0. 741 Return current read pointer. */ 742 CHAR_TYPE lowc = TOLOWER (c); 743 if (!((c >= L_('0') && c <= L_('9')) 744 || (base == 16 && lowc >= L_('a') && lowc <= L_('f')) 745 || ( 746 #ifdef USE_WIDE_CHAR 747 c == (wint_t) decimal 748 #else 749 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) 750 if (decimal[cnt] != cp[cnt]) 751 break; 752 decimal[cnt] == '\0'; }) 753 #endif 754 /* '0x.' alone is not a valid hexadecimal number. 755 '.' alone is not valid either, but that has been checked 756 already earlier. */ 757 && (base != 16 758 || cp != start_of_digits 759 || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9')) 760 || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]); 761 lo >= L_('a') && lo <= L_('f'); }))) 762 || (base == 16 && (cp != start_of_digits 763 && lowc == L_('p'))) 764 || (base != 16 && lowc == L_('e')))) 765 { 766 #ifdef USE_WIDE_CHAR 767 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, 768 grouping); 769 #else 770 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, 771 grouping); 772 #endif 773 /* If TP is at the start of the digits, there was no correctly 774 grouped prefix of the string; so no number found. */ 775 RETURN (negative ? -0.0 : 0.0, 776 tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp); 777 } 778 779 /* Remember first significant digit and read following characters until the 780 decimal point, exponent character or any non-FP number character. */ 781 startp = cp; 782 dig_no = 0; 783 while (1) 784 { 785 if ((c >= L_('0') && c <= L_('9')) 786 || (base == 16 787 && ({ CHAR_TYPE lo = TOLOWER (c); 788 lo >= L_('a') && lo <= L_('f'); }))) 789 ++dig_no; 790 else 791 { 792 #ifdef USE_WIDE_CHAR 793 if (__builtin_expect ((wint_t) thousands == L'\0', 1) 794 || c != (wint_t) thousands) 795 /* Not a digit or separator: end of the integer part. */ 796 break; 797 #else 798 if (__builtin_expect (thousands == NULL, 1)) 799 break; 800 else 801 { 802 for (cnt = 0; thousands[cnt] != '\0'; ++cnt) 803 if (thousands[cnt] != cp[cnt]) 804 break; 805 if (thousands[cnt] != '\0') 806 break; 807 cp += cnt - 1; 808 } 809 #endif 810 } 811 c = *++cp; 812 } 813 814 if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits) 815 { 816 /* Check the grouping of the digits. */ 817 #ifdef USE_WIDE_CHAR 818 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, 819 grouping); 820 #else 821 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, 822 grouping); 823 #endif 824 if (cp != tp) 825 { 826 /* Less than the entire string was correctly grouped. */ 827 828 if (tp == start_of_digits) 829 /* No valid group of numbers at all: no valid number. */ 830 RETURN (0.0, nptr); 831 832 if (tp < startp) 833 /* The number is validly grouped, but consists 834 only of zeroes. The whole value is zero. */ 835 RETURN (negative ? -0.0 : 0.0, tp); 836 837 /* Recompute DIG_NO so we won't read more digits than 838 are properly grouped. */ 839 cp = tp; 840 dig_no = 0; 841 for (tp = startp; tp < cp; ++tp) 842 if (*tp >= L_('0') && *tp <= L_('9')) 843 ++dig_no; 844 845 int_no = dig_no; 846 lead_zero = 0; 847 848 goto number_parsed; 849 } 850 } 851 852 /* We have the number of digits in the integer part. Whether these 853 are all or any is really a fractional digit will be decided 854 later. */ 855 int_no = dig_no; 856 lead_zero = int_no == 0 ? (size_t) -1 : 0; 857 858 /* Read the fractional digits. A special case are the 'american 859 style' numbers like `16.' i.e. with decimal point but without 860 trailing digits. */ 861 if ( 862 #ifdef USE_WIDE_CHAR 863 c == (wint_t) decimal 864 #else 865 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) 866 if (decimal[cnt] != cp[cnt]) 867 break; 868 decimal[cnt] == '\0'; }) 869 #endif 870 ) 871 { 872 cp += decimal_len; 873 c = *cp; 874 while ((c >= L_('0') && c <= L_('9')) || 875 (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c); 876 lo >= L_('a') && lo <= L_('f'); }))) 877 { 878 if (c != L_('0') && lead_zero == (size_t) -1) 879 lead_zero = dig_no - int_no; 880 ++dig_no; 881 c = *++cp; 882 } 883 } 884 assert (dig_no <= (uintmax_t) INTMAX_MAX); 885 886 /* Remember start of exponent (if any). */ 887 expp = cp; 888 889 /* Read exponent. */ 890 lowc = TOLOWER (c); 891 if ((base == 16 && lowc == L_('p')) 892 || (base != 16 && lowc == L_('e'))) 893 { 894 int exp_negative = 0; 895 896 c = *++cp; 897 if (c == L_('-')) 898 { 899 exp_negative = 1; 900 c = *++cp; 901 } 902 else if (c == L_('+')) 903 c = *++cp; 904 905 if (c >= L_('0') && c <= L_('9')) 906 { 907 intmax_t exp_limit; 908 909 /* Get the exponent limit. */ 910 if (base == 16) 911 { 912 if (exp_negative) 913 { 914 assert (int_no <= (uintmax_t) (INTMAX_MAX 915 + MIN_EXP - MANT_DIG) / 4); 916 exp_limit = -MIN_EXP + MANT_DIG + 4 * (intmax_t) int_no; 917 } 918 else 919 { 920 if (int_no) 921 { 922 assert (lead_zero == 0 923 && int_no <= (uintmax_t) INTMAX_MAX / 4); 924 exp_limit = MAX_EXP - 4 * (intmax_t) int_no + 3; 925 } 926 else if (lead_zero == (size_t) -1) 927 { 928 /* The number is zero and this limit is 929 arbitrary. */ 930 exp_limit = MAX_EXP + 3; 931 } 932 else 933 { 934 assert (lead_zero 935 <= (uintmax_t) (INTMAX_MAX - MAX_EXP - 3) / 4); 936 exp_limit = (MAX_EXP 937 + 4 * (intmax_t) lead_zero 938 + 3); 939 } 940 } 941 } 942 else 943 { 944 if (exp_negative) 945 { 946 assert (int_no 947 <= (uintmax_t) (INTMAX_MAX + MIN_10_EXP - MANT_DIG)); 948 exp_limit = -MIN_10_EXP + MANT_DIG + (intmax_t) int_no; 949 } 950 else 951 { 952 if (int_no) 953 { 954 assert (lead_zero == 0 955 && int_no <= (uintmax_t) INTMAX_MAX); 956 exp_limit = MAX_10_EXP - (intmax_t) int_no + 1; 957 } 958 else if (lead_zero == (size_t) -1) 959 { 960 /* The number is zero and this limit is 961 arbitrary. */ 962 exp_limit = MAX_10_EXP + 1; 963 } 964 else 965 { 966 assert (lead_zero 967 <= (uintmax_t) (INTMAX_MAX - MAX_10_EXP - 1)); 968 exp_limit = MAX_10_EXP + (intmax_t) lead_zero + 1; 969 } 970 } 971 } 972 973 if (exp_limit < 0) 974 exp_limit = 0; 975 976 do 977 { 978 if (__builtin_expect ((exponent > exp_limit / 10 979 || (exponent == exp_limit / 10 980 && c - L_('0') > exp_limit % 10)), 0)) 981 /* The exponent is too large/small to represent a valid 982 number. */ 983 { 984 FLOAT result; 985 986 /* We have to take care for special situation: a joker 987 might have written "0.0e100000" which is in fact 988 zero. */ 989 if (lead_zero == (size_t) -1) 990 result = negative ? -0.0 : 0.0; 991 else 992 { 993 /* Overflow or underflow. */ 994 #if defined HAVE_ERRNO_H && defined ERANGE 995 errno = ERANGE; 996 #endif 997 result = (exp_negative ? (negative ? -0.0 : 0.0) : 998 negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL); 999 } 1000 1001 /* Accept all following digits as part of the exponent. */ 1002 do 1003 ++cp; 1004 while (*cp >= L_('0') && *cp <= L_('9')); 1005 1006 RETURN (result, cp); 1007 /* NOTREACHED */ 1008 } 1009 1010 exponent *= 10; 1011 exponent += c - L_('0'); 1012 1013 c = *++cp; 1014 } 1015 while (c >= L_('0') && c <= L_('9')); 1016 1017 if (exp_negative) 1018 exponent = -exponent; 1019 } 1020 else 1021 cp = expp; 1022 } 1023 1024 /* We don't want to have to work with trailing zeroes after the radix. */ 1025 if (dig_no > int_no) 1026 { 1027 while (expp[-1] == L_('0')) 1028 { 1029 --expp; 1030 --dig_no; 1031 } 1032 assert (dig_no >= int_no); 1033 } 1034 1035 if (dig_no == int_no && dig_no > 0 && exponent < 0) 1036 do 1037 { 1038 while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1]))) 1039 --expp; 1040 1041 if (expp[-1] != L_('0')) 1042 break; 1043 1044 --expp; 1045 --dig_no; 1046 --int_no; 1047 exponent += base == 16 ? 4 : 1; 1048 } 1049 while (dig_no > 0 && exponent < 0); 1050 1051 number_parsed: 1052 1053 /* The whole string is parsed. Store the address of the next character. */ 1054 if (endptr) 1055 *endptr = (STRING_TYPE *) cp; 1056 1057 if (dig_no == 0) 1058 return negative ? -0.0 : 0.0; 1059 1060 if (lead_zero) 1061 { 1062 /* Find the decimal point */ 1063 #ifdef USE_WIDE_CHAR 1064 while (*startp != decimal) 1065 ++startp; 1066 #else 1067 while (1) 1068 { 1069 if (*startp == decimal[0]) 1070 { 1071 for (cnt = 1; decimal[cnt] != '\0'; ++cnt) 1072 if (decimal[cnt] != startp[cnt]) 1073 break; 1074 if (decimal[cnt] == '\0') 1075 break; 1076 } 1077 ++startp; 1078 } 1079 #endif 1080 startp += lead_zero + decimal_len; 1081 assert (lead_zero <= (base == 16 1082 ? (uintmax_t) INTMAX_MAX / 4 1083 : (uintmax_t) INTMAX_MAX)); 1084 assert (lead_zero <= (base == 16 1085 ? ((uintmax_t) exponent 1086 - (uintmax_t) INTMAX_MIN) / 4 1087 : ((uintmax_t) exponent - (uintmax_t) INTMAX_MIN))); 1088 exponent -= base == 16 ? 4 * (intmax_t) lead_zero : (intmax_t) lead_zero; 1089 dig_no -= lead_zero; 1090 } 1091 1092 /* If the BASE is 16 we can use a simpler algorithm. */ 1093 if (base == 16) 1094 { 1095 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3, 1096 4, 4, 4, 4, 4, 4, 4, 4 }; 1097 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB; 1098 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB; 1099 mp_limb_t val; 1100 1101 while (!ISXDIGIT (*startp)) 1102 ++startp; 1103 while (*startp == L_('0')) 1104 ++startp; 1105 if (ISDIGIT (*startp)) 1106 val = *startp++ - L_('0'); 1107 else 1108 val = 10 + TOLOWER (*startp++) - L_('a'); 1109 bits = nbits[val]; 1110 /* We cannot have a leading zero. */ 1111 assert (bits != 0); 1112 1113 if (pos + 1 >= 4 || pos + 1 >= bits) 1114 { 1115 /* We don't have to care for wrapping. This is the normal 1116 case so we add the first clause in the `if' expression as 1117 an optimization. It is a compile-time constant and so does 1118 not cost anything. */ 1119 retval[idx] = val << (pos - bits + 1); 1120 pos -= bits; 1121 } 1122 else 1123 { 1124 retval[idx--] = val >> (bits - pos - 1); 1125 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1)); 1126 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1); 1127 } 1128 1129 /* Adjust the exponent for the bits we are shifting in. */ 1130 assert (int_no <= (uintmax_t) (exponent < 0 1131 ? (INTMAX_MAX - bits + 1) / 4 1132 : (INTMAX_MAX - exponent - bits + 1) / 4)); 1133 exponent += bits - 1 + ((intmax_t) int_no - 1) * 4; 1134 1135 while (--dig_no > 0 && idx >= 0) 1136 { 1137 if (!ISXDIGIT (*startp)) 1138 startp += decimal_len; 1139 if (ISDIGIT (*startp)) 1140 val = *startp++ - L_('0'); 1141 else 1142 val = 10 + TOLOWER (*startp++) - L_('a'); 1143 1144 if (pos + 1 >= 4) 1145 { 1146 retval[idx] |= val << (pos - 4 + 1); 1147 pos -= 4; 1148 } 1149 else 1150 { 1151 retval[idx--] |= val >> (4 - pos - 1); 1152 val <<= BITS_PER_MP_LIMB - (4 - pos - 1); 1153 if (idx < 0) 1154 { 1155 int rest_nonzero = 0; 1156 while (--dig_no > 0) 1157 { 1158 if (*startp != L_('0')) 1159 { 1160 rest_nonzero = 1; 1161 break; 1162 } 1163 startp++; 1164 } 1165 return round_and_return (retval, exponent, negative, val, 1166 BITS_PER_MP_LIMB - 1, rest_nonzero); 1167 } 1168 1169 retval[idx] = val; 1170 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1); 1171 } 1172 } 1173 1174 /* We ran out of digits. */ 1175 MPN_ZERO (retval, idx); 1176 1177 return round_and_return (retval, exponent, negative, 0, 0, 0); 1178 } 1179 1180 /* Now we have the number of digits in total and the integer digits as well 1181 as the exponent and its sign. We can decide whether the read digits are 1182 really integer digits or belong to the fractional part; i.e. we normalize 1183 123e-2 to 1.23. */ 1184 { 1185 register intmax_t incr = (exponent < 0 1186 ? MAX (-(intmax_t) int_no, exponent) 1187 : MIN ((intmax_t) dig_no - (intmax_t) int_no, 1188 exponent)); 1189 int_no += incr; 1190 exponent -= incr; 1191 } 1192 1193 if (__builtin_expect (exponent > MAX_10_EXP + 1 - (intmax_t) int_no, 0)) 1194 return overflow_value (negative); 1195 1196 if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0)) 1197 return underflow_value (negative); 1198 1199 if (int_no > 0) 1200 { 1201 /* Read the integer part as a multi-precision number to NUM. */ 1202 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent 1203 #ifndef USE_WIDE_CHAR 1204 , decimal, decimal_len, thousands 1205 #endif 1206 ); 1207 1208 if (exponent > 0) 1209 { 1210 /* We now multiply the gained number by the given power of ten. */ 1211 mp_limb_t *psrc = num; 1212 mp_limb_t *pdest = den; 1213 int expbit = 1; 1214 const struct mp_power *ttab = &_fpioconst_pow10[0]; 1215 1216 do 1217 { 1218 if ((exponent & expbit) != 0) 1219 { 1220 size_t size = ttab->arraysize - _FPIO_CONST_OFFSET; 1221 mp_limb_t cy; 1222 exponent ^= expbit; 1223 1224 /* FIXME: not the whole multiplication has to be 1225 done. If we have the needed number of bits we 1226 only need the information whether more non-zero 1227 bits follow. */ 1228 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET) 1229 cy = mpn_mul (pdest, psrc, numsize, 1230 &__tens[ttab->arrayoff 1231 + _FPIO_CONST_OFFSET], 1232 size); 1233 else 1234 cy = mpn_mul (pdest, &__tens[ttab->arrayoff 1235 + _FPIO_CONST_OFFSET], 1236 size, psrc, numsize); 1237 numsize += size; 1238 if (cy == 0) 1239 --numsize; 1240 (void) SWAP (psrc, pdest); 1241 } 1242 expbit <<= 1; 1243 ++ttab; 1244 } 1245 while (exponent != 0); 1246 1247 if (psrc == den) 1248 memcpy (num, den, numsize * sizeof (mp_limb_t)); 1249 } 1250 1251 /* Determine how many bits of the result we already have. */ 1252 count_leading_zeros (bits, num[numsize - 1]); 1253 bits = numsize * BITS_PER_MP_LIMB - bits; 1254 1255 /* Now we know the exponent of the number in base two. 1256 Check it against the maximum possible exponent. */ 1257 if (__builtin_expect (bits > MAX_EXP, 0)) 1258 return overflow_value (negative); 1259 1260 /* We have already the first BITS bits of the result. Together with 1261 the information whether more non-zero bits follow this is enough 1262 to determine the result. */ 1263 if (bits > MANT_DIG) 1264 { 1265 int i; 1266 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB; 1267 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB; 1268 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1 1269 : least_idx; 1270 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1 1271 : least_bit - 1; 1272 1273 if (least_bit == 0) 1274 memcpy (retval, &num[least_idx], 1275 RETURN_LIMB_SIZE * sizeof (mp_limb_t)); 1276 else 1277 { 1278 for (i = least_idx; i < numsize - 1; ++i) 1279 retval[i - least_idx] = (num[i] >> least_bit) 1280 | (num[i + 1] 1281 << (BITS_PER_MP_LIMB - least_bit)); 1282 if (i - least_idx < RETURN_LIMB_SIZE) 1283 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit; 1284 } 1285 1286 /* Check whether any limb beside the ones in RETVAL are non-zero. */ 1287 for (i = 0; num[i] == 0; ++i) 1288 ; 1289 1290 return round_and_return (retval, bits - 1, negative, 1291 num[round_idx], round_bit, 1292 int_no < dig_no || i < round_idx); 1293 /* NOTREACHED */ 1294 } 1295 else if (dig_no == int_no) 1296 { 1297 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; 1298 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB; 1299 1300 if (target_bit == is_bit) 1301 { 1302 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num, 1303 numsize * sizeof (mp_limb_t)); 1304 /* FIXME: the following loop can be avoided if we assume a 1305 maximal MANT_DIG value. */ 1306 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); 1307 } 1308 else if (target_bit > is_bit) 1309 { 1310 (void) mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize], 1311 num, numsize, target_bit - is_bit); 1312 /* FIXME: the following loop can be avoided if we assume a 1313 maximal MANT_DIG value. */ 1314 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); 1315 } 1316 else 1317 { 1318 mp_limb_t cy; 1319 assert (numsize < RETURN_LIMB_SIZE); 1320 1321 cy = mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize], 1322 num, numsize, is_bit - target_bit); 1323 retval[RETURN_LIMB_SIZE - numsize - 1] = cy; 1324 /* FIXME: the following loop can be avoided if we assume a 1325 maximal MANT_DIG value. */ 1326 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1); 1327 } 1328 1329 return round_and_return (retval, bits - 1, negative, 0, 0, 0); 1330 /* NOTREACHED */ 1331 } 1332 1333 /* Store the bits we already have. */ 1334 memcpy (retval, num, numsize * sizeof (mp_limb_t)); 1335 #if RETURN_LIMB_SIZE > 1 1336 if (numsize < RETURN_LIMB_SIZE) 1337 # if RETURN_LIMB_SIZE == 2 1338 retval[numsize] = 0; 1339 # else 1340 MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize); 1341 # endif 1342 #endif 1343 } 1344 1345 /* We have to compute at least some of the fractional digits. */ 1346 { 1347 /* We construct a fraction and the result of the division gives us 1348 the needed digits. The denominator is 1.0 multiplied by the 1349 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and 1350 123e-6 gives 123 / 1000000. */ 1351 1352 int expbit; 1353 int neg_exp; 1354 int more_bits; 1355 int need_frac_digits; 1356 mp_limb_t cy; 1357 mp_limb_t *psrc = den; 1358 mp_limb_t *pdest = num; 1359 const struct mp_power *ttab = &_fpioconst_pow10[0]; 1360 1361 assert (dig_no > int_no 1362 && exponent <= 0 1363 && exponent >= MIN_10_EXP - (DIG + 1)); 1364 1365 /* We need to compute MANT_DIG - BITS fractional bits that lie 1366 within the mantissa of the result, the following bit for 1367 rounding, and to know whether any subsequent bit is 0. 1368 Computing a bit with value 2^-n means looking at n digits after 1369 the decimal point. */ 1370 if (bits > 0) 1371 { 1372 /* The bits required are those immediately after the point. */ 1373 assert (int_no > 0 && exponent == 0); 1374 need_frac_digits = 1 + MANT_DIG - bits; 1375 } 1376 else 1377 { 1378 /* The number is in the form .123eEXPONENT. */ 1379 assert (int_no == 0 && *startp != L_('0')); 1380 /* The number is at least 10^(EXPONENT-1), and 10^3 < 1381 2^10. */ 1382 int neg_exp_2 = ((1 - exponent) * 10) / 3 + 1; 1383 /* The number is at least 2^-NEG_EXP_2. We need up to 1384 MANT_DIG bits following that bit. */ 1385 need_frac_digits = neg_exp_2 + MANT_DIG; 1386 /* However, we never need bits beyond 1/4 ulp of the smallest 1387 representable value. (That 1/4 ulp bit is only needed to 1388 determine tinyness on machines where tinyness is determined 1389 after rounding.) */ 1390 if (need_frac_digits > MANT_DIG - MIN_EXP + 2) 1391 need_frac_digits = MANT_DIG - MIN_EXP + 2; 1392 /* At this point, NEED_FRAC_DIGITS is the total number of 1393 digits needed after the point, but some of those may be 1394 leading 0s. */ 1395 need_frac_digits += exponent; 1396 /* Any cases underflowing enough that none of the fractional 1397 digits are needed should have been caught earlier (such 1398 cases are on the order of 10^-n or smaller where 2^-n is 1399 the least subnormal). */ 1400 assert (need_frac_digits > 0); 1401 } 1402 1403 if (need_frac_digits > (intmax_t) dig_no - (intmax_t) int_no) 1404 need_frac_digits = (intmax_t) dig_no - (intmax_t) int_no; 1405 1406 if ((intmax_t) dig_no > (intmax_t) int_no + need_frac_digits) 1407 { 1408 dig_no = int_no + need_frac_digits; 1409 more_bits = 1; 1410 } 1411 else 1412 more_bits = 0; 1413 1414 neg_exp = (intmax_t) dig_no - (intmax_t) int_no - exponent; 1415 1416 /* Construct the denominator. */ 1417 densize = 0; 1418 expbit = 1; 1419 do 1420 { 1421 if ((neg_exp & expbit) != 0) 1422 { 1423 mp_limb_t cy; 1424 neg_exp ^= expbit; 1425 1426 if (densize == 0) 1427 { 1428 densize = ttab->arraysize - _FPIO_CONST_OFFSET; 1429 memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET], 1430 densize * sizeof (mp_limb_t)); 1431 } 1432 else 1433 { 1434 cy = mpn_mul (pdest, &__tens[ttab->arrayoff 1435 + _FPIO_CONST_OFFSET], 1436 ttab->arraysize - _FPIO_CONST_OFFSET, 1437 psrc, densize); 1438 densize += ttab->arraysize - _FPIO_CONST_OFFSET; 1439 if (cy == 0) 1440 --densize; 1441 (void) SWAP (psrc, pdest); 1442 } 1443 } 1444 expbit <<= 1; 1445 ++ttab; 1446 } 1447 while (neg_exp != 0); 1448 1449 if (psrc == num) 1450 memcpy (den, num, densize * sizeof (mp_limb_t)); 1451 1452 /* Read the fractional digits from the string. */ 1453 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent 1454 #ifndef USE_WIDE_CHAR 1455 , decimal, decimal_len, thousands 1456 #endif 1457 ); 1458 1459 /* We now have to shift both numbers so that the highest bit in the 1460 denominator is set. In the same process we copy the numerator to 1461 a high place in the array so that the division constructs the wanted 1462 digits. This is done by a "quasi fix point" number representation. 1463 1464 num: ddddddddddd . 0000000000000000000000 1465 |--- m ---| 1466 den: ddddddddddd n >= m 1467 |--- n ---| 1468 */ 1469 1470 count_leading_zeros (cnt, den[densize - 1]); 1471 1472 if (cnt > 0) 1473 { 1474 /* Don't call `mpn_shift' with a count of zero since the specification 1475 does not allow this. */ 1476 (void) mpn_lshift (den, den, densize, cnt); 1477 cy = mpn_lshift (num, num, numsize, cnt); 1478 if (cy != 0) 1479 num[numsize++] = cy; 1480 } 1481 1482 /* Now we are ready for the division. But it is not necessary to 1483 do a full multi-precision division because we only need a small 1484 number of bits for the result. So we do not use mpn_divmod 1485 here but instead do the division here by hand and stop whenever 1486 the needed number of bits is reached. The code itself comes 1487 from the GNU MP Library by Torbj\"orn Granlund. */ 1488 1489 exponent = bits; 1490 1491 switch (densize) 1492 { 1493 case 1: 1494 { 1495 mp_limb_t d, n, quot; 1496 int used = 0; 1497 1498 n = num[0]; 1499 d = den[0]; 1500 assert (numsize == 1 && n < d); 1501 1502 do 1503 { 1504 udiv_qrnnd (quot, n, n, 0, d); 1505 1506 #define got_limb \ 1507 if (bits == 0) \ 1508 { \ 1509 register int cnt; \ 1510 if (quot == 0) \ 1511 cnt = BITS_PER_MP_LIMB; \ 1512 else \ 1513 count_leading_zeros (cnt, quot); \ 1514 exponent -= cnt; \ 1515 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \ 1516 { \ 1517 used = MANT_DIG + cnt; \ 1518 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \ 1519 bits = MANT_DIG + 1; \ 1520 } \ 1521 else \ 1522 { \ 1523 /* Note that we only clear the second element. */ \ 1524 /* The conditional is determined at compile time. */ \ 1525 if (RETURN_LIMB_SIZE > 1) \ 1526 retval[1] = 0; \ 1527 retval[0] = quot; \ 1528 bits = -cnt; \ 1529 } \ 1530 } \ 1531 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \ 1532 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \ 1533 quot); \ 1534 else \ 1535 { \ 1536 used = MANT_DIG - bits; \ 1537 if (used > 0) \ 1538 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \ 1539 } \ 1540 bits += BITS_PER_MP_LIMB 1541 1542 got_limb; 1543 } 1544 while (bits <= MANT_DIG); 1545 1546 return round_and_return (retval, exponent - 1, negative, 1547 quot, BITS_PER_MP_LIMB - 1 - used, 1548 more_bits || n != 0); 1549 } 1550 case 2: 1551 { 1552 mp_limb_t d0, d1, n0, n1; 1553 mp_limb_t quot = 0; 1554 int used = 0; 1555 1556 d0 = den[0]; 1557 d1 = den[1]; 1558 1559 if (numsize < densize) 1560 { 1561 if (num[0] >= d1) 1562 { 1563 /* The numerator of the number occupies fewer bits than 1564 the denominator but the one limb is bigger than the 1565 high limb of the numerator. */ 1566 n1 = 0; 1567 n0 = num[0]; 1568 } 1569 else 1570 { 1571 if (bits <= 0) 1572 exponent -= BITS_PER_MP_LIMB; 1573 else 1574 { 1575 if (bits + BITS_PER_MP_LIMB <= MANT_DIG) 1576 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, 1577 BITS_PER_MP_LIMB, 0); 1578 else 1579 { 1580 used = MANT_DIG - bits; 1581 if (used > 0) 1582 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); 1583 } 1584 bits += BITS_PER_MP_LIMB; 1585 } 1586 n1 = num[0]; 1587 n0 = 0; 1588 } 1589 } 1590 else 1591 { 1592 n1 = num[1]; 1593 n0 = num[0]; 1594 } 1595 1596 while (bits <= MANT_DIG) 1597 { 1598 mp_limb_t r; 1599 1600 if (n1 == d1) 1601 { 1602 /* QUOT should be either 111..111 or 111..110. We need 1603 special treatment of this rare case as normal division 1604 would give overflow. */ 1605 quot = ~(mp_limb_t) 0; 1606 1607 r = n0 + d1; 1608 if (r < d1) /* Carry in the addition? */ 1609 { 1610 add_ssaaaa (n1, n0, r - d0, 0, 0, d0); 1611 goto have_quot; 1612 } 1613 n1 = d0 - (d0 != 0); 1614 n0 = -d0; 1615 } 1616 else 1617 { 1618 udiv_qrnnd (quot, r, n1, n0, d1); 1619 umul_ppmm (n1, n0, d0, quot); 1620 } 1621 1622 q_test: 1623 if (n1 > r || (n1 == r && n0 > 0)) 1624 { 1625 /* The estimated QUOT was too large. */ 1626 --quot; 1627 1628 sub_ddmmss (n1, n0, n1, n0, 0, d0); 1629 r += d1; 1630 if (r >= d1) /* If not carry, test QUOT again. */ 1631 goto q_test; 1632 } 1633 sub_ddmmss (n1, n0, r, 0, n1, n0); 1634 1635 have_quot: 1636 got_limb; 1637 } 1638 1639 return round_and_return (retval, exponent - 1, negative, 1640 quot, BITS_PER_MP_LIMB - 1 - used, 1641 more_bits || n1 != 0 || n0 != 0); 1642 } 1643 default: 1644 { 1645 int i; 1646 mp_limb_t cy, dX, d1, n0, n1; 1647 mp_limb_t quot = 0; 1648 int used = 0; 1649 1650 dX = den[densize - 1]; 1651 d1 = den[densize - 2]; 1652 1653 /* The division does not work if the upper limb of the two-limb 1654 numerator is greater than the denominator. */ 1655 if (mpn_cmp (num, &den[densize - numsize], numsize) > 0) 1656 num[numsize++] = 0; 1657 1658 if (numsize < densize) 1659 { 1660 mp_size_t empty = densize - numsize; 1661 register int i; 1662 1663 if (bits <= 0) 1664 exponent -= empty * BITS_PER_MP_LIMB; 1665 else 1666 { 1667 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG) 1668 { 1669 /* We make a difference here because the compiler 1670 cannot optimize the `else' case that good and 1671 this reflects all currently used FLOAT types 1672 and GMP implementations. */ 1673 #if RETURN_LIMB_SIZE <= 2 1674 assert (empty == 1); 1675 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, 1676 BITS_PER_MP_LIMB, 0); 1677 #else 1678 for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i) 1679 retval[i] = retval[i - empty]; 1680 while (i >= 0) 1681 retval[i--] = 0; 1682 #endif 1683 } 1684 else 1685 { 1686 used = MANT_DIG - bits; 1687 if (used >= BITS_PER_MP_LIMB) 1688 { 1689 register int i; 1690 (void) mpn_lshift (&retval[used 1691 / BITS_PER_MP_LIMB], 1692 retval, 1693 (RETURN_LIMB_SIZE 1694 - used / BITS_PER_MP_LIMB), 1695 used % BITS_PER_MP_LIMB); 1696 for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i) 1697 retval[i] = 0; 1698 } 1699 else if (used > 0) 1700 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); 1701 } 1702 bits += empty * BITS_PER_MP_LIMB; 1703 } 1704 for (i = numsize; i > 0; --i) 1705 num[i + empty] = num[i - 1]; 1706 MPN_ZERO (num, empty + 1); 1707 } 1708 else 1709 { 1710 int i; 1711 assert (numsize == densize); 1712 for (i = numsize; i > 0; --i) 1713 num[i] = num[i - 1]; 1714 num[0] = 0; 1715 } 1716 1717 den[densize] = 0; 1718 n0 = num[densize]; 1719 1720 while (bits <= MANT_DIG) 1721 { 1722 if (n0 == dX) 1723 /* This might over-estimate QUOT, but it's probably not 1724 worth the extra code here to find out. */ 1725 quot = ~(mp_limb_t) 0; 1726 else 1727 { 1728 mp_limb_t r; 1729 1730 udiv_qrnnd (quot, r, n0, num[densize - 1], dX); 1731 umul_ppmm (n1, n0, d1, quot); 1732 1733 while (n1 > r || (n1 == r && n0 > num[densize - 2])) 1734 { 1735 --quot; 1736 r += dX; 1737 if (r < dX) /* I.e. "carry in previous addition?" */ 1738 break; 1739 n1 -= n0 < d1; 1740 n0 -= d1; 1741 } 1742 } 1743 1744 /* Possible optimization: We already have (q * n0) and (1 * n1) 1745 after the calculation of QUOT. Taking advantage of this, we 1746 could make this loop make two iterations less. */ 1747 1748 cy = mpn_submul_1 (num, den, densize + 1, quot); 1749 1750 if (num[densize] != cy) 1751 { 1752 cy = mpn_add_n (num, num, den, densize); 1753 assert (cy != 0); 1754 --quot; 1755 } 1756 n0 = num[densize] = num[densize - 1]; 1757 for (i = densize - 1; i > 0; --i) 1758 num[i] = num[i - 1]; 1759 num[0] = 0; 1760 1761 got_limb; 1762 } 1763 1764 for (i = densize; num[i] == 0 && i >= 0; --i) 1765 ; 1766 return round_and_return (retval, exponent - 1, negative, 1767 quot, BITS_PER_MP_LIMB - 1 - used, 1768 more_bits || i >= 0); 1769 } 1770 } 1771 } 1772 1773 /* NOTREACHED */ 1774 } 1775