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