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