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