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