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