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