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