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