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