stdlib/strtod_l.c

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