Back to index

glibc  2.9
strtod_l.c
Go to the documentation of this file.
00001 /* Convert string representing a number to float value, using given locale.
00002    Copyright (C) 1997,1998,2002,2004,2005,2006,2007,2008
00003    Free Software Foundation, Inc.
00004    This file is part of the GNU C Library.
00005    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
00006 
00007    The GNU C Library is free software; you can redistribute it and/or
00008    modify it under the terms of the GNU Lesser General Public
00009    License as published by the Free Software Foundation; either
00010    version 2.1 of the License, or (at your option) any later version.
00011 
00012    The GNU C Library is distributed in the hope that it will be useful,
00013    but WITHOUT ANY WARRANTY; without even the implied warranty of
00014    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015    Lesser General Public License for more details.
00016 
00017    You should have received a copy of the GNU Lesser General Public
00018    License along with the GNU C Library; if not, write to the Free
00019    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00020    02111-1307 USA.  */
00021 
00022 #include <xlocale.h>
00023 
00024 extern double ____strtod_l_internal (const char *, char **, int, __locale_t);
00025 extern unsigned long long int ____strtoull_l_internal (const char *, char **,
00026                                                  int, int, __locale_t);
00027 
00028 /* Configuration part.  These macros are defined by `strtold.c',
00029    `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
00030    `long double' and `float' versions of the reader.  */
00031 #ifndef FLOAT
00032 # include <math_ldbl_opt.h>
00033 # define FLOAT              double
00034 # define FLT         DBL
00035 # ifdef USE_WIDE_CHAR
00036 #  define STRTOF     wcstod_l
00037 #  define __STRTOF   __wcstod_l
00038 # else
00039 #  define STRTOF     strtod_l
00040 #  define __STRTOF   __strtod_l
00041 # endif
00042 # define MPN2FLOAT   __mpn_construct_double
00043 # define FLOAT_HUGE_VAL     HUGE_VAL
00044 # define SET_MANTISSA(flt, mant) \
00045   do { union ieee754_double u;                                              \
00046        u.d = (flt);                                                  \
00047        if ((mant & 0xfffffffffffffULL) == 0)                                \
00048         mant = 0x8000000000000ULL;                                   \
00049        u.ieee.mantissa0 = ((mant) >> 32) & 0xfffff;                         \
00050        u.ieee.mantissa1 = (mant) & 0xffffffff;                              \
00051        (flt) = u.d;                                                  \
00052   } while (0)
00053 #endif
00054 /* End of configuration part.  */
00055 
00056 #include <ctype.h>
00057 #include <errno.h>
00058 #include <float.h>
00059 #include <ieee754.h>
00060 #include "../locale/localeinfo.h"
00061 #include <locale.h>
00062 #include <math.h>
00063 #include <stdlib.h>
00064 #include <string.h>
00065 
00066 /* The gmp headers need some configuration frobs.  */
00067 #define HAVE_ALLOCA 1
00068 
00069 /* Include gmp-mparam.h first, such that definitions of _SHORT_LIMB
00070    and _LONG_LONG_LIMB in it can take effect into gmp.h.  */
00071 #include <gmp-mparam.h>
00072 #include <gmp.h>
00073 #include "gmp-impl.h"
00074 #include "longlong.h"
00075 #include "fpioconst.h"
00076 
00077 #define NDEBUG 1
00078 #include <assert.h>
00079 
00080 
00081 /* We use this code for the extended locale handling where the
00082    function gets as an additional argument the locale which has to be
00083    used.  To access the values we have to redefine the _NL_CURRENT and
00084    _NL_CURRENT_WORD macros.  */
00085 #undef _NL_CURRENT
00086 #define _NL_CURRENT(category, item) \
00087   (current->values[_NL_ITEM_INDEX (item)].string)
00088 #undef _NL_CURRENT_WORD
00089 #define _NL_CURRENT_WORD(category, item) \
00090   ((uint32_t) current->values[_NL_ITEM_INDEX (item)].word)
00091 
00092 #if defined _LIBC || defined HAVE_WCHAR_H
00093 # include <wchar.h>
00094 #endif
00095 
00096 #ifdef USE_WIDE_CHAR
00097 # include <wctype.h>
00098 # define STRING_TYPE wchar_t
00099 # define CHAR_TYPE wint_t
00100 # define L_(Ch) L##Ch
00101 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
00102 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
00103 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
00104 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
00105 # define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
00106 # define STRNCASECMP(S1, S2, N) \
00107   __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
00108 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
00109 #else
00110 # define STRING_TYPE char
00111 # define CHAR_TYPE char
00112 # define L_(Ch) Ch
00113 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
00114 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
00115 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
00116 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
00117 # define TOLOWER_C(Ch) __tolower_l ((Ch), _nl_C_locobj_ptr)
00118 # define STRNCASECMP(S1, S2, N) \
00119   __strncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
00120 # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc)
00121 #endif
00122 
00123 
00124 /* Constants we need from float.h; select the set for the FLOAT precision.  */
00125 #define MANT_DIG     PASTE(FLT,_MANT_DIG)
00126 #define       DIG           PASTE(FLT,_DIG)
00127 #define       MAX_EXP              PASTE(FLT,_MAX_EXP)
00128 #define       MIN_EXP              PASTE(FLT,_MIN_EXP)
00129 #define MAX_10_EXP   PASTE(FLT,_MAX_10_EXP)
00130 #define MIN_10_EXP   PASTE(FLT,_MIN_10_EXP)
00131 
00132 /* Extra macros required to get FLT expanded before the pasting.  */
00133 #define PASTE(a,b)   PASTE1(a,b)
00134 #define PASTE1(a,b)  a##b
00135 
00136 /* Function to construct a floating point number from an MP integer
00137    containing the fraction bits, a base 2 exponent, and a sign flag.  */
00138 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
00139 
00140 /* Definitions according to limb size used.  */
00141 #if    BITS_PER_MP_LIMB == 32
00142 # define MAX_DIG_PER_LIMB   9
00143 # define MAX_FAC_PER_LIMB   1000000000UL
00144 #elif  BITS_PER_MP_LIMB == 64
00145 # define MAX_DIG_PER_LIMB   19
00146 # define MAX_FAC_PER_LIMB   10000000000000000000ULL
00147 #else
00148 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
00149 #endif
00150 
00151 extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1];
00152 
00153 #ifndef       howmany
00154 #define       howmany(x,y)         (((x)+((y)-1))/(y))
00155 #endif
00156 #define SWAP(x, y)          ({ typeof(x) _tmp = x; x = y; y = _tmp; })
00157 
00158 #define NDIG                (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
00159 #define HEXNDIG                    ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
00160 #define       RETURN_LIMB_SIZE            howmany (MANT_DIG, BITS_PER_MP_LIMB)
00161 
00162 #define RETURN(val,end)                                                     \
00163     do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end);               \
00164         return val; } while (0)
00165 
00166 /* Maximum size necessary for mpn integers to hold floating point numbers.  */
00167 #define       MPNSIZE              (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
00168                       + 2)
00169 /* Declare an mpn integer variable that big.  */
00170 #define       MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
00171 /* Copy an mpn integer value.  */
00172 #define MPN_ASSIGN(dst, src) \
00173        memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
00174 
00175 
00176 /* Return a floating point number of the needed type according to the given
00177    multi-precision number after possible rounding.  */
00178 static FLOAT
00179 round_and_return (mp_limb_t *retval, int exponent, int negative,
00180                 mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
00181 {
00182   if (exponent < MIN_EXP - 1)
00183     {
00184       mp_size_t shift = MIN_EXP - 1 - exponent;
00185 
00186       if (shift > MANT_DIG)
00187        {
00188          __set_errno (EDOM);
00189          return 0.0;
00190        }
00191 
00192       more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
00193       if (shift == MANT_DIG)
00194        /* This is a special case to handle the very seldom case where
00195           the mantissa will be empty after the shift.  */
00196        {
00197          int i;
00198 
00199          round_limb = retval[RETURN_LIMB_SIZE - 1];
00200          round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
00201          for (i = 0; i < RETURN_LIMB_SIZE; ++i)
00202            more_bits |= retval[i] != 0;
00203          MPN_ZERO (retval, RETURN_LIMB_SIZE);
00204        }
00205       else if (shift >= BITS_PER_MP_LIMB)
00206        {
00207          int i;
00208 
00209          round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
00210          round_bit = (shift - 1) % BITS_PER_MP_LIMB;
00211          for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
00212            more_bits |= retval[i] != 0;
00213          more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
00214                      != 0);
00215 
00216          (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
00217                                RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
00218                                shift % BITS_PER_MP_LIMB);
00219           MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
00220                     shift / BITS_PER_MP_LIMB);
00221        }
00222       else if (shift > 0)
00223        {
00224           round_limb = retval[0];
00225           round_bit = shift - 1;
00226          (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
00227        }
00228       /* This is a hook for the m68k long double format, where the
00229         exponent bias is the same for normalized and denormalized
00230         numbers.  */
00231 #ifndef DENORM_EXP
00232 # define DENORM_EXP (MIN_EXP - 2)
00233 #endif
00234       exponent = DENORM_EXP;
00235     }
00236 
00237   if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
00238       && (more_bits || (retval[0] & 1) != 0
00239           || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
00240     {
00241       mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
00242 
00243       if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
00244           ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
00245            (retval[RETURN_LIMB_SIZE - 1]
00246             & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
00247        {
00248          ++exponent;
00249          (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
00250          retval[RETURN_LIMB_SIZE - 1]
00251            |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
00252        }
00253       else if (exponent == DENORM_EXP
00254               && (retval[RETURN_LIMB_SIZE - 1]
00255                  & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
00256               != 0)
00257          /* The number was denormalized but now normalized.  */
00258        exponent = MIN_EXP - 1;
00259     }
00260 
00261   if (exponent > MAX_EXP)
00262     return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
00263 
00264   return MPN2FLOAT (retval, exponent, negative);
00265 }
00266 
00267 
00268 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
00269    into N.  Return the size of the number limbs in NSIZE at the first
00270    character od the string that is not part of the integer as the function
00271    value.  If the EXPONENT is small enough to be taken as an additional
00272    factor for the resulting number (see code) multiply by it.  */
00273 static const STRING_TYPE *
00274 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
00275            int *exponent
00276 #ifndef USE_WIDE_CHAR
00277            , const char *decimal, size_t decimal_len, const char *thousands
00278 #endif
00279 
00280            )
00281 {
00282   /* Number of digits for actual limb.  */
00283   int cnt = 0;
00284   mp_limb_t low = 0;
00285   mp_limb_t start;
00286 
00287   *nsize = 0;
00288   assert (digcnt > 0);
00289   do
00290     {
00291       if (cnt == MAX_DIG_PER_LIMB)
00292        {
00293          if (*nsize == 0)
00294            {
00295              n[0] = low;
00296              *nsize = 1;
00297            }
00298          else
00299            {
00300              mp_limb_t cy;
00301              cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
00302              cy += __mpn_add_1 (n, n, *nsize, low);
00303              if (cy != 0)
00304               {
00305                 n[*nsize] = cy;
00306                 ++(*nsize);
00307               }
00308            }
00309          cnt = 0;
00310          low = 0;
00311        }
00312 
00313       /* There might be thousands separators or radix characters in
00314         the string.  But these all can be ignored because we know the
00315         format of the number is correct and we have an exact number
00316         of characters to read.  */
00317 #ifdef USE_WIDE_CHAR
00318       if (*str < L'0' || *str > L'9')
00319        ++str;
00320 #else
00321       if (*str < '0' || *str > '9')
00322        {
00323          int inner = 0;
00324          if (thousands != NULL && *str == *thousands
00325              && ({ for (inner = 1; thousands[inner] != '\0'; ++inner)
00326                     if (thousands[inner] != str[inner])
00327                      break;
00328                   thousands[inner] == '\0'; }))
00329            str += inner;
00330          else
00331            str += decimal_len;
00332        }
00333 #endif
00334       low = low * 10 + *str++ - L_('0');
00335       ++cnt;
00336     }
00337   while (--digcnt > 0);
00338 
00339   if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
00340     {
00341       low *= _tens_in_limb[*exponent];
00342       start = _tens_in_limb[cnt + *exponent];
00343       *exponent = 0;
00344     }
00345   else
00346     start = _tens_in_limb[cnt];
00347 
00348   if (*nsize == 0)
00349     {
00350       n[0] = low;
00351       *nsize = 1;
00352     }
00353   else
00354     {
00355       mp_limb_t cy;
00356       cy = __mpn_mul_1 (n, n, *nsize, start);
00357       cy += __mpn_add_1 (n, n, *nsize, low);
00358       if (cy != 0)
00359        n[(*nsize)++] = cy;
00360     }
00361 
00362   return str;
00363 }
00364 
00365 
00366 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
00367    with the COUNT most significant bits of LIMB.
00368 
00369    Tege doesn't like this function so I have to write it here myself. :)
00370    --drepper */
00371 static inline void
00372 __attribute ((always_inline))
00373 __mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
00374               mp_limb_t limb)
00375 {
00376   if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB)
00377     {
00378       /* Optimize the case of shifting by exactly a word:
00379         just copy words, with no actual bit-shifting.  */
00380       mp_size_t i;
00381       for (i = size - 1; i > 0; --i)
00382        ptr[i] = ptr[i - 1];
00383       ptr[0] = limb;
00384     }
00385   else
00386     {
00387       (void) __mpn_lshift (ptr, ptr, size, count);
00388       ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
00389     }
00390 }
00391 
00392 
00393 #define INTERNAL(x) INTERNAL1(x)
00394 #define INTERNAL1(x) __##x##_internal
00395 #ifndef ____STRTOF_INTERNAL
00396 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
00397 #endif
00398 
00399 /* This file defines a function to check for correct grouping.  */
00400 #include "grouping.h"
00401 
00402 
00403 /* Return a floating point number with the value of the given string NPTR.
00404    Set *ENDPTR to the character after the last used one.  If the number is
00405    smaller than the smallest representable number, set `errno' to ERANGE and
00406    return 0.0.  If the number is too big to be represented, set `errno' to
00407    ERANGE and return HUGE_VAL with the appropriate sign.  */
00408 FLOAT
00409 ____STRTOF_INTERNAL (nptr, endptr, group, loc)
00410      const STRING_TYPE *nptr;
00411      STRING_TYPE **endptr;
00412      int group;
00413      __locale_t loc;
00414 {
00415   int negative;                    /* The sign of the number.  */
00416   MPN_VAR (num);            /* MP representation of the number.  */
00417   int exponent;                    /* Exponent of the number.  */
00418 
00419   /* Numbers starting `0X' or `0x' have to be processed with base 16.  */
00420   int base = 10;
00421 
00422   /* When we have to compute fractional digits we form a fraction with a
00423      second multi-precision number (and we sometimes need a second for
00424      temporary results).  */
00425   MPN_VAR (den);
00426 
00427   /* Representation for the return value.  */
00428   mp_limb_t retval[RETURN_LIMB_SIZE];
00429   /* Number of bits currently in result value.  */
00430   int bits;
00431 
00432   /* Running pointer after the last character processed in the string.  */
00433   const STRING_TYPE *cp, *tp;
00434   /* Start of significant part of the number.  */
00435   const STRING_TYPE *startp, *start_of_digits;
00436   /* Points at the character following the integer and fractional digits.  */
00437   const STRING_TYPE *expp;
00438   /* Total number of digit and number of digits in integer part.  */
00439   int dig_no, int_no, lead_zero;
00440   /* Contains the last character read.  */
00441   CHAR_TYPE c;
00442 
00443 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
00444    there.  So define it ourselves if it remains undefined.  */
00445 #ifndef _WINT_T
00446   typedef unsigned int wint_t;
00447 #endif
00448   /* The radix character of the current locale.  */
00449 #ifdef USE_WIDE_CHAR
00450   wchar_t decimal;
00451 #else
00452   const char *decimal;
00453   size_t decimal_len;
00454 #endif
00455   /* The thousands character of the current locale.  */
00456 #ifdef USE_WIDE_CHAR
00457   wchar_t thousands = L'\0';
00458 #else
00459   const char *thousands = NULL;
00460 #endif
00461   /* The numeric grouping specification of the current locale,
00462      in the format described in <locale.h>.  */
00463   const char *grouping;
00464   /* Used in several places.  */
00465   int cnt;
00466 
00467   struct locale_data *current = loc->__locales[LC_NUMERIC];
00468 
00469   if (__builtin_expect (group, 0))
00470     {
00471       grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
00472       if (*grouping <= 0 || *grouping == CHAR_MAX)
00473        grouping = NULL;
00474       else
00475        {
00476          /* Figure out the thousands separator character.  */
00477 #ifdef USE_WIDE_CHAR
00478          thousands = _NL_CURRENT_WORD (LC_NUMERIC,
00479                                    _NL_NUMERIC_THOUSANDS_SEP_WC);
00480          if (thousands == L'\0')
00481            grouping = NULL;
00482 #else
00483          thousands = _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
00484          if (*thousands == '\0')
00485            {
00486              thousands = NULL;
00487              grouping = NULL;
00488            }
00489 #endif
00490        }
00491     }
00492   else
00493     grouping = NULL;
00494 
00495   /* Find the locale's decimal point character.  */
00496 #ifdef USE_WIDE_CHAR
00497   decimal = _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_DECIMAL_POINT_WC);
00498   assert (decimal != L'\0');
00499 # define decimal_len 1
00500 #else
00501   decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
00502   decimal_len = strlen (decimal);
00503   assert (decimal_len > 0);
00504 #endif
00505 
00506   /* Prepare number representation.  */
00507   exponent = 0;
00508   negative = 0;
00509   bits = 0;
00510 
00511   /* Parse string to get maximal legal prefix.  We need the number of
00512      characters of the integer part, the fractional part and the exponent.  */
00513   cp = nptr - 1;
00514   /* Ignore leading white space.  */
00515   do
00516     c = *++cp;
00517   while (ISSPACE (c));
00518 
00519   /* Get sign of the result.  */
00520   if (c == L_('-'))
00521     {
00522       negative = 1;
00523       c = *++cp;
00524     }
00525   else if (c == L_('+'))
00526     c = *++cp;
00527 
00528   /* Return 0.0 if no legal string is found.
00529      No character is used even if a sign was found.  */
00530 #ifdef USE_WIDE_CHAR
00531   if (c == (wint_t) decimal
00532       && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9')
00533     {
00534       /* We accept it.  This funny construct is here only to indent
00535         the code correctly.  */
00536     }
00537 #else
00538   for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
00539     if (cp[cnt] != decimal[cnt])
00540       break;
00541   if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9')
00542     {
00543       /* We accept it.  This funny construct is here only to indent
00544         the code correctly.  */
00545     }
00546 #endif
00547   else if (c < L_('0') || c > L_('9'))
00548     {
00549       /* Check for `INF' or `INFINITY'.  */
00550       CHAR_TYPE lowc = TOLOWER_C (c);
00551 
00552       if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0)
00553        {
00554          /* Return +/- infinity.  */
00555          if (endptr != NULL)
00556            *endptr = (STRING_TYPE *)
00557                     (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0
00558                           ? 8 : 3));
00559 
00560          return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
00561        }
00562 
00563       if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
00564        {
00565          /* Return NaN.  */
00566          FLOAT retval = NAN;
00567 
00568          cp += 3;
00569 
00570          /* Match `(n-char-sequence-digit)'.  */
00571          if (*cp == L_('('))
00572            {
00573              const STRING_TYPE *startp = cp;
00574              do
00575               ++cp;
00576              while ((*cp >= L_('0') && *cp <= L_('9'))
00577                    || ({ CHAR_TYPE lo = TOLOWER (*cp);
00578                         lo >= L_('a') && lo <= L_('z'); })
00579                    || *cp == L_('_'));
00580 
00581              if (*cp != L_(')'))
00582               /* The closing brace is missing.  Only match the NAN
00583                  part.  */
00584               cp = startp;
00585              else
00586               {
00587                 /* This is a system-dependent way to specify the
00588                    bitmask used for the NaN.  We expect it to be
00589                    a number which is put in the mantissa of the
00590                    number.  */
00591                 STRING_TYPE *endp;
00592                 unsigned long long int mant;
00593 
00594                 mant = STRTOULL (startp + 1, &endp, 0);
00595                 if (endp == cp)
00596                   SET_MANTISSA (retval, mant);
00597 
00598                 /* Consume the closing brace.  */
00599                 ++cp;
00600               }
00601            }
00602 
00603          if (endptr != NULL)
00604            *endptr = (STRING_TYPE *) cp;
00605 
00606          return retval;
00607        }
00608 
00609       /* It is really a text we do not recognize.  */
00610       RETURN (0.0, nptr);
00611     }
00612 
00613   /* First look whether we are faced with a hexadecimal number.  */
00614   if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
00615     {
00616       /* Okay, it is a hexa-decimal number.  Remember this and skip
00617         the characters.  BTW: hexadecimal numbers must not be
00618         grouped.  */
00619       base = 16;
00620       cp += 2;
00621       c = *cp;
00622       grouping = NULL;
00623     }
00624 
00625   /* Record the start of the digits, in case we will check their grouping.  */
00626   start_of_digits = startp = cp;
00627 
00628   /* Ignore leading zeroes.  This helps us to avoid useless computations.  */
00629 #ifdef USE_WIDE_CHAR
00630   while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands))
00631     c = *++cp;
00632 #else
00633   if (__builtin_expect (thousands == NULL, 1))
00634     while (c == '0')
00635       c = *++cp;
00636   else
00637     {
00638       /* We also have the multibyte thousands string.  */
00639       while (1)
00640        {
00641          if (c != '0')
00642            {
00643              for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
00644               if (thousands[cnt] != cp[cnt])
00645                 break;
00646              if (thousands[cnt] != '\0')
00647               break;
00648              cp += cnt - 1;
00649            }
00650          c = *++cp;
00651        }
00652     }
00653 #endif
00654 
00655   /* If no other digit but a '0' is found the result is 0.0.
00656      Return current read pointer.  */
00657   CHAR_TYPE lowc = TOLOWER (c);
00658   if (!((c >= L_('0') && c <= L_('9'))
00659        || (base == 16 && lowc >= L_('a') && lowc <= L_('f'))
00660        || (
00661 #ifdef USE_WIDE_CHAR
00662            c == (wint_t) decimal
00663 #else
00664            ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
00665                if (decimal[cnt] != cp[cnt])
00666                  break;
00667               decimal[cnt] == '\0'; })
00668 #endif
00669            /* '0x.' alone is not a valid hexadecimal number.
00670               '.' alone is not valid either, but that has been checked
00671               already earlier.  */
00672            && (base != 16
00673               || cp != start_of_digits
00674               || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9'))
00675               || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]);
00676                     lo >= L_('a') && lo <= L_('f'); })))
00677        || (base == 16 && (cp != start_of_digits
00678                         && lowc == L_('p')))
00679        || (base != 16 && lowc == L_('e'))))
00680     {
00681 #ifdef USE_WIDE_CHAR
00682       tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
00683                                     grouping);
00684 #else
00685       tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
00686                                     grouping);
00687 #endif
00688       /* If TP is at the start of the digits, there was no correctly
00689         grouped prefix of the string; so no number found.  */
00690       RETURN (negative ? -0.0 : 0.0,
00691              tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
00692     }
00693 
00694   /* Remember first significant digit and read following characters until the
00695      decimal point, exponent character or any non-FP number character.  */
00696   startp = cp;
00697   dig_no = 0;
00698   while (1)
00699     {
00700       if ((c >= L_('0') && c <= L_('9'))
00701          || (base == 16
00702              && ({ CHAR_TYPE lo = TOLOWER (c);
00703                   lo >= L_('a') && lo <= L_('f'); })))
00704        ++dig_no;
00705       else
00706        {
00707 #ifdef USE_WIDE_CHAR
00708          if (__builtin_expect ((wint_t) thousands == L'\0', 1)
00709              || c != (wint_t) thousands)
00710            /* Not a digit or separator: end of the integer part.  */
00711            break;
00712 #else
00713          if (__builtin_expect (thousands == NULL, 1))
00714            break;
00715          else
00716            {
00717              for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
00718               if (thousands[cnt] != cp[cnt])
00719                 break;
00720              if (thousands[cnt] != '\0')
00721               break;
00722              cp += cnt - 1;
00723            }
00724 #endif
00725        }
00726       c = *++cp;
00727     }
00728 
00729   if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits)
00730     {
00731       /* Check the grouping of the digits.  */
00732 #ifdef USE_WIDE_CHAR
00733       tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
00734                                     grouping);
00735 #else
00736       tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
00737                                     grouping);
00738 #endif
00739       if (cp != tp)
00740         {
00741          /* Less than the entire string was correctly grouped.  */
00742 
00743          if (tp == start_of_digits)
00744            /* No valid group of numbers at all: no valid number.  */
00745            RETURN (0.0, nptr);
00746 
00747          if (tp < startp)
00748            /* The number is validly grouped, but consists
00749               only of zeroes.  The whole value is zero.  */
00750            RETURN (negative ? -0.0 : 0.0, tp);
00751 
00752          /* Recompute DIG_NO so we won't read more digits than
00753             are properly grouped.  */
00754          cp = tp;
00755          dig_no = 0;
00756          for (tp = startp; tp < cp; ++tp)
00757            if (*tp >= L_('0') && *tp <= L_('9'))
00758              ++dig_no;
00759 
00760          int_no = dig_no;
00761          lead_zero = 0;
00762 
00763          goto number_parsed;
00764        }
00765     }
00766 
00767   /* We have the number of digits in the integer part.  Whether these
00768      are all or any is really a fractional digit will be decided
00769      later.  */
00770   int_no = dig_no;
00771   lead_zero = int_no == 0 ? -1 : 0;
00772 
00773   /* Read the fractional digits.  A special case are the 'american
00774      style' numbers like `16.' i.e. with decimal point but without
00775      trailing digits.  */
00776   if (
00777 #ifdef USE_WIDE_CHAR
00778       c == (wint_t) decimal
00779 #else
00780       ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
00781           if (decimal[cnt] != cp[cnt])
00782             break;
00783         decimal[cnt] == '\0'; })
00784 #endif
00785       )
00786     {
00787       cp += decimal_len;
00788       c = *cp;
00789       while ((c >= L_('0') && c <= L_('9')) ||
00790             (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c);
00791                             lo >= L_('a') && lo <= L_('f'); })))
00792        {
00793          if (c != L_('0') && lead_zero == -1)
00794            lead_zero = dig_no - int_no;
00795          ++dig_no;
00796          c = *++cp;
00797        }
00798     }
00799 
00800   /* Remember start of exponent (if any).  */
00801   expp = cp;
00802 
00803   /* Read exponent.  */
00804   lowc = TOLOWER (c);
00805   if ((base == 16 && lowc == L_('p'))
00806       || (base != 16 && lowc == L_('e')))
00807     {
00808       int exp_negative = 0;
00809 
00810       c = *++cp;
00811       if (c == L_('-'))
00812        {
00813          exp_negative = 1;
00814          c = *++cp;
00815        }
00816       else if (c == L_('+'))
00817        c = *++cp;
00818 
00819       if (c >= L_('0') && c <= L_('9'))
00820        {
00821          int exp_limit;
00822 
00823          /* Get the exponent limit. */
00824          if (base == 16)
00825            exp_limit = (exp_negative ?
00826                       -MIN_EXP + MANT_DIG + 4 * int_no :
00827                       MAX_EXP - 4 * int_no + 4 * lead_zero + 3);
00828          else
00829            exp_limit = (exp_negative ?
00830                       -MIN_10_EXP + MANT_DIG + int_no :
00831                       MAX_10_EXP - int_no + lead_zero + 1);
00832 
00833          do
00834            {
00835              exponent *= 10;
00836              exponent += c - L_('0');
00837 
00838              if (__builtin_expect (exponent > exp_limit, 0))
00839               /* The exponent is too large/small to represent a valid
00840                  number.  */
00841               {
00842                 FLOAT result;
00843 
00844                 /* We have to take care for special situation: a joker
00845                    might have written "0.0e100000" which is in fact
00846                    zero.  */
00847                 if (lead_zero == -1)
00848                   result = negative ? -0.0 : 0.0;
00849                 else
00850                   {
00851                     /* Overflow or underflow.  */
00852                     __set_errno (ERANGE);
00853                     result = (exp_negative ? (negative ? -0.0 : 0.0) :
00854                             negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
00855                   }
00856 
00857                 /* Accept all following digits as part of the exponent.  */
00858                 do
00859                   ++cp;
00860                 while (*cp >= L_('0') && *cp <= L_('9'));
00861 
00862                 RETURN (result, cp);
00863                 /* NOTREACHED */
00864               }
00865 
00866              c = *++cp;
00867            }
00868          while (c >= L_('0') && c <= L_('9'));
00869 
00870          if (exp_negative)
00871            exponent = -exponent;
00872        }
00873       else
00874        cp = expp;
00875     }
00876 
00877   /* We don't want to have to work with trailing zeroes after the radix.  */
00878   if (dig_no > int_no)
00879     {
00880       while (expp[-1] == L_('0'))
00881        {
00882          --expp;
00883          --dig_no;
00884        }
00885       assert (dig_no >= int_no);
00886     }
00887 
00888   if (dig_no == int_no && dig_no > 0 && exponent < 0)
00889     do
00890       {
00891        while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1])))
00892          --expp;
00893 
00894        if (expp[-1] != L_('0'))
00895          break;
00896 
00897        --expp;
00898        --dig_no;
00899        --int_no;
00900        exponent += base == 16 ? 4 : 1;
00901       }
00902     while (dig_no > 0 && exponent < 0);
00903 
00904  number_parsed:
00905 
00906   /* The whole string is parsed.  Store the address of the next character.  */
00907   if (endptr)
00908     *endptr = (STRING_TYPE *) cp;
00909 
00910   if (dig_no == 0)
00911     return negative ? -0.0 : 0.0;
00912 
00913   if (lead_zero)
00914     {
00915       /* Find the decimal point */
00916 #ifdef USE_WIDE_CHAR
00917       while (*startp != decimal)
00918        ++startp;
00919 #else
00920       while (1)
00921        {
00922          if (*startp == decimal[0])
00923            {
00924              for (cnt = 1; decimal[cnt] != '\0'; ++cnt)
00925               if (decimal[cnt] != startp[cnt])
00926                 break;
00927              if (decimal[cnt] == '\0')
00928               break;
00929            }
00930          ++startp;
00931        }
00932 #endif
00933       startp += lead_zero + decimal_len;
00934       exponent -= base == 16 ? 4 * lead_zero : lead_zero;
00935       dig_no -= lead_zero;
00936     }
00937 
00938   /* If the BASE is 16 we can use a simpler algorithm.  */
00939   if (base == 16)
00940     {
00941       static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
00942                                  4, 4, 4, 4, 4, 4, 4, 4 };
00943       int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
00944       int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
00945       mp_limb_t val;
00946 
00947       while (!ISXDIGIT (*startp))
00948        ++startp;
00949       while (*startp == L_('0'))
00950        ++startp;
00951       if (ISDIGIT (*startp))
00952        val = *startp++ - L_('0');
00953       else
00954        val = 10 + TOLOWER (*startp++) - L_('a');
00955       bits = nbits[val];
00956       /* We cannot have a leading zero.  */
00957       assert (bits != 0);
00958 
00959       if (pos + 1 >= 4 || pos + 1 >= bits)
00960        {
00961          /* We don't have to care for wrapping.  This is the normal
00962             case so we add the first clause in the `if' expression as
00963             an optimization.  It is a compile-time constant and so does
00964             not cost anything.  */
00965          retval[idx] = val << (pos - bits + 1);
00966          pos -= bits;
00967        }
00968       else
00969        {
00970          retval[idx--] = val >> (bits - pos - 1);
00971          retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
00972          pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
00973        }
00974 
00975       /* Adjust the exponent for the bits we are shifting in.  */
00976       exponent += bits - 1 + (int_no - 1) * 4;
00977 
00978       while (--dig_no > 0 && idx >= 0)
00979        {
00980          if (!ISXDIGIT (*startp))
00981            startp += decimal_len;
00982          if (ISDIGIT (*startp))
00983            val = *startp++ - L_('0');
00984          else
00985            val = 10 + TOLOWER (*startp++) - L_('a');
00986 
00987          if (pos + 1 >= 4)
00988            {
00989              retval[idx] |= val << (pos - 4 + 1);
00990              pos -= 4;
00991            }
00992          else
00993            {
00994              retval[idx--] |= val >> (4 - pos - 1);
00995              val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
00996              if (idx < 0)
00997               return round_and_return (retval, exponent, negative, val,
00998                                     BITS_PER_MP_LIMB - 1, dig_no > 0);
00999 
01000              retval[idx] = val;
01001              pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
01002            }
01003        }
01004 
01005       /* We ran out of digits.  */
01006       MPN_ZERO (retval, idx);
01007 
01008       return round_and_return (retval, exponent, negative, 0, 0, 0);
01009     }
01010 
01011   /* Now we have the number of digits in total and the integer digits as well
01012      as the exponent and its sign.  We can decide whether the read digits are
01013      really integer digits or belong to the fractional part; i.e. we normalize
01014      123e-2 to 1.23.  */
01015   {
01016     register int incr = (exponent < 0 ? MAX (-int_no, exponent)
01017                       : MIN (dig_no - int_no, exponent));
01018     int_no += incr;
01019     exponent -= incr;
01020   }
01021 
01022   if (__builtin_expect (int_no + exponent > MAX_10_EXP + 1, 0))
01023     {
01024       __set_errno (ERANGE);
01025       return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
01026     }
01027 
01028   if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0))
01029     {
01030       __set_errno (ERANGE);
01031       return negative ? -0.0 : 0.0;
01032     }
01033 
01034   if (int_no > 0)
01035     {
01036       /* Read the integer part as a multi-precision number to NUM.  */
01037       startp = str_to_mpn (startp, int_no, num, &numsize, &exponent
01038 #ifndef USE_WIDE_CHAR
01039                         , decimal, decimal_len, thousands
01040 #endif
01041                         );
01042 
01043       if (exponent > 0)
01044        {
01045          /* We now multiply the gained number by the given power of ten.  */
01046          mp_limb_t *psrc = num;
01047          mp_limb_t *pdest = den;
01048          int expbit = 1;
01049          const struct mp_power *ttab = &_fpioconst_pow10[0];
01050 
01051          do
01052            {
01053              if ((exponent & expbit) != 0)
01054               {
01055                 size_t size = ttab->arraysize - _FPIO_CONST_OFFSET;
01056                 mp_limb_t cy;
01057                 exponent ^= expbit;
01058 
01059                 /* FIXME: not the whole multiplication has to be
01060                    done.  If we have the needed number of bits we
01061                    only need the information whether more non-zero
01062                    bits follow.  */
01063                 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
01064                   cy = __mpn_mul (pdest, psrc, numsize,
01065                                 &__tens[ttab->arrayoff
01066                                       + _FPIO_CONST_OFFSET],
01067                                 size);
01068                 else
01069                   cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
01070                                             + _FPIO_CONST_OFFSET],
01071                                 size, psrc, numsize);
01072                 numsize += size;
01073                 if (cy == 0)
01074                   --numsize;
01075                 (void) SWAP (psrc, pdest);
01076               }
01077              expbit <<= 1;
01078              ++ttab;
01079            }
01080          while (exponent != 0);
01081 
01082          if (psrc == den)
01083            memcpy (num, den, numsize * sizeof (mp_limb_t));
01084        }
01085 
01086       /* Determine how many bits of the result we already have.  */
01087       count_leading_zeros (bits, num[numsize - 1]);
01088       bits = numsize * BITS_PER_MP_LIMB - bits;
01089 
01090       /* Now we know the exponent of the number in base two.
01091         Check it against the maximum possible exponent.  */
01092       if (__builtin_expect (bits > MAX_EXP, 0))
01093        {
01094          __set_errno (ERANGE);
01095          return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
01096        }
01097 
01098       /* We have already the first BITS bits of the result.  Together with
01099         the information whether more non-zero bits follow this is enough
01100         to determine the result.  */
01101       if (bits > MANT_DIG)
01102        {
01103          int i;
01104          const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
01105          const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
01106          const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
01107                                                : least_idx;
01108          const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
01109                                                : least_bit - 1;
01110 
01111          if (least_bit == 0)
01112            memcpy (retval, &num[least_idx],
01113                   RETURN_LIMB_SIZE * sizeof (mp_limb_t));
01114          else
01115             {
01116               for (i = least_idx; i < numsize - 1; ++i)
01117                 retval[i - least_idx] = (num[i] >> least_bit)
01118                                         | (num[i + 1]
01119                                            << (BITS_PER_MP_LIMB - least_bit));
01120               if (i - least_idx < RETURN_LIMB_SIZE)
01121                 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
01122             }
01123 
01124          /* Check whether any limb beside the ones in RETVAL are non-zero.  */
01125          for (i = 0; num[i] == 0; ++i)
01126            ;
01127 
01128          return round_and_return (retval, bits - 1, negative,
01129                                num[round_idx], round_bit,
01130                                int_no < dig_no || i < round_idx);
01131          /* NOTREACHED */
01132        }
01133       else if (dig_no == int_no)
01134        {
01135          const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
01136          const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
01137 
01138          if (target_bit == is_bit)
01139            {
01140              memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
01141                     numsize * sizeof (mp_limb_t));
01142              /* FIXME: the following loop can be avoided if we assume a
01143                maximal MANT_DIG value.  */
01144              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
01145            }
01146          else if (target_bit > is_bit)
01147            {
01148              (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
01149                                num, numsize, target_bit - is_bit);
01150              /* FIXME: the following loop can be avoided if we assume a
01151                maximal MANT_DIG value.  */
01152              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
01153            }
01154          else
01155            {
01156              mp_limb_t cy;
01157              assert (numsize < RETURN_LIMB_SIZE);
01158 
01159              cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
01160                              num, numsize, is_bit - target_bit);
01161              retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
01162              /* FIXME: the following loop can be avoided if we assume a
01163                maximal MANT_DIG value.  */
01164              MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
01165            }
01166 
01167          return round_and_return (retval, bits - 1, negative, 0, 0, 0);
01168          /* NOTREACHED */
01169        }
01170 
01171       /* Store the bits we already have.  */
01172       memcpy (retval, num, numsize * sizeof (mp_limb_t));
01173 #if RETURN_LIMB_SIZE > 1
01174       if (numsize < RETURN_LIMB_SIZE)
01175 # if RETURN_LIMB_SIZE == 2
01176         retval[numsize] = 0;
01177 # else
01178        MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize);
01179 # endif
01180 #endif
01181     }
01182 
01183   /* We have to compute at least some of the fractional digits.  */
01184   {
01185     /* We construct a fraction and the result of the division gives us
01186        the needed digits.  The denominator is 1.0 multiplied by the
01187        exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
01188        123e-6 gives 123 / 1000000.  */
01189 
01190     int expbit;
01191     int neg_exp;
01192     int more_bits;
01193     mp_limb_t cy;
01194     mp_limb_t *psrc = den;
01195     mp_limb_t *pdest = num;
01196     const struct mp_power *ttab = &_fpioconst_pow10[0];
01197 
01198     assert (dig_no > int_no && exponent <= 0);
01199 
01200 
01201     /* For the fractional part we need not process too many digits.  One
01202        decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
01203                         ceil(BITS / 3) =: N
01204        digits we should have enough bits for the result.  The remaining
01205        decimal digits give us the information that more bits are following.
01206        This can be used while rounding.  (Two added as a safety margin.)  */
01207     if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 2)
01208       {
01209         dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 2;
01210         more_bits = 1;
01211       }
01212     else
01213       more_bits = 0;
01214 
01215     neg_exp = dig_no - int_no - exponent;
01216 
01217     /* Construct the denominator.  */
01218     densize = 0;
01219     expbit = 1;
01220     do
01221       {
01222        if ((neg_exp & expbit) != 0)
01223          {
01224            mp_limb_t cy;
01225            neg_exp ^= expbit;
01226 
01227            if (densize == 0)
01228              {
01229               densize = ttab->arraysize - _FPIO_CONST_OFFSET;
01230               memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET],
01231                      densize * sizeof (mp_limb_t));
01232              }
01233            else
01234              {
01235               cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
01236                                          + _FPIO_CONST_OFFSET],
01237                             ttab->arraysize - _FPIO_CONST_OFFSET,
01238                             psrc, densize);
01239               densize += ttab->arraysize - _FPIO_CONST_OFFSET;
01240               if (cy == 0)
01241                 --densize;
01242               (void) SWAP (psrc, pdest);
01243              }
01244          }
01245        expbit <<= 1;
01246        ++ttab;
01247       }
01248     while (neg_exp != 0);
01249 
01250     if (psrc == num)
01251       memcpy (den, num, densize * sizeof (mp_limb_t));
01252 
01253     /* Read the fractional digits from the string.  */
01254     (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent
01255 #ifndef USE_WIDE_CHAR
01256                      , decimal, decimal_len, thousands
01257 #endif
01258                      );
01259 
01260     /* We now have to shift both numbers so that the highest bit in the
01261        denominator is set.  In the same process we copy the numerator to
01262        a high place in the array so that the division constructs the wanted
01263        digits.  This is done by a "quasi fix point" number representation.
01264 
01265        num:   ddddddddddd . 0000000000000000000000
01266               |--- m ---|
01267        den:                            ddddddddddd      n >= m
01268                                        |--- n ---|
01269      */
01270 
01271     count_leading_zeros (cnt, den[densize - 1]);
01272 
01273     if (cnt > 0)
01274       {
01275        /* Don't call `mpn_shift' with a count of zero since the specification
01276           does not allow this.  */
01277        (void) __mpn_lshift (den, den, densize, cnt);
01278        cy = __mpn_lshift (num, num, numsize, cnt);
01279        if (cy != 0)
01280          num[numsize++] = cy;
01281       }
01282 
01283     /* Now we are ready for the division.  But it is not necessary to
01284        do a full multi-precision division because we only need a small
01285        number of bits for the result.  So we do not use __mpn_divmod
01286        here but instead do the division here by hand and stop whenever
01287        the needed number of bits is reached.  The code itself comes
01288        from the GNU MP Library by Torbj\"orn Granlund.  */
01289 
01290     exponent = bits;
01291 
01292     switch (densize)
01293       {
01294       case 1:
01295        {
01296          mp_limb_t d, n, quot;
01297          int used = 0;
01298 
01299          n = num[0];
01300          d = den[0];
01301          assert (numsize == 1 && n < d);
01302 
01303          do
01304            {
01305              udiv_qrnnd (quot, n, n, 0, d);
01306 
01307 #define got_limb                                                     \
01308              if (bits == 0)                                          \
01309               {                                                      \
01310                 register int cnt;                                    \
01311                 if (quot == 0)                                       \
01312                   cnt = BITS_PER_MP_LIMB;                            \
01313                 else                                                 \
01314                   count_leading_zeros (cnt, quot);                          \
01315                 exponent -= cnt;                                     \
01316                 if (BITS_PER_MP_LIMB - cnt > MANT_DIG)               \
01317                   {                                                  \
01318                     used = MANT_DIG + cnt;                                  \
01319                     retval[0] = quot >> (BITS_PER_MP_LIMB - used);          \
01320                     bits = MANT_DIG + 1;                             \
01321                   }                                                  \
01322                 else                                                 \
01323                   {                                                  \
01324                     /* Note that we only clear the second element.  */      \
01325                     /* The conditional is determined at compile time.  */   \
01326                     if (RETURN_LIMB_SIZE > 1)                               \
01327                      retval[1] = 0;                                         \
01328                     retval[0] = quot;                                       \
01329                     bits = -cnt;                                     \
01330                   }                                                  \
01331               }                                                      \
01332              else if (bits + BITS_PER_MP_LIMB <= MANT_DIG)                  \
01333               __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,   \
01334                             quot);                                   \
01335              else                                                    \
01336               {                                                      \
01337                 used = MANT_DIG - bits;                              \
01338                 if (used > 0)                                               \
01339                   __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);    \
01340               }                                                      \
01341              bits += BITS_PER_MP_LIMB
01342 
01343              got_limb;
01344            }
01345          while (bits <= MANT_DIG);
01346 
01347          return round_and_return (retval, exponent - 1, negative,
01348                                quot, BITS_PER_MP_LIMB - 1 - used,
01349                                more_bits || n != 0);
01350        }
01351       case 2:
01352        {
01353          mp_limb_t d0, d1, n0, n1;
01354          mp_limb_t quot = 0;
01355          int used = 0;
01356 
01357          d0 = den[0];
01358          d1 = den[1];
01359 
01360          if (numsize < densize)
01361            {
01362              if (num[0] >= d1)
01363               {
01364                 /* The numerator of the number occupies fewer bits than
01365                    the denominator but the one limb is bigger than the
01366                    high limb of the numerator.  */
01367                 n1 = 0;
01368                 n0 = num[0];
01369               }
01370              else
01371               {
01372                 if (bits <= 0)
01373                   exponent -= BITS_PER_MP_LIMB;
01374                 else
01375                   {
01376                     if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
01377                      __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
01378                                    BITS_PER_MP_LIMB, 0);
01379                     else
01380                      {
01381                        used = MANT_DIG - bits;
01382                        if (used > 0)
01383                          __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
01384                      }
01385                     bits += BITS_PER_MP_LIMB;
01386                   }
01387                 n1 = num[0];
01388                 n0 = 0;
01389               }
01390            }
01391          else
01392            {
01393              n1 = num[1];
01394              n0 = num[0];
01395            }
01396 
01397          while (bits <= MANT_DIG)
01398            {
01399              mp_limb_t r;
01400 
01401              if (n1 == d1)
01402               {
01403                 /* QUOT should be either 111..111 or 111..110.  We need
01404                    special treatment of this rare case as normal division
01405                    would give overflow.  */
01406                 quot = ~(mp_limb_t) 0;
01407 
01408                 r = n0 + d1;
01409                 if (r < d1) /* Carry in the addition?  */
01410                   {
01411                     add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
01412                     goto have_quot;
01413                   }
01414                 n1 = d0 - (d0 != 0);
01415                 n0 = -d0;
01416               }
01417              else
01418               {
01419                 udiv_qrnnd (quot, r, n1, n0, d1);
01420                 umul_ppmm (n1, n0, d0, quot);
01421               }
01422 
01423            q_test:
01424              if (n1 > r || (n1 == r && n0 > 0))
01425               {
01426                 /* The estimated QUOT was too large.  */
01427                 --quot;
01428 
01429                 sub_ddmmss (n1, n0, n1, n0, 0, d0);
01430                 r += d1;
01431                 if (r >= d1)       /* If not carry, test QUOT again.  */
01432                   goto q_test;
01433               }
01434              sub_ddmmss (n1, n0, r, 0, n1, n0);
01435 
01436            have_quot:
01437              got_limb;
01438            }
01439 
01440          return round_and_return (retval, exponent - 1, negative,
01441                                quot, BITS_PER_MP_LIMB - 1 - used,
01442                                more_bits || n1 != 0 || n0 != 0);
01443        }
01444       default:
01445        {
01446          int i;
01447          mp_limb_t cy, dX, d1, n0, n1;
01448          mp_limb_t quot = 0;
01449          int used = 0;
01450 
01451          dX = den[densize - 1];
01452          d1 = den[densize - 2];
01453 
01454          /* The division does not work if the upper limb of the two-limb
01455             numerator is greater than the denominator.  */
01456          if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
01457            num[numsize++] = 0;
01458 
01459          if (numsize < densize)
01460            {
01461              mp_size_t empty = densize - numsize;
01462              register int i;
01463 
01464              if (bits <= 0)
01465               exponent -= empty * BITS_PER_MP_LIMB;
01466              else
01467               {
01468                 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
01469                   {
01470                     /* We make a difference here because the compiler
01471                       cannot optimize the `else' case that good and
01472                       this reflects all currently used FLOAT types
01473                       and GMP implementations.  */
01474 #if RETURN_LIMB_SIZE <= 2
01475                     assert (empty == 1);
01476                     __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
01477                                   BITS_PER_MP_LIMB, 0);
01478 #else
01479                     for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i)
01480                      retval[i] = retval[i - empty];
01481                     while (i >= 0)
01482                      retval[i--] = 0;
01483 #endif
01484                   }
01485                 else
01486                   {
01487                     used = MANT_DIG - bits;
01488                     if (used >= BITS_PER_MP_LIMB)
01489                      {
01490                        register int i;
01491                        (void) __mpn_lshift (&retval[used
01492                                                  / BITS_PER_MP_LIMB],
01493                                           retval, RETURN_LIMB_SIZE,
01494                                           used % BITS_PER_MP_LIMB);
01495                        for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i)
01496                          retval[i] = 0;
01497                      }
01498                     else if (used > 0)
01499                      __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
01500                   }
01501                 bits += empty * BITS_PER_MP_LIMB;
01502               }
01503              for (i = numsize; i > 0; --i)
01504               num[i + empty] = num[i - 1];
01505              MPN_ZERO (num, empty + 1);
01506            }
01507          else
01508            {
01509              int i;
01510              assert (numsize == densize);
01511              for (i = numsize; i > 0; --i)
01512               num[i] = num[i - 1];
01513            }
01514 
01515          den[densize] = 0;
01516          n0 = num[densize];
01517 
01518          while (bits <= MANT_DIG)
01519            {
01520              if (n0 == dX)
01521               /* This might over-estimate QUOT, but it's probably not
01522                  worth the extra code here to find out.  */
01523               quot = ~(mp_limb_t) 0;
01524              else
01525               {
01526                 mp_limb_t r;
01527 
01528                 udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
01529                 umul_ppmm (n1, n0, d1, quot);
01530 
01531                 while (n1 > r || (n1 == r && n0 > num[densize - 2]))
01532                   {
01533                     --quot;
01534                     r += dX;
01535                     if (r < dX) /* I.e. "carry in previous addition?" */
01536                      break;
01537                     n1 -= n0 < d1;
01538                     n0 -= d1;
01539                   }
01540               }
01541 
01542              /* Possible optimization: We already have (q * n0) and (1 * n1)
01543                after the calculation of QUOT.  Taking advantage of this, we
01544                could make this loop make two iterations less.  */
01545 
01546              cy = __mpn_submul_1 (num, den, densize + 1, quot);
01547 
01548              if (num[densize] != cy)
01549               {
01550                 cy = __mpn_add_n (num, num, den, densize);
01551                 assert (cy != 0);
01552                 --quot;
01553               }
01554              n0 = num[densize] = num[densize - 1];
01555              for (i = densize - 1; i > 0; --i)
01556               num[i] = num[i - 1];
01557 
01558              got_limb;
01559            }
01560 
01561          for (i = densize; num[i] == 0 && i >= 0; --i)
01562            ;
01563          return round_and_return (retval, exponent - 1, negative,
01564                                quot, BITS_PER_MP_LIMB - 1 - used,
01565                                more_bits || i >= 0);
01566        }
01567       }
01568   }
01569 
01570   /* NOTREACHED */
01571 }
01572 #if defined _LIBC && !defined USE_WIDE_CHAR
01573 libc_hidden_def (____STRTOF_INTERNAL)
01574 #endif
01575 
01576 /* External user entry point.  */
01577 
01578 FLOAT
01579 #ifdef weak_function
01580 weak_function
01581 #endif
01582 __STRTOF (nptr, endptr, loc)
01583      const STRING_TYPE *nptr;
01584      STRING_TYPE **endptr;
01585      __locale_t loc;
01586 {
01587   return ____STRTOF_INTERNAL (nptr, endptr, 0, loc);
01588 }
01589 #if defined _LIBC
01590 libc_hidden_def (__STRTOF)
01591 libc_hidden_ver (__STRTOF, STRTOF)
01592 #endif
01593 weak_alias (__STRTOF, STRTOF)
01594 
01595 #ifdef LONG_DOUBLE_COMPAT
01596 # if LONG_DOUBLE_COMPAT(libc, GLIBC_2_1)
01597 #  ifdef USE_WIDE_CHAR
01598 compat_symbol (libc, __wcstod_l, __wcstold_l, GLIBC_2_1);
01599 #  else
01600 compat_symbol (libc, __strtod_l, __strtold_l, GLIBC_2_1);
01601 #  endif
01602 # endif
01603 # if LONG_DOUBLE_COMPAT(libc, GLIBC_2_3)
01604 #  ifdef USE_WIDE_CHAR
01605 compat_symbol (libc, wcstod_l, wcstold_l, GLIBC_2_3);
01606 #  else
01607 compat_symbol (libc, strtod_l, strtold_l, GLIBC_2_3);
01608 #  endif
01609 # endif
01610 #endif