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