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 && (c < TOLOWER (L_('a')) || c > TOLOWER (L_('f')))) 679 #ifdef USE_WIDE_CHAR 680 && c != decimal 681 #else 682 && ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) 683 if (decimal[cnt] != cp[cnt]) 684 break; 685 decimal[cnt] != '\0'; }) 686 #endif 687 && (base == 16 && (cp == start_of_digits || TOLOWER (c) != L_('p'))) 688 && (base != 16 && TOLOWER (c) != L_('e'))) 689 { 690 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping); 691 /* If TP is at the start of the digits, there was no correctly 692 grouped prefix of the string; so no number found. */ 693 RETURN (0.0, tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp); 694 } 695 696 /* Remember first significant digit and read following characters until the 697 decimal point, exponent character or any non-FP number character. */ 698 startp = cp; 699 dig_no = 0; 700 while (1) 701 { 702 if ((c >= L_('0') && c <= L_('9')) 703 || (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f'))) 704 ++dig_no; 705 else 706 { 707 #ifdef USE_WIDE_CHAR 708 if (thousands == L'\0' || c != thousands) 709 /* Not a digit or separator: end of the integer part. */ 710 break; 711 #else 712 if (thousands == NULL) 713 break; 714 else 715 { 716 for (cnt = 0; thousands[cnt] != '\0'; ++cnt) 717 if (thousands[cnt] != cp[cnt]) 718 break; 719 if (thousands[cnt] != '\0') 720 break; 721 } 722 #endif 723 } 724 c = *++cp; 725 } 726 727 if (grouping && dig_no > 0) 728 { 729 /* Check the grouping of the digits. */ 730 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping); 731 if (cp != tp) 732 { 733 /* Less than the entire string was correctly grouped. */ 734 735 if (tp == start_of_digits) 736 /* No valid group of numbers at all: no valid number. */ 737 RETURN (0.0, nptr); 738 739 if (tp < startp) 740 /* The number is validly grouped, but consists 741 only of zeroes. The whole value is zero. */ 742 RETURN (0.0, tp); 743 744 /* Recompute DIG_NO so we won't read more digits than 745 are properly grouped. */ 746 cp = tp; 747 dig_no = 0; 748 for (tp = startp; tp < cp; ++tp) 749 if (*tp >= L_('0') && *tp <= L_('9')) 750 ++dig_no; 751 752 int_no = dig_no; 753 lead_zero = 0; 754 755 goto number_parsed; 756 } 757 } 758 759 /* We have the number digits in the integer part. Whether these are all or 760 any is really a fractional digit will be decided later. */ 761 int_no = dig_no; 762 lead_zero = int_no == 0 ? -1 : 0; 763 764 /* Read the fractional digits. A special case are the 'american style' 765 numbers like `16.' i.e. with decimal but without trailing digits. */ 766 if ( 767 #ifdef USE_WIDE_CHAR 768 c == decimal 769 #else 770 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) 771 if (decimal[cnt] != cp[cnt]) 772 break; 773 decimal[cnt] == '\0'; }) 774 #endif 775 ) 776 { 777 cp += decimal_len; 778 c = *cp; 779 while ((c >= L_('0') && c <= L_('9')) || 780 (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f'))) 781 { 782 if (c != L_('0') && lead_zero == -1) 783 lead_zero = dig_no - int_no; 784 ++dig_no; 785 c = *++cp; 786 } 787 } 788 789 /* Remember start of exponent (if any). */ 790 expp = cp; 791 792 /* Read exponent. */ 793 if ((base == 16 && TOLOWER (c) == L_('p')) 794 || (base != 16 && TOLOWER (c) == L_('e'))) 795 { 796 int exp_negative = 0; 797 798 c = *++cp; 799 if (c == L_('-')) 800 { 801 exp_negative = 1; 802 c = *++cp; 803 } 804 else if (c == L_('+')) 805 c = *++cp; 806 807 if (c >= L_('0') && c <= L_('9')) 808 { 809 int exp_limit; 810 811 /* Get the exponent limit. */ 812 if (base == 16) 813 exp_limit = (exp_negative ? 814 -MIN_EXP + MANT_DIG + 4 * int_no : 815 MAX_EXP - 4 * int_no + lead_zero); 816 else 817 exp_limit = (exp_negative ? 818 -MIN_10_EXP + MANT_DIG + int_no : 819 MAX_10_EXP - int_no + lead_zero); 820 821 do 822 { 823 exponent *= 10; 824 825 if (exponent > exp_limit) 826 /* The exponent is too large/small to represent a valid 827 number. */ 828 { 829 FLOAT result; 830 831 /* We have to take care for special situation: a joker 832 might have written "0.0e100000" which is in fact 833 zero. */ 834 if (lead_zero == -1) 835 result = negative ? -0.0 : 0.0; 836 else 837 { 838 /* Overflow or underflow. */ 839 __set_errno (ERANGE); 840 result = (exp_negative ? 0.0 : 841 negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL); 842 } 843 844 /* Accept all following digits as part of the exponent. */ 845 do 846 ++cp; 847 while (*cp >= L_('0') && *cp <= L_('9')); 848 849 RETURN (result, cp); 850 /* NOTREACHED */ 851 } 852 853 exponent += c - L_('0'); 854 c = *++cp; 855 } 856 while (c >= L_('0') && c <= L_('9')); 857 858 if (exp_negative) 859 exponent = -exponent; 860 } 861 else 862 cp = expp; 863 } 864 865 /* We don't want to have to work with trailing zeroes after the radix. */ 866 if (dig_no > int_no) 867 { 868 while (expp[-1] == L_('0')) 869 { 870 --expp; 871 --dig_no; 872 } 873 assert (dig_no >= int_no); 874 } 875 876 if (dig_no == int_no && dig_no > 0 && exponent < 0) 877 do 878 { 879 while (expp[-1] < L_('0') || expp[-1] > L_('9')) 880 --expp; 881 882 if (expp[-1] != L_('0')) 883 break; 884 885 --expp; 886 --dig_no; 887 --int_no; 888 ++exponent; 889 } 890 while (dig_no > 0 && exponent < 0); 891 892 number_parsed: 893 894 /* The whole string is parsed. Store the address of the next character. */ 895 if (endptr) 896 *endptr = (STRING_TYPE *) cp; 897 898 if (dig_no == 0) 899 return negative ? -0.0 : 0.0; 900 901 if (lead_zero) 902 { 903 /* Find the decimal point */ 904 #ifdef USE_WIDE_CHAR 905 while (*startp != decimal) 906 ++startp; 907 #else 908 while (1) 909 { 910 if (*startp == decimal[0]) 911 { 912 for (cnt = 1; decimal[cnt] != '\0'; ++cnt) 913 if (decimal[cnt] != startp[cnt]) 914 break; 915 if (decimal[cnt] == '\0') 916 break; 917 } 918 ++startp; 919 } 920 #endif 921 startp += lead_zero + decimal_len; 922 exponent -= base == 16 ? 4 * lead_zero : lead_zero; 923 dig_no -= lead_zero; 924 } 925 926 /* If the BASE is 16 we can use a simpler algorithm. */ 927 if (base == 16) 928 { 929 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3, 930 4, 4, 4, 4, 4, 4, 4, 4 }; 931 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB; 932 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB; 933 mp_limb_t val; 934 935 while (!ISXDIGIT (*startp)) 936 ++startp; 937 while (*startp == L_('0')) 938 ++startp; 939 if (ISDIGIT (*startp)) 940 val = *startp++ - L_('0'); 941 else 942 val = 10 + TOLOWER (*startp++) - L_('a'); 943 bits = nbits[val]; 944 /* We cannot have a leading zero. */ 945 assert (bits != 0); 946 947 if (pos + 1 >= 4 || pos + 1 >= bits) 948 { 949 /* We don't have to care for wrapping. This is the normal 950 case so we add the first clause in the `if' expression as 951 an optimization. It is a compile-time constant and so does 952 not cost anything. */ 953 retval[idx] = val << (pos - bits + 1); 954 pos -= bits; 955 } 956 else 957 { 958 retval[idx--] = val >> (bits - pos - 1); 959 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1)); 960 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1); 961 } 962 963 /* Adjust the exponent for the bits we are shifting in. */ 964 exponent += bits - 1 + (int_no - 1) * 4; 965 966 while (--dig_no > 0 && idx >= 0) 967 { 968 if (!ISXDIGIT (*startp)) 969 startp += decimal_len; 970 if (ISDIGIT (*startp)) 971 val = *startp++ - L_('0'); 972 else 973 val = 10 + TOLOWER (*startp++) - L_('a'); 974 975 if (pos + 1 >= 4) 976 { 977 retval[idx] |= val << (pos - 4 + 1); 978 pos -= 4; 979 } 980 else 981 { 982 retval[idx--] |= val >> (4 - pos - 1); 983 val <<= BITS_PER_MP_LIMB - (4 - pos - 1); 984 if (idx < 0) 985 return round_and_return (retval, exponent, negative, val, 986 BITS_PER_MP_LIMB - 1, dig_no > 0); 987 988 retval[idx] = val; 989 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1); 990 } 991 } 992 993 /* We ran out of digits. */ 994 MPN_ZERO (retval, idx); 995 996 return round_and_return (retval, exponent, negative, 0, 0, 0); 997 } 998 999 /* Now we have the number of digits in total and the integer digits as well 1000 as the exponent and its sign. We can decide whether the read digits are 1001 really integer digits or belong to the fractional part; i.e. we normalize 1002 123e-2 to 1.23. */ 1003 { 1004 register int incr = (exponent < 0 ? MAX (-int_no, exponent) 1005 : MIN (dig_no - int_no, exponent)); 1006 int_no += incr; 1007 exponent -= incr; 1008 } 1009 1010 if (int_no + exponent > MAX_10_EXP + 1) 1011 { 1012 __set_errno (ERANGE); 1013 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; 1014 } 1015 1016 if (exponent < MIN_10_EXP - (DIG + 1)) 1017 { 1018 __set_errno (ERANGE); 1019 return 0.0; 1020 } 1021 1022 if (int_no > 0) 1023 { 1024 /* Read the integer part as a multi-precision number to NUM. */ 1025 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent 1026 #ifndef USE_WIDE_CHAR 1027 , decimal, decimal_len, thousands 1028 #endif 1029 ); 1030 1031 if (exponent > 0) 1032 { 1033 /* We now multiply the gained number by the given power of ten. */ 1034 mp_limb_t *psrc = num; 1035 mp_limb_t *pdest = den; 1036 int expbit = 1; 1037 const struct mp_power *ttab = &_fpioconst_pow10[0]; 1038 1039 do 1040 { 1041 if ((exponent & expbit) != 0) 1042 { 1043 size_t size = ttab->arraysize - _FPIO_CONST_OFFSET; 1044 mp_limb_t cy; 1045 exponent ^= expbit; 1046 1047 /* FIXME: not the whole multiplication has to be 1048 done. If we have the needed number of bits we 1049 only need the information whether more non-zero 1050 bits follow. */ 1051 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET) 1052 cy = __mpn_mul (pdest, psrc, numsize, 1053 &__tens[ttab->arrayoff 1054 + _FPIO_CONST_OFFSET], 1055 size); 1056 else 1057 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff 1058 + _FPIO_CONST_OFFSET], 1059 size, psrc, numsize); 1060 numsize += size; 1061 if (cy == 0) 1062 --numsize; 1063 (void) SWAP (psrc, pdest); 1064 } 1065 expbit <<= 1; 1066 ++ttab; 1067 } 1068 while (exponent != 0); 1069 1070 if (psrc == den) 1071 memcpy (num, den, numsize * sizeof (mp_limb_t)); 1072 } 1073 1074 /* Determine how many bits of the result we already have. */ 1075 count_leading_zeros (bits, num[numsize - 1]); 1076 bits = numsize * BITS_PER_MP_LIMB - bits; 1077 1078 /* Now we know the exponent of the number in base two. 1079 Check it against the maximum possible exponent. */ 1080 if (bits > MAX_EXP) 1081 { 1082 __set_errno (ERANGE); 1083 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; 1084 } 1085 1086 /* We have already the first BITS bits of the result. Together with 1087 the information whether more non-zero bits follow this is enough 1088 to determine the result. */ 1089 if (bits > MANT_DIG) 1090 { 1091 int i; 1092 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB; 1093 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB; 1094 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1 1095 : least_idx; 1096 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1 1097 : least_bit - 1; 1098 1099 if (least_bit == 0) 1100 memcpy (retval, &num[least_idx], 1101 RETURN_LIMB_SIZE * sizeof (mp_limb_t)); 1102 else 1103 { 1104 for (i = least_idx; i < numsize - 1; ++i) 1105 retval[i - least_idx] = (num[i] >> least_bit) 1106 | (num[i + 1] 1107 << (BITS_PER_MP_LIMB - least_bit)); 1108 if (i - least_idx < RETURN_LIMB_SIZE) 1109 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit; 1110 } 1111 1112 /* Check whether any limb beside the ones in RETVAL are non-zero. */ 1113 for (i = 0; num[i] == 0; ++i) 1114 ; 1115 1116 return round_and_return (retval, bits - 1, negative, 1117 num[round_idx], round_bit, 1118 int_no < dig_no || i < round_idx); 1119 /* NOTREACHED */ 1120 } 1121 else if (dig_no == int_no) 1122 { 1123 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; 1124 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB; 1125 1126 if (target_bit == is_bit) 1127 { 1128 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num, 1129 numsize * sizeof (mp_limb_t)); 1130 /* FIXME: the following loop can be avoided if we assume a 1131 maximal MANT_DIG value. */ 1132 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); 1133 } 1134 else if (target_bit > is_bit) 1135 { 1136 (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize], 1137 num, numsize, target_bit - is_bit); 1138 /* FIXME: the following loop can be avoided if we assume a 1139 maximal MANT_DIG value. */ 1140 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); 1141 } 1142 else 1143 { 1144 mp_limb_t cy; 1145 assert (numsize < RETURN_LIMB_SIZE); 1146 1147 cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize], 1148 num, numsize, is_bit - target_bit); 1149 retval[RETURN_LIMB_SIZE - numsize - 1] = cy; 1150 /* FIXME: the following loop can be avoided if we assume a 1151 maximal MANT_DIG value. */ 1152 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1); 1153 } 1154 1155 return round_and_return (retval, bits - 1, negative, 0, 0, 0); 1156 /* NOTREACHED */ 1157 } 1158 1159 /* Store the bits we already have. */ 1160 memcpy (retval, num, numsize * sizeof (mp_limb_t)); 1161 #if RETURN_LIMB_SIZE > 1 1162 if (numsize < RETURN_LIMB_SIZE) 1163 retval[numsize] = 0; 1164 #endif 1165 } 1166 1167 /* We have to compute at least some of the fractional digits. */ 1168 { 1169 /* We construct a fraction and the result of the division gives us 1170 the needed digits. The denominator is 1.0 multiplied by the 1171 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and 1172 123e-6 gives 123 / 1000000. */ 1173 1174 int expbit; 1175 int neg_exp; 1176 int more_bits; 1177 mp_limb_t cy; 1178 mp_limb_t *psrc = den; 1179 mp_limb_t *pdest = num; 1180 const struct mp_power *ttab = &_fpioconst_pow10[0]; 1181 1182 assert (dig_no > int_no && exponent <= 0); 1183 1184 1185 /* For the fractional part we need not process too many digits. One 1186 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute 1187 ceil(BITS / 3) =: N 1188 digits we should have enough bits for the result. The remaining 1189 decimal digits give us the information that more bits are following. 1190 This can be used while rounding. (One added as a safety margin.) */ 1191 if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1) 1192 { 1193 dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1; 1194 more_bits = 1; 1195 } 1196 else 1197 more_bits = 0; 1198 1199 neg_exp = dig_no - int_no - exponent; 1200 1201 /* Construct the denominator. */ 1202 densize = 0; 1203 expbit = 1; 1204 do 1205 { 1206 if ((neg_exp & expbit) != 0) 1207 { 1208 mp_limb_t cy; 1209 neg_exp ^= expbit; 1210 1211 if (densize == 0) 1212 { 1213 densize = ttab->arraysize - _FPIO_CONST_OFFSET; 1214 memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET], 1215 densize * sizeof (mp_limb_t)); 1216 } 1217 else 1218 { 1219 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff 1220 + _FPIO_CONST_OFFSET], 1221 ttab->arraysize - _FPIO_CONST_OFFSET, 1222 psrc, densize); 1223 densize += ttab->arraysize - _FPIO_CONST_OFFSET; 1224 if (cy == 0) 1225 --densize; 1226 (void) SWAP (psrc, pdest); 1227 } 1228 } 1229 expbit <<= 1; 1230 ++ttab; 1231 } 1232 while (neg_exp != 0); 1233 1234 if (psrc == num) 1235 memcpy (den, num, densize * sizeof (mp_limb_t)); 1236 1237 /* Read the fractional digits from the string. */ 1238 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent 1239 #ifndef USE_WIDE_CHAR 1240 , decimal, decimal_len, thousands 1241 #endif 1242 ); 1243 1244 /* We now have to shift both numbers so that the highest bit in the 1245 denominator is set. In the same process we copy the numerator to 1246 a high place in the array so that the division constructs the wanted 1247 digits. This is done by a "quasi fix point" number representation. 1248 1249 num: ddddddddddd . 0000000000000000000000 1250 |--- m ---| 1251 den: ddddddddddd n >= m 1252 |--- n ---| 1253 */ 1254 1255 count_leading_zeros (cnt, den[densize - 1]); 1256 1257 if (cnt > 0) 1258 { 1259 /* Don't call `mpn_shift' with a count of zero since the specification 1260 does not allow this. */ 1261 (void) __mpn_lshift (den, den, densize, cnt); 1262 cy = __mpn_lshift (num, num, numsize, cnt); 1263 if (cy != 0) 1264 num[numsize++] = cy; 1265 } 1266 1267 /* Now we are ready for the division. But it is not necessary to 1268 do a full multi-precision division because we only need a small 1269 number of bits for the result. So we do not use __mpn_divmod 1270 here but instead do the division here by hand and stop whenever 1271 the needed number of bits is reached. The code itself comes 1272 from the GNU MP Library by Torbj\"orn Granlund. */ 1273 1274 exponent = bits; 1275 1276 switch (densize) 1277 { 1278 case 1: 1279 { 1280 mp_limb_t d, n, quot; 1281 int used = 0; 1282 1283 n = num[0]; 1284 d = den[0]; 1285 assert (numsize == 1 && n < d); 1286 1287 do 1288 { 1289 udiv_qrnnd (quot, n, n, 0, d); 1290 1291 #define got_limb \ 1292 if (bits == 0) \ 1293 { \ 1294 register int cnt; \ 1295 if (quot == 0) \ 1296 cnt = BITS_PER_MP_LIMB; \ 1297 else \ 1298 count_leading_zeros (cnt, quot); \ 1299 exponent -= cnt; \ 1300 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \ 1301 { \ 1302 used = MANT_DIG + cnt; \ 1303 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \ 1304 bits = MANT_DIG + 1; \ 1305 } \ 1306 else \ 1307 { \ 1308 /* Note that we only clear the second element. */ \ 1309 /* The conditional is determined at compile time. */ \ 1310 if (RETURN_LIMB_SIZE > 1) \ 1311 retval[1] = 0; \ 1312 retval[0] = quot; \ 1313 bits = -cnt; \ 1314 } \ 1315 } \ 1316 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \ 1317 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \ 1318 quot); \ 1319 else \ 1320 { \ 1321 used = MANT_DIG - bits; \ 1322 if (used > 0) \ 1323 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \ 1324 } \ 1325 bits += BITS_PER_MP_LIMB 1326 1327 got_limb; 1328 } 1329 while (bits <= MANT_DIG); 1330 1331 return round_and_return (retval, exponent - 1, negative, 1332 quot, BITS_PER_MP_LIMB - 1 - used, 1333 more_bits || n != 0); 1334 } 1335 case 2: 1336 { 1337 mp_limb_t d0, d1, n0, n1; 1338 mp_limb_t quot = 0; 1339 int used = 0; 1340 1341 d0 = den[0]; 1342 d1 = den[1]; 1343 1344 if (numsize < densize) 1345 { 1346 if (num[0] >= d1) 1347 { 1348 /* The numerator of the number occupies fewer bits than 1349 the denominator but the one limb is bigger than the 1350 high limb of the numerator. */ 1351 n1 = 0; 1352 n0 = num[0]; 1353 } 1354 else 1355 { 1356 if (bits <= 0) 1357 exponent -= BITS_PER_MP_LIMB; 1358 else 1359 { 1360 if (bits + BITS_PER_MP_LIMB <= MANT_DIG) 1361 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, 1362 BITS_PER_MP_LIMB, 0); 1363 else 1364 { 1365 used = MANT_DIG - bits; 1366 if (used > 0) 1367 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); 1368 } 1369 bits += BITS_PER_MP_LIMB; 1370 } 1371 n1 = num[0]; 1372 n0 = 0; 1373 } 1374 } 1375 else 1376 { 1377 n1 = num[1]; 1378 n0 = num[0]; 1379 } 1380 1381 while (bits <= MANT_DIG) 1382 { 1383 mp_limb_t r; 1384 1385 if (n1 == d1) 1386 { 1387 /* QUOT should be either 111..111 or 111..110. We need 1388 special treatment of this rare case as normal division 1389 would give overflow. */ 1390 quot = ~(mp_limb_t) 0; 1391 1392 r = n0 + d1; 1393 if (r < d1) /* Carry in the addition? */ 1394 { 1395 add_ssaaaa (n1, n0, r - d0, 0, 0, d0); 1396 goto have_quot; 1397 } 1398 n1 = d0 - (d0 != 0); 1399 n0 = -d0; 1400 } 1401 else 1402 { 1403 udiv_qrnnd (quot, r, n1, n0, d1); 1404 umul_ppmm (n1, n0, d0, quot); 1405 } 1406 1407 q_test: 1408 if (n1 > r || (n1 == r && n0 > 0)) 1409 { 1410 /* The estimated QUOT was too large. */ 1411 --quot; 1412 1413 sub_ddmmss (n1, n0, n1, n0, 0, d0); 1414 r += d1; 1415 if (r >= d1) /* If not carry, test QUOT again. */ 1416 goto q_test; 1417 } 1418 sub_ddmmss (n1, n0, r, 0, n1, n0); 1419 1420 have_quot: 1421 got_limb; 1422 } 1423 1424 return round_and_return (retval, exponent - 1, negative, 1425 quot, BITS_PER_MP_LIMB - 1 - used, 1426 more_bits || n1 != 0 || n0 != 0); 1427 } 1428 default: 1429 { 1430 int i; 1431 mp_limb_t cy, dX, d1, n0, n1; 1432 mp_limb_t quot = 0; 1433 int used = 0; 1434 1435 dX = den[densize - 1]; 1436 d1 = den[densize - 2]; 1437 1438 /* The division does not work if the upper limb of the two-limb 1439 numerator is greater than the denominator. */ 1440 if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0) 1441 num[numsize++] = 0; 1442 1443 if (numsize < densize) 1444 { 1445 mp_size_t empty = densize - numsize; 1446 1447 if (bits <= 0) 1448 { 1449 register int i; 1450 for (i = numsize; i > 0; --i) 1451 num[i + empty] = num[i - 1]; 1452 MPN_ZERO (num, empty + 1); 1453 exponent -= empty * BITS_PER_MP_LIMB; 1454 } 1455 else 1456 { 1457 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG) 1458 { 1459 /* We make a difference here because the compiler 1460 cannot optimize the `else' case that good and 1461 this reflects all currently used FLOAT types 1462 and GMP implementations. */ 1463 register int i; 1464 #if RETURN_LIMB_SIZE <= 2 1465 assert (empty == 1); 1466 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, 1467 BITS_PER_MP_LIMB, 0); 1468 #else 1469 for (i = RETURN_LIMB_SIZE; i > empty; --i) 1470 retval[i] = retval[i - empty]; 1471 #endif 1472 #if RETURN_LIMB_SIZE > 1 1473 retval[1] = 0; 1474 #endif 1475 for (i = numsize; i > 0; --i) 1476 num[i + empty] = num[i - 1]; 1477 MPN_ZERO (num, empty + 1); 1478 } 1479 else 1480 { 1481 used = MANT_DIG - bits; 1482 if (used >= BITS_PER_MP_LIMB) 1483 { 1484 register int i; 1485 (void) __mpn_lshift (&retval[used 1486 / BITS_PER_MP_LIMB], 1487 retval, RETURN_LIMB_SIZE, 1488 used % BITS_PER_MP_LIMB); 1489 for (i = used / BITS_PER_MP_LIMB; i >= 0; --i) 1490 retval[i] = 0; 1491 } 1492 else if (used > 0) 1493 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); 1494 } 1495 bits += empty * BITS_PER_MP_LIMB; 1496 } 1497 } 1498 else 1499 { 1500 int i; 1501 assert (numsize == densize); 1502 for (i = numsize; i > 0; --i) 1503 num[i] = num[i - 1]; 1504 } 1505 1506 den[densize] = 0; 1507 n0 = num[densize]; 1508 1509 while (bits <= MANT_DIG) 1510 { 1511 if (n0 == dX) 1512 /* This might over-estimate QUOT, but it's probably not 1513 worth the extra code here to find out. */ 1514 quot = ~(mp_limb_t) 0; 1515 else 1516 { 1517 mp_limb_t r; 1518 1519 udiv_qrnnd (quot, r, n0, num[densize - 1], dX); 1520 umul_ppmm (n1, n0, d1, quot); 1521 1522 while (n1 > r || (n1 == r && n0 > num[densize - 2])) 1523 { 1524 --quot; 1525 r += dX; 1526 if (r < dX) /* I.e. "carry in previous addition?" */ 1527 break; 1528 n1 -= n0 < d1; 1529 n0 -= d1; 1530 } 1531 } 1532 1533 /* Possible optimization: We already have (q * n0) and (1 * n1) 1534 after the calculation of QUOT. Taking advantage of this, we 1535 could make this loop make two iterations less. */ 1536 1537 cy = __mpn_submul_1 (num, den, densize + 1, quot); 1538 1539 if (num[densize] != cy) 1540 { 1541 cy = __mpn_add_n (num, num, den, densize); 1542 assert (cy != 0); 1543 --quot; 1544 } 1545 n0 = num[densize] = num[densize - 1]; 1546 for (i = densize - 1; i > 0; --i) 1547 num[i] = num[i - 1]; 1548 1549 got_limb; 1550 } 1551 1552 for (i = densize; num[i] == 0 && i >= 0; --i) 1553 ; 1554 return round_and_return (retval, exponent - 1, negative, 1555 quot, BITS_PER_MP_LIMB - 1 - used, 1556 more_bits || i >= 0); 1557 } 1558 } 1559 } 1560 1561 /* NOTREACHED */ 1562 } 1563 1564 /* External user entry point. */ 1565 1566 FLOAT 1567 #ifdef weak_function 1568 weak_function 1569 #endif 1570 STRTOF (nptr, endptr LOCALE_PARAM) 1571 const STRING_TYPE *nptr; 1572 STRING_TYPE **endptr; 1573 LOCALE_PARAM_DECL 1574 { 1575 return INTERNAL (STRTOF) (nptr, endptr, 0 LOCALE_PARAM); 1576 } 1577