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