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