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