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