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