3818c81ecb0c821f0899d25583d1c7fafd038887
[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     {
485       if (isdigit (cp[1]))
486         {
487           c = *++cp;
488           do
489             {
490               if (c != '0' && lead_zero == -1)
491                 lead_zero = dig_no - int_no;
492               ++dig_no;
493               c = *++cp;
494             }
495           while (isdigit (c));
496         }
497     }
498
499   /* Remember start of exponent (if any).  */
500   expp = cp;
501
502   /* Read exponent.  */
503   if (tolower (c) == 'e')
504     {
505       int exp_negative = 0;
506
507       c = *++cp;
508       if (c == '-')
509         {
510           exp_negative = 1;
511           c = *++cp;
512         }
513       else if (c == '+')
514         c = *++cp;
515
516       if (isdigit (c))
517         {
518           int exp_limit;
519
520           /* Get the exponent limit. */
521           exp_limit = exp_negative ?
522                 -MIN_10_EXP + MANT_DIG - int_no :
523                 MAX_10_EXP - int_no + lead_zero;
524
525           do
526             {
527               exponent *= 10;
528
529               if (exponent > exp_limit)
530                 /* The exponent is too large/small to represent a valid
531                    number.  */
532                 {
533                   FLOAT retval;
534
535                   /* Overflow or underflow.  */
536                   errno = ERANGE;
537                   retval = (exp_negative ? 0.0 :
538                             negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
539
540                   /* Accept all following digits as part of the exponent.  */
541                   do
542                     ++cp;
543                   while (isdigit (*cp));
544
545                   RETURN (retval, cp);
546                   /* NOTREACHED */
547                 }
548
549               exponent += c - '0';
550               c = *++cp;
551             }
552           while (isdigit (c));
553
554           if (exp_negative)
555             exponent = -exponent;
556         }
557       else
558         cp = expp;
559     }
560
561   /* We don't want to have to work with trailing zeroes after the radix.  */
562   if (dig_no > int_no)
563     {
564       while (expp[-1] == '0')
565         {
566           --expp;
567           --dig_no;
568         }
569       assert (dig_no >= int_no);
570     }
571
572  number_parsed:
573
574   /* The whole string is parsed.  Store the address of the next character.  */
575   if (endptr)
576     *endptr = (char *) cp;
577
578   if (dig_no == 0)
579     return 0.0;
580
581   if (lead_zero)
582     {
583       /* Find the decimal point */
584       while (*startp != decimal) startp++;
585       startp += lead_zero + 1;
586       exponent -= lead_zero;
587       dig_no -= lead_zero;
588     }
589
590   /* Now we have the number of digits in total and the integer digits as well
591      as the exponent and its sign.  We can decide whether the read digits are
592      really integer digits or belong to the fractional part; i.e. we normalize
593      123e-2 to 1.23.  */
594   {
595     register int incr = exponent < 0 ? MAX (-int_no, exponent)
596                                      : MIN (dig_no - int_no, exponent);
597     int_no += incr;
598     exponent -= incr;
599   }
600
601   if (int_no + exponent > MAX_10_EXP + 1)
602     {
603       errno = ERANGE;
604       return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
605     }
606
607   if (exponent < MIN_10_EXP - (DIG + 1))
608     {
609       errno = ERANGE;
610       return 0.0;
611     }
612
613   if (int_no > 0)
614     {
615       /* Read the integer part as a multi-precision number to NUM.  */
616       startp = str_to_mpn (startp, int_no, num, &numsize, &exponent);
617
618       if (exponent > 0)
619         {
620           /* We now multiply the gained number by the given power of ten.  */
621           mp_limb *psrc = num;
622           mp_limb *pdest = den;
623           int expbit = 1;
624           const struct mp_power *ttab = &_fpioconst_pow10[0];
625
626           do
627             {
628               if ((exponent & expbit) != 0)
629                 {
630                   mp_limb cy;
631                   exponent ^= expbit;
632
633                   /* FIXME: not the whole multiplication has to be done.
634                      If we have the needed number of bits we only need the
635                      information whether more non-zero bits follow.  */
636                   if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
637                     cy = __mpn_mul (pdest, psrc, numsize,
638                                     &ttab->array[_FPIO_CONST_OFFSET],
639                                     ttab->arraysize - _FPIO_CONST_OFFSET);
640                   else
641                     cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
642                                     ttab->arraysize - _FPIO_CONST_OFFSET,
643                                     psrc, numsize);
644                   numsize += ttab->arraysize - _FPIO_CONST_OFFSET;
645                   if (cy == 0)
646                     --numsize;
647                   SWAP (psrc, pdest);
648                 }
649               expbit <<= 1;
650               ++ttab;
651             }
652           while (exponent != 0);
653
654           if (psrc == den)
655             memcpy (num, den, numsize * sizeof (mp_limb));
656         }
657
658       /* Determine how many bits of the result we already have.  */
659       count_leading_zeros (bits, num[numsize - 1]);
660       bits = numsize * BITS_PER_MP_LIMB - bits;
661
662       /* Now we know the exponent of the number in base two.
663          Check it against the maximum possible exponent.  */
664       if (bits > MAX_EXP)
665         {
666           errno = ERANGE;
667           return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
668         }
669
670       /* We have already the first BITS bits of the result.  Together with
671          the information whether more non-zero bits follow this is enough
672          to determine the result.  */
673       if (bits > MANT_DIG)
674         {
675           int i;
676           const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
677           const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
678           const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
679                                                      : least_idx;
680           const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
681                                                      : least_bit - 1;
682
683           if (least_bit == 0)
684             memcpy (retval, &num[least_idx],
685                     RETURN_LIMB_SIZE * sizeof (mp_limb));
686           else
687             {
688               for (i = least_idx; i < numsize - 1; ++i)
689                 retval[i - least_idx] = (num[i] >> least_bit)
690                                         | (num[i + 1]
691                                            << (BITS_PER_MP_LIMB - least_bit));
692               if (i - least_idx < RETURN_LIMB_SIZE)
693                 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
694             }
695
696           /* Check whether any limb beside the ones in RETVAL are non-zero.  */
697           for (i = 0; num[i] == 0; ++i)
698             ;
699
700           return round_and_return (retval, bits - 1, negative,
701                                    num[round_idx], round_bit,
702                                    int_no < dig_no || i < round_idx);
703           /* NOTREACHED */
704         }
705       else if (dig_no == int_no)
706         {
707           const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
708           const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
709
710           if (target_bit == is_bit)
711             {
712               memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
713                       numsize * sizeof (mp_limb));
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 if (target_bit > is_bit)
719             {
720               (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
721                                    num, numsize, target_bit - is_bit);
722               /* FIXME: the following loop can be avoided if we assume a
723                  maximal MANT_DIG value.  */
724               MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
725             }
726           else
727             {
728               mp_limb cy;
729               assert (numsize < RETURN_LIMB_SIZE);
730
731               cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
732                                  num, numsize, is_bit - target_bit);
733               retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
734               /* FIXME: the following loop can be avoided if we assume a
735                  maximal MANT_DIG value.  */
736               MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
737             }
738
739           return round_and_return (retval, bits - 1, negative, 0, 0, 0);
740           /* NOTREACHED */
741         }
742
743       /* Store the bits we already have.  */
744       memcpy (retval, num, numsize * sizeof (mp_limb));
745 #if RETURN_LIMB_SIZE > 1
746       if (numsize < RETURN_LIMB_SIZE)
747         retval[numsize] = 0;
748 #endif
749     }
750
751   /* We have to compute at least some of the fractional digits.  */
752   {
753     /* We construct a fraction and the result of the division gives us
754        the needed digits.  The denominator is 1.0 multiplied by the
755        exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
756        123e-6 gives 123 / 1000000.  */
757
758     int expbit;
759     int cnt;
760     int neg_exp;
761     int more_bits;
762     mp_limb cy;
763     mp_limb *psrc = den;
764     mp_limb *pdest = num;
765     const struct mp_power *ttab = &_fpioconst_pow10[0];
766
767     assert (dig_no > int_no && exponent <= 0);
768
769
770     /* For the fractional part we need not process too much digits.  One
771        decimal digits gives us log_2(10) ~ 3.32 bits.  If we now compute
772                         ceil(BITS / 3) =: N
773        digits we should have enough bits for the result.  The remaining
774        decimal digits give us the information that more bits are following.
775        This can be used while rounding.  (One added as a safety margin.)  */
776     if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
777       {
778         dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
779         more_bits = 1;
780       }
781     else
782       more_bits = 0;
783
784     neg_exp = dig_no - int_no - exponent;
785
786     /* Construct the denominator.  */
787     densize = 0;
788     expbit = 1;
789     do
790       {
791         if ((neg_exp & expbit) != 0)
792           {
793             mp_limb cy;
794             neg_exp ^= expbit;
795
796             if (densize == 0)
797               {
798                 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
799                 memcpy (psrc, &ttab->array[_FPIO_CONST_OFFSET],
800                         densize * sizeof (mp_limb));
801               }
802             else
803               {
804                 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
805                                 ttab->arraysize - _FPIO_CONST_OFFSET,
806                                 psrc, densize);
807                 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
808                 if (cy == 0)
809                   --densize;
810                 SWAP (psrc, pdest);
811               }
812           }
813         expbit <<= 1;
814         ++ttab;
815       }
816     while (neg_exp != 0);
817
818     if (psrc == num)
819       memcpy (den, num, densize * sizeof (mp_limb));
820
821     /* Read the fractional digits from the string.  */
822     (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent);
823
824
825     /* We now have to shift both numbers so that the highest bit in the
826        denominator is set.  In the same process we copy the numerator to
827        a high place in the array so that the division constructs the wanted
828        digits.  This is done by a "quasi fix point" number representation.
829
830        num:   ddddddddddd . 0000000000000000000000
831               |--- m ---|
832        den:                            ddddddddddd      n >= m
833                                        |--- n ---|
834      */
835
836     count_leading_zeros (cnt, den[densize - 1]);
837
838     (void) __mpn_lshift (den, den, densize, cnt);
839     cy = __mpn_lshift (num, num, numsize, cnt);
840     if (cy != 0)
841       num[numsize++] = cy;
842
843     /* Now we are ready for the division.  But it is not necessary to
844        do a full multi-precision division because we only need a small
845        number of bits for the result.  So we do not use __mpn_divmod
846        here but instead do the division here by hand and stop whenever
847        the needed number of bits is reached.  The code itself comes
848        from the GNU MP Library by Torbj\"orn Granlund.  */
849
850     exponent = bits;
851
852     switch (densize)
853       {
854       case 1:
855         {
856           mp_limb d, n, quot;
857           int used = 0;
858
859           n = num[0];
860           d = den[0];
861           assert (numsize == 1 && n < d);
862
863           do
864             {
865               udiv_qrnnd (quot, n, n, 0, d);
866
867 #define got_limb                                                              \
868               if (bits == 0)                                                  \
869                 {                                                             \
870                   register int cnt;                                           \
871                   if (quot == 0)                                              \
872                     cnt = BITS_PER_MP_LIMB;                                   \
873                   else                                                        \
874                     count_leading_zeros (cnt, quot);                          \
875                   exponent -= cnt;                                            \
876                   if (BITS_PER_MP_LIMB - cnt > MANT_DIG)                      \
877                     {                                                         \
878                       used = MANT_DIG + cnt;                                  \
879                       retval[0] = quot >> (BITS_PER_MP_LIMB - used);          \
880                       bits = MANT_DIG + 1;                                    \
881                     }                                                         \
882                   else                                                        \
883                     {                                                         \
884                       /* Note that we only clear the second element.  */      \
885                       /* The conditional is determined at compile time.  */   \
886                       if (RETURN_LIMB_SIZE > 1)                               \
887                         retval[1] = 0;                                        \
888                       retval[0] = quot;                                       \
889                       bits = -cnt;                                            \
890                     }                                                         \
891                 }                                                             \
892               else if (bits + BITS_PER_MP_LIMB <= MANT_DIG)                   \
893                 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB,   \
894                                 quot);                                        \
895               else                                                            \
896                 {                                                             \
897                   used = MANT_DIG - bits;                                     \
898                   if (used > 0)                                               \
899                     __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot);    \
900                 }                                                             \
901               bits += BITS_PER_MP_LIMB
902
903               got_limb;
904             }
905           while (bits <= MANT_DIG);
906
907           return round_and_return (retval, exponent - 1, negative,
908                                    quot, BITS_PER_MP_LIMB - 1 - used,
909                                    more_bits || n != 0);
910         }
911       case 2:
912         {
913           mp_limb d0, d1, n0, n1;
914           mp_limb quot = 0;
915           int used = 0;
916
917           d0 = den[0];
918           d1 = den[1];
919
920           if (numsize < densize)
921             {
922               if (num[0] >= d1)
923                 {
924                   /* The numerator of the number occupies fewer bits than
925                      the denominator but the one limb is bigger than the
926                      high limb of the numerator.  */
927                   n1 = 0;
928                   n0 = num[0];
929                 }
930               else
931                 {
932                   if (bits <= 0)
933                     exponent -= BITS_PER_MP_LIMB;
934                   else
935                     {
936                       if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
937                         __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
938                                         BITS_PER_MP_LIMB, 0);
939                       else
940                         {
941                           used = MANT_DIG - bits;
942                           if (used > 0)
943                             __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
944                         }
945                       bits += BITS_PER_MP_LIMB;
946                     }
947                   n1 = num[0];
948                   n0 = 0;
949                 }
950             }
951           else
952             {
953               n1 = num[1];
954               n0 = num[0];
955             }
956
957           while (bits <= MANT_DIG)
958             {
959               mp_limb r;
960
961               if (n1 == d1)
962                 {
963                   /* QUOT should be either 111..111 or 111..110.  We need
964                      special treatment of this rare case as normal division
965                      would give overflow.  */
966                   quot = ~(mp_limb) 0;
967
968                   r = n0 + d1;
969                   if (r < d1)   /* Carry in the addition?  */
970                     {
971                       add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
972                       goto have_quot;
973                     }
974                   n1 = d0 - (d0 != 0);
975                   n0 = -d0;
976                 }
977               else
978                 {
979                   udiv_qrnnd (quot, r, n1, n0, d1);
980                   umul_ppmm (n1, n0, d0, quot);
981                 }
982
983             q_test:
984               if (n1 > r || (n1 == r && n0 > 0))
985                 {
986                   /* The estimated QUOT was too large.  */
987                   --quot;
988
989                   sub_ddmmss (n1, n0, n1, n0, 0, d0);
990                   r += d1;
991                   if (r >= d1)  /* If not carry, test QUOT again.  */
992                     goto q_test;
993                 }
994               sub_ddmmss (n1, n0, r, 0, n1, n0);
995
996             have_quot:
997               got_limb;
998             }
999
1000           return round_and_return (retval, exponent - 1, negative,
1001                                    quot, BITS_PER_MP_LIMB - 1 - used,
1002                                    more_bits || n1 != 0 || n0 != 0);
1003         }
1004       default:
1005         {
1006           int i;
1007           mp_limb cy, dX, d1, n0, n1;
1008           mp_limb quot = 0;
1009           int used = 0;
1010
1011           dX = den[densize - 1];
1012           d1 = den[densize - 2];
1013
1014           /* The division does not work if the upper limb of the two-limb
1015              numerator is greater than the denominator.  */
1016           if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1017             num[numsize++] = 0;
1018
1019           if (numsize < densize)
1020             {
1021               mp_size_t empty = densize - numsize;
1022
1023               if (bits <= 0)
1024                 {
1025                   register int i;
1026                   for (i = numsize; i > 0; --i)
1027                     num[i + empty] = num[i - 1];
1028                   MPN_ZERO (num, empty + 1);
1029                   exponent -= empty * BITS_PER_MP_LIMB;
1030                 }
1031               else
1032                 {
1033                   if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1034                     {
1035                       /* We make a difference here because the compiler
1036                          cannot optimize the `else' case that good and
1037                          this reflects all currently used FLOAT types
1038                          and GMP implementations.  */
1039                       register int i;
1040 #if RETURN_LIMB_SIZE <= 2
1041                       assert (empty == 1);
1042                       __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1043                                       BITS_PER_MP_LIMB, 0);
1044 #else
1045                       for (i = RETURN_LIMB_SIZE; i > empty; --i)
1046                         retval[i] = retval[i - empty];
1047 #endif
1048                       retval[1] = 0;
1049                       for (i = numsize; i > 0; --i)
1050                         num[i + empty] = num[i - 1];
1051                       MPN_ZERO (num, empty + 1);
1052                     }
1053                   else
1054                     {
1055                       used = MANT_DIG - bits;
1056                       if (used >= BITS_PER_MP_LIMB)
1057                         {
1058                           register int i;
1059                           (void) __mpn_lshift (&retval[used
1060                                                        / BITS_PER_MP_LIMB],
1061                                                retval, RETURN_LIMB_SIZE,
1062                                                used % BITS_PER_MP_LIMB);
1063                           for (i = used / BITS_PER_MP_LIMB; i >= 0; --i)
1064                             retval[i] = 0;
1065                         }
1066                       else if (used > 0)
1067                         __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1068                     }
1069                   bits += empty * BITS_PER_MP_LIMB;
1070                 }
1071             }
1072           else
1073             {
1074               int i;
1075               assert (numsize == densize);
1076               for (i = numsize; i > 0; --i)
1077                 num[i] = num[i - 1];
1078             }
1079
1080           den[densize] = 0;
1081           n0 = num[densize];
1082
1083           while (bits <= MANT_DIG)
1084             {
1085               if (n0 == dX)
1086                 /* This might over-estimate QUOT, but it's probably not
1087                    worth the extra code here to find out.  */
1088                 quot = ~(mp_limb) 0;
1089               else
1090                 {
1091                   mp_limb r;
1092
1093                   udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1094                   umul_ppmm (n1, n0, d1, quot);
1095
1096                   while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1097                     {
1098                       --quot;
1099                       r += dX;
1100                       if (r < dX) /* I.e. "carry in previous addition?" */
1101                         break;
1102                       n1 -= n0 < d1;
1103                       n0 -= d1;
1104                     }
1105                 }
1106
1107               /* Possible optimization: We already have (q * n0) and (1 * n1)
1108                  after the calculation of QUOT.  Taking advantage of this, we
1109                  could make this loop make two iterations less.  */
1110
1111               cy = __mpn_submul_1 (num, den, densize + 1, quot);
1112
1113               if (num[densize] != cy)
1114                 {
1115                   cy = __mpn_add_n (num, num, den, densize);
1116                   assert (cy != 0);
1117                   --quot;
1118                 }
1119               n0 = num[densize] = num[densize - 1];
1120               for (i = densize - 1; i > 0; --i)
1121                 num[i] = num[i - 1];
1122
1123               got_limb;
1124             }
1125
1126           for (i = densize; num[i] == 0 && i >= 0; --i)
1127             ;
1128           return round_and_return (retval, exponent - 1, negative,
1129                                    quot, BITS_PER_MP_LIMB - 1 - used,
1130                                    more_bits || i >= 0);
1131         }
1132       }
1133   }
1134
1135   /* NOTREACHED */
1136 }
1137 \f
1138 /* External user entry point.  */
1139
1140 FLOAT
1141 STRTOF (nptr, endptr)
1142      const char *nptr;
1143      char **endptr;
1144 {
1145   return INTERNAL (STRTOF) (nptr, endptr, 0);
1146 }
1147
1148 #define weak_this(x) weak_symbol(x)
1149 weak_this (STRTOF)