OSDN Git Service

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