Update.
[glibc.git] / stdlib / strtod.c
blob061cedc98a0a80db5e02d8f09cd1b23a4edac02e
1 /* Read decimal floating point numbers.
2 This file is part of the GNU C Library.
3 Copyright (C) 1995, 1996, 1997, 1998 Free Software Foundation, Inc.
4 Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Library General Public License as
8 published by the Free Software Foundation; either version 2 of the
9 License, or (at your option) any later version.
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Library General Public License for more details.
16 You should have received a copy of the GNU Library General Public
17 License along with the GNU C Library; see the file COPYING.LIB. If not,
18 write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19 Boston, MA 02111-1307, USA. */
21 /* Configuration part. These macros are defined by `strtold.c',
22 `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
23 `long double' and `float' versions of the reader. */
24 #ifndef FLOAT
25 # define FLOAT double
26 # define FLT DBL
27 # ifdef USE_WIDE_CHAR
28 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
29 # define STRTOF __wcstod_l
30 # else
31 # define STRTOF wcstod
32 # endif
33 # else
34 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
35 # define STRTOF __strtod_l
36 # else
37 # define STRTOF strtod
38 # endif
39 # endif
40 # define MPN2FLOAT __mpn_construct_double
41 # define FLOAT_HUGE_VAL HUGE_VAL
42 # define SET_MANTISSA(flt, mant) \
43 do { union ieee754_double u; \
44 u.d = (flt); \
45 if ((mant & 0xfffffffffffffULL) == 0) \
46 mant = 0x8000000000000ULL; \
47 u.ieee.mantissa0 = ((mant) >> 32) & 0xfffff; \
48 u.ieee.mantissa1 = (mant) & 0xffffffff; \
49 (flt) = u.d; \
50 } while (0)
51 #endif
52 /* End of configuration part. */
54 #include <ctype.h>
55 #include <errno.h>
56 #include <float.h>
57 #include <ieee754.h>
58 #include "../locale/localeinfo.h"
59 #include <math.h>
60 #include <stdlib.h>
61 #include <string.h>
63 /* The gmp headers need some configuration frobs. */
64 #define HAVE_ALLOCA 1
66 #include <gmp.h>
67 #include <gmp-impl.h>
68 #include <gmp-mparam.h>
69 #include <longlong.h>
70 #include "fpioconst.h"
72 #define NDEBUG 1
73 #include <assert.h>
76 /* We use this code also for the extended locale handling where the
77 function gets as an additional argument the locale which has to be
78 used. To access the values we have to redefine the _NL_CURRENT
79 macro. */
80 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
81 # undef _NL_CURRENT
82 # define _NL_CURRENT(category, item) \
83 (current->values[_NL_ITEM_INDEX (item)].string)
84 # define LOCALE_PARAM , loc
85 # define LOCALE_PARAM_DECL __locale_t loc;
86 #else
87 # define LOCALE_PARAM
88 # define LOCALE_PARAM_DECL
89 #endif
92 #ifdef USE_WIDE_CHAR
93 # include <wctype.h>
94 # include <wchar.h>
95 # define STRING_TYPE wchar_t
96 # define CHAR_TYPE wint_t
97 # define L_(Ch) L##Ch
98 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
99 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
100 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
101 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
102 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
103 # define STRNCASECMP(S1, S2, N) __wcsncasecmp_l ((S1), (S2), (N), loc)
104 # define STRTOULL(S, E, B) __wcstoull_l ((S), (E), (B), loc)
105 # else
106 # define ISSPACE(Ch) iswspace (Ch)
107 # define ISDIGIT(Ch) iswdigit (Ch)
108 # define ISXDIGIT(Ch) iswxdigit (Ch)
109 # define TOLOWER(Ch) towlower (Ch)
110 # define STRNCASECMP(S1, S2, N) __wcsncasecmp ((S1), (S2), (N))
111 # define STRTOULL(S, E, B) wcstoull ((S), (E), (B))
112 # endif
113 #else
114 # define STRING_TYPE char
115 # define CHAR_TYPE char
116 # define L_(Ch) Ch
117 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
118 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
119 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
120 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
121 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
122 # define STRNCASECMP(S1, S2, N) __strncasecmp_l ((S1), (S2), (N), loc)
123 # define STRTOULL(S, E, B) __strtoull_l ((S), (E), (B), loc)
124 # else
125 # define ISSPACE(Ch) isspace (Ch)
126 # define ISDIGIT(Ch) isdigit (Ch)
127 # define ISXDIGIT(Ch) isxdigit (Ch)
128 # define TOLOWER(Ch) tolower (Ch)
129 # define STRNCASECMP(S1, S2, N) __strncasecmp ((S1), (S2), (N))
130 # define STRTOULL(S, E, B) strtoull ((S), (E), (B))
131 # endif
132 #endif
135 /* Constants we need from float.h; select the set for the FLOAT precision. */
136 #define MANT_DIG PASTE(FLT,_MANT_DIG)
137 #define DIG PASTE(FLT,_DIG)
138 #define MAX_EXP PASTE(FLT,_MAX_EXP)
139 #define MIN_EXP PASTE(FLT,_MIN_EXP)
140 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
141 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
143 /* Extra macros required to get FLT expanded before the pasting. */
144 #define PASTE(a,b) PASTE1(a,b)
145 #define PASTE1(a,b) a##b
147 /* Function to construct a floating point number from an MP integer
148 containing the fraction bits, a base 2 exponent, and a sign flag. */
149 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
151 /* Definitions according to limb size used. */
152 #if BITS_PER_MP_LIMB == 32
153 # define MAX_DIG_PER_LIMB 9
154 # define MAX_FAC_PER_LIMB 1000000000UL
155 #elif BITS_PER_MP_LIMB == 64
156 # define MAX_DIG_PER_LIMB 19
157 # define MAX_FAC_PER_LIMB 10000000000000000000UL
158 #else
159 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
160 #endif
163 /* Local data structure. */
164 static const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] =
165 { 0, 10, 100,
166 1000, 10000, 100000,
167 1000000, 10000000, 100000000,
168 1000000000
169 #if BITS_PER_MP_LIMB > 32
170 , 10000000000U, 100000000000U,
171 1000000000000U, 10000000000000U, 100000000000000U,
172 1000000000000000U, 10000000000000000U, 100000000000000000U,
173 1000000000000000000U, 10000000000000000000U
174 #endif
175 #if BITS_PER_MP_LIMB > 64
176 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
177 #endif
180 #ifndef howmany
181 #define howmany(x,y) (((x)+((y)-1))/(y))
182 #endif
183 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
185 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
186 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
187 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
189 #define RETURN(val,end) \
190 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
191 return val; } while (0)
193 /* Maximum size necessary for mpn integers to hold floating point numbers. */
194 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
195 + 2)
196 /* Declare an mpn integer variable that big. */
197 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
198 /* Copy an mpn integer value. */
199 #define MPN_ASSIGN(dst, src) \
200 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
203 /* Return a floating point number of the needed type according to the given
204 multi-precision number after possible rounding. */
205 static inline FLOAT
206 round_and_return (mp_limb_t *retval, int exponent, int negative,
207 mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
209 if (exponent < MIN_EXP - 1)
211 mp_size_t shift = MIN_EXP - 1 - exponent;
213 if (shift > MANT_DIG)
215 __set_errno (EDOM);
216 return 0.0;
219 more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
220 if (shift == MANT_DIG)
221 /* This is a special case to handle the very seldom case where
222 the mantissa will be empty after the shift. */
224 int i;
226 round_limb = retval[RETURN_LIMB_SIZE - 1];
227 round_bit = BITS_PER_MP_LIMB - 1;
228 for (i = 0; i < RETURN_LIMB_SIZE; ++i)
229 more_bits |= retval[i] != 0;
230 MPN_ZERO (retval, RETURN_LIMB_SIZE);
232 else if (shift >= BITS_PER_MP_LIMB)
234 int i;
236 round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
237 round_bit = (shift - 1) % BITS_PER_MP_LIMB;
238 for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
239 more_bits |= retval[i] != 0;
240 more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
241 != 0);
243 (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
244 RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
245 shift % BITS_PER_MP_LIMB);
246 MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
247 shift / BITS_PER_MP_LIMB);
249 else if (shift > 0)
251 round_limb = retval[0];
252 round_bit = shift - 1;
253 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
255 exponent = MIN_EXP - 2;
258 if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
259 && (more_bits || (retval[0] & 1) != 0
260 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
262 mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
264 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
265 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
266 (retval[RETURN_LIMB_SIZE - 1]
267 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
269 ++exponent;
270 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
271 retval[RETURN_LIMB_SIZE - 1]
272 |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
274 else if (exponent == MIN_EXP - 2
275 && (retval[RETURN_LIMB_SIZE - 1]
276 & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
277 != 0)
278 /* The number was denormalized but now normalized. */
279 exponent = MIN_EXP - 1;
282 if (exponent > MAX_EXP)
283 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
285 return MPN2FLOAT (retval, exponent, negative);
289 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
290 into N. Return the size of the number limbs in NSIZE at the first
291 character od the string that is not part of the integer as the function
292 value. If the EXPONENT is small enough to be taken as an additional
293 factor for the resulting number (see code) multiply by it. */
294 static inline const STRING_TYPE *
295 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
296 int *exponent)
298 /* Number of digits for actual limb. */
299 int cnt = 0;
300 mp_limb_t low = 0;
301 mp_limb_t start;
303 *nsize = 0;
304 assert (digcnt > 0);
307 if (cnt == MAX_DIG_PER_LIMB)
309 if (*nsize == 0)
310 n[0] = low;
311 else
313 mp_limb_t cy;
314 cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
315 cy += __mpn_add_1 (n, n, *nsize, low);
316 if (cy != 0)
317 n[*nsize] = cy;
319 ++(*nsize);
320 cnt = 0;
321 low = 0;
324 /* There might be thousands separators or radix characters in
325 the string. But these all can be ignored because we know the
326 format of the number is correct and we have an exact number
327 of characters to read. */
328 while (*str < L_('0') || *str > L_('9'))
329 ++str;
330 low = low * 10 + *str++ - L_('0');
331 ++cnt;
333 while (--digcnt > 0);
335 if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
337 low *= _tens_in_limb[*exponent];
338 start = _tens_in_limb[cnt + *exponent];
339 *exponent = 0;
341 else
342 start = _tens_in_limb[cnt];
344 if (*nsize == 0)
346 n[0] = low;
347 *nsize = 1;
349 else
351 mp_limb_t cy;
352 cy = __mpn_mul_1 (n, n, *nsize, start);
353 cy += __mpn_add_1 (n, n, *nsize, low);
354 if (cy != 0)
355 n[(*nsize)++] = cy;
358 return str;
362 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
363 with the COUNT most significant bits of LIMB.
365 Tege doesn't like this function so I have to write it here myself. :)
366 --drepper */
367 static inline void
368 __mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
369 mp_limb_t limb)
371 if (count == BITS_PER_MP_LIMB)
373 /* Optimize the case of shifting by exactly a word:
374 just copy words, with no actual bit-shifting. */
375 mp_size_t i;
376 for (i = size - 1; i > 0; --i)
377 ptr[i] = ptr[i - 1];
378 ptr[0] = limb;
380 else
382 (void) __mpn_lshift (ptr, ptr, size, count);
383 ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
388 #define INTERNAL(x) INTERNAL1(x)
389 #define INTERNAL1(x) __##x##_internal
391 /* This file defines a function to check for correct grouping. */
392 #include "grouping.h"
395 /* Return a floating point number with the value of the given string NPTR.
396 Set *ENDPTR to the character after the last used one. If the number is
397 smaller than the smallest representable number, set `errno' to ERANGE and
398 return 0.0. If the number is too big to be represented, set `errno' to
399 ERANGE and return HUGE_VAL with the appropriate sign. */
400 FLOAT
401 INTERNAL (STRTOF) (nptr, endptr, group LOCALE_PARAM)
402 const STRING_TYPE *nptr;
403 STRING_TYPE **endptr;
404 int group;
405 LOCALE_PARAM_DECL
407 int negative; /* The sign of the number. */
408 MPN_VAR (num); /* MP representation of the number. */
409 int exponent; /* Exponent of the number. */
411 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
412 int base = 10;
414 /* When we have to compute fractional digits we form a fraction with a
415 second multi-precision number (and we sometimes need a second for
416 temporary results). */
417 MPN_VAR (den);
419 /* Representation for the return value. */
420 mp_limb_t retval[RETURN_LIMB_SIZE];
421 /* Number of bits currently in result value. */
422 int bits;
424 /* Running pointer after the last character processed in the string. */
425 const STRING_TYPE *cp, *tp;
426 /* Start of significant part of the number. */
427 const STRING_TYPE *startp, *start_of_digits;
428 /* Points at the character following the integer and fractional digits. */
429 const STRING_TYPE *expp;
430 /* Total number of digit and number of digits in integer part. */
431 int dig_no, int_no, lead_zero;
432 /* Contains the last character read. */
433 CHAR_TYPE c;
435 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
436 there. So define it ourselves if it remains undefined. */
437 #ifndef _WINT_T
438 typedef unsigned int wint_t;
439 #endif
440 /* The radix character of the current locale. */
441 wchar_t decimal;
442 /* The thousands character of the current locale. */
443 wchar_t thousands;
444 /* The numeric grouping specification of the current locale,
445 in the format described in <locale.h>. */
446 const char *grouping;
448 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
449 struct locale_data *current = loc->__locales[LC_NUMERIC];
450 #endif
452 if (group)
454 grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
455 if (*grouping <= 0 || *grouping == CHAR_MAX)
456 grouping = NULL;
457 else
459 /* Figure out the thousands separator character. */
460 if (mbtowc (&thousands, _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP),
461 strlen (_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP))) <= 0)
462 thousands = (wchar_t) *_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
463 if (thousands == L'\0')
464 grouping = NULL;
467 else
469 grouping = NULL;
470 thousands = L'\0';
473 /* Find the locale's decimal point character. */
474 if (mbtowc ((wchar_t *) &decimal, _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT),
475 strlen (_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT))) <= 0)
476 decimal = (wchar_t) *_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
477 assert (decimal != L'\0');
479 /* Prepare number representation. */
480 exponent = 0;
481 negative = 0;
482 bits = 0;
484 /* Parse string to get maximal legal prefix. We need the number of
485 characters of the integer part, the fractional part and the exponent. */
486 cp = nptr - 1;
487 /* Ignore leading white space. */
489 c = *++cp;
490 while (ISSPACE (c));
492 /* Get sign of the result. */
493 if (c == L_('-'))
495 negative = 1;
496 c = *++cp;
498 else if (c == L_('+'))
499 c = *++cp;
501 /* Return 0.0 if no legal string is found.
502 No character is used even if a sign was found. */
503 if ((c < L_('0') || c > L_('9'))
504 && ((wchar_t) c != decimal || cp[1] < L_('0') || cp[1] > L_('9')))
506 int matched = 0;
507 /* Check for `INF' or `INFINITY'. */
508 if (TOLOWER (c) == L_('i')
509 && ((STRNCASECMP (cp, L_("inf"), 3) == 0 && (matched = 3))
510 || (STRNCASECMP (cp, L_("infinity"), 8) == 0 && (matched = 8))))
512 /* Return +/- infinity. */
513 if (endptr != NULL)
514 *endptr = (STRING_TYPE *) (cp + matched);
516 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
519 if (TOLOWER (c) == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
521 /* Return NaN. */
522 FLOAT retval = NAN;
524 cp += 3;
526 /* Match `(n-char-sequence-digit)'. */
527 if (*cp == L_('('))
529 const STRING_TYPE *startp = cp;
531 ++cp;
532 while ((*cp >= L_('0') && *cp <= L_('9'))
533 || (TOLOWER (*cp) >= L_('a') && TOLOWER (*cp) <= L_('z'))
534 || *cp == L_('_'));
536 if (*cp != L_(')'))
537 /* The closing brace is missing. Only match the NAN
538 part. */
539 cp = startp;
540 else
542 /* This is a system-dependent way to specify the
543 bitmask used for the NaN. We expect it to be
544 a number which is put in the mantissa of the
545 number. */
546 STRING_TYPE *endp;
547 unsigned long long int mant;
549 mant = STRTOULL (startp + 1, &endp, 0);
550 if (endp == cp)
551 SET_MANTISSA (retval, mant);
555 if (endptr != NULL)
556 *endptr = (STRING_TYPE *) cp;
558 return retval;
561 /* It is really a text we do not recognize. */
562 RETURN (0.0, nptr);
565 /* First look whether we are faced with a hexadecimal number. */
566 if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
568 /* Okay, it is a hexa-decimal number. Remember this and skip
569 the characters. BTW: hexadecimal numbers must not be
570 grouped. */
571 base = 16;
572 cp += 2;
573 c = *cp;
574 grouping = NULL;
577 /* Record the start of the digits, in case we will check their grouping. */
578 start_of_digits = startp = cp;
580 /* Ignore leading zeroes. This helps us to avoid useless computations. */
581 while (c == L_('0') || (thousands != L'\0' && (wchar_t) c == thousands))
582 c = *++cp;
584 /* If no other digit but a '0' is found the result is 0.0.
585 Return current read pointer. */
586 if ((c < L_('0') || c > L_('9')) &&
587 (base == 16 && (c < TOLOWER (L_('a')) || c > TOLOWER (L_('f')))) &&
588 (wchar_t) c != decimal &&
589 (base == 16 && (cp == start_of_digits || TOLOWER (c) != L_('p'))) &&
590 (base != 16 && TOLOWER (c) != L_('e')))
592 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
593 /* If TP is at the start of the digits, there was no correctly
594 grouped prefix of the string; so no number found. */
595 RETURN (0.0, tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
598 /* Remember first significant digit and read following characters until the
599 decimal point, exponent character or any non-FP number character. */
600 startp = cp;
601 dig_no = 0;
602 while (dig_no < (base == 16 ? HEXNDIG : NDIG) ||
603 /* If parsing grouping info, keep going past useful digits
604 so we can check all the grouping separators. */
605 grouping)
607 if ((c >= L_('0') && c <= L_('9'))
608 || (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f')))
609 ++dig_no;
610 else if (thousands == L'\0' || (wchar_t) c != thousands)
611 /* Not a digit or separator: end of the integer part. */
612 break;
613 c = *++cp;
616 if (grouping && dig_no > 0)
618 /* Check the grouping of the digits. */
619 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
620 if (cp != tp)
622 /* Less than the entire string was correctly grouped. */
624 if (tp == start_of_digits)
625 /* No valid group of numbers at all: no valid number. */
626 RETURN (0.0, nptr);
628 if (tp < startp)
629 /* The number is validly grouped, but consists
630 only of zeroes. The whole value is zero. */
631 RETURN (0.0, tp);
633 /* Recompute DIG_NO so we won't read more digits than
634 are properly grouped. */
635 cp = tp;
636 dig_no = 0;
637 for (tp = startp; tp < cp; ++tp)
638 if (*tp >= L_('0') && *tp <= L_('9'))
639 ++dig_no;
641 int_no = dig_no;
642 lead_zero = 0;
644 goto number_parsed;
648 if (dig_no >= (base == 16 ? HEXNDIG : NDIG))
649 /* Too many digits to be representable. Assigning this to EXPONENT
650 allows us to read the full number but return HUGE_VAL after parsing. */
651 exponent = MAX_10_EXP;
653 /* We have the number digits in the integer part. Whether these are all or
654 any is really a fractional digit will be decided later. */
655 int_no = dig_no;
656 lead_zero = int_no == 0 ? -1 : 0;
658 /* Read the fractional digits. A special case are the 'american style'
659 numbers like `16.' i.e. with decimal but without trailing digits. */
660 if ((wchar_t) c == decimal)
662 c = *++cp;
663 while (c >= L_('0') && c <= L_('9') ||
664 (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f')))
666 if (c != L_('0') && lead_zero == -1)
667 lead_zero = dig_no - int_no;
668 ++dig_no;
669 c = *++cp;
673 /* Remember start of exponent (if any). */
674 expp = cp;
676 /* Read exponent. */
677 if ((base == 16 && TOLOWER (c) == L_('p'))
678 || (base != 16 && TOLOWER (c) == L_('e')))
680 int exp_negative = 0;
682 c = *++cp;
683 if (c == L_('-'))
685 exp_negative = 1;
686 c = *++cp;
688 else if (c == L_('+'))
689 c = *++cp;
691 if (c >= L_('0') && c <= L_('9'))
693 int exp_limit;
695 /* Get the exponent limit. */
696 if (base == 16)
697 exp_limit = (exp_negative ?
698 -MIN_EXP + MANT_DIG - 4 * int_no :
699 MAX_EXP - 4 * int_no + lead_zero);
700 else
701 exp_limit = (exp_negative ?
702 -MIN_10_EXP + MANT_DIG - int_no :
703 MAX_10_EXP - int_no + lead_zero);
707 exponent *= 10;
709 if (exponent > exp_limit)
710 /* The exponent is too large/small to represent a valid
711 number. */
713 FLOAT result;
715 /* We have to take care for special situation: a joker
716 might have written "0.0e100000" which is in fact
717 zero. */
718 if (lead_zero == -1)
719 result = negative ? -0.0 : 0.0;
720 else
722 /* Overflow or underflow. */
723 __set_errno (ERANGE);
724 result = (exp_negative ? 0.0 :
725 negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
728 /* Accept all following digits as part of the exponent. */
730 ++cp;
731 while (*cp >= L_('0') && *cp <= L_('9'));
733 RETURN (result, cp);
734 /* NOTREACHED */
737 exponent += c - L_('0');
738 c = *++cp;
740 while (c >= L_('0') && c <= L_('9'));
742 if (exp_negative)
743 exponent = -exponent;
745 else
746 cp = expp;
749 /* We don't want to have to work with trailing zeroes after the radix. */
750 if (dig_no > int_no)
752 while (expp[-1] == L_('0'))
754 --expp;
755 --dig_no;
757 assert (dig_no >= int_no);
760 number_parsed:
762 /* The whole string is parsed. Store the address of the next character. */
763 if (endptr)
764 *endptr = (STRING_TYPE *) cp;
766 if (dig_no == 0)
767 return negative ? -0.0 : 0.0;
769 if (lead_zero)
771 /* Find the decimal point */
772 while ((wchar_t) *startp != decimal)
773 ++startp;
774 startp += lead_zero + 1;
775 exponent -= base == 16 ? 4 * lead_zero : lead_zero;
776 dig_no -= lead_zero;
779 /* If the BASE is 16 we can use a simpler algorithm. */
780 if (base == 16)
782 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
783 4, 4, 4, 4, 4, 4, 4, 4 };
784 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
785 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
786 mp_limb_t val;
788 while (!ISXDIGIT (*startp))
789 ++startp;
790 if (ISDIGIT (*startp))
791 val = *startp++ - L_('0');
792 else
793 val = 10 + TOLOWER (*startp++) - L_('a');
794 bits = nbits[val];
796 if (pos + 1 >= 4 || pos + 1 >= bits)
798 /* We don't have to care for wrapping. This is the normal
799 case so we add the first clause in the `if' expression as
800 an optimization. It is a compile-time constant and so does
801 not cost anything. */
802 retval[idx] = val << (pos - bits + 1);
803 pos -= bits;
805 else
807 retval[idx--] = val >> (bits - pos - 1);
808 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
809 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
812 while (--dig_no > 0 && idx >= 0)
814 while (!ISXDIGIT (*startp))
815 ++startp;
816 if (ISDIGIT (*startp))
817 val = *startp++ - L_('0');
818 else
819 val = 10 + TOLOWER (*startp++) - L_('a');
821 if (pos + 1 >= 4)
823 retval[idx] |= val << (pos - 4 + 1);
824 pos -= 4;
826 else
828 retval[idx--] |= val >> (4 - pos - 1);
829 val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
830 if (idx < 0)
831 return round_and_return (retval, exponent, negative, val,
832 BITS_PER_MP_LIMB - 1, dig_no > 0);
834 retval[idx] = val;
835 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
839 /* We ran out of digits. */
840 MPN_ZERO (retval, idx);
842 return round_and_return (retval, exponent, negative, 0, 0, 0);
845 /* Now we have the number of digits in total and the integer digits as well
846 as the exponent and its sign. We can decide whether the read digits are
847 really integer digits or belong to the fractional part; i.e. we normalize
848 123e-2 to 1.23. */
850 register int incr = (exponent < 0 ? MAX (-int_no, exponent)
851 : MIN (dig_no - int_no, exponent));
852 int_no += incr;
853 exponent -= incr;
856 if (int_no + exponent > MAX_10_EXP + 1)
858 __set_errno (ERANGE);
859 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
862 if (exponent < MIN_10_EXP - (DIG + 1))
864 __set_errno (ERANGE);
865 return 0.0;
868 if (int_no > 0)
870 /* Read the integer part as a multi-precision number to NUM. */
871 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent);
873 if (exponent > 0)
875 /* We now multiply the gained number by the given power of ten. */
876 mp_limb_t *psrc = num;
877 mp_limb_t *pdest = den;
878 int expbit = 1;
879 const struct mp_power *ttab = &_fpioconst_pow10[0];
883 if ((exponent & expbit) != 0)
885 mp_limb_t cy;
886 exponent ^= expbit;
888 /* FIXME: not the whole multiplication has to be
889 done. If we have the needed number of bits we
890 only need the information whether more non-zero
891 bits follow. */
892 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
893 cy = __mpn_mul (pdest, psrc, numsize,
894 &ttab->array[_FPIO_CONST_OFFSET],
895 ttab->arraysize - _FPIO_CONST_OFFSET);
896 else
897 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
898 ttab->arraysize - _FPIO_CONST_OFFSET,
899 psrc, numsize);
900 numsize += ttab->arraysize - _FPIO_CONST_OFFSET;
901 if (cy == 0)
902 --numsize;
903 SWAP (psrc, pdest);
905 expbit <<= 1;
906 ++ttab;
908 while (exponent != 0);
910 if (psrc == den)
911 memcpy (num, den, numsize * sizeof (mp_limb_t));
914 /* Determine how many bits of the result we already have. */
915 count_leading_zeros (bits, num[numsize - 1]);
916 bits = numsize * BITS_PER_MP_LIMB - bits;
918 /* Now we know the exponent of the number in base two.
919 Check it against the maximum possible exponent. */
920 if (bits > MAX_EXP)
922 __set_errno (ERANGE);
923 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
926 /* We have already the first BITS bits of the result. Together with
927 the information whether more non-zero bits follow this is enough
928 to determine the result. */
929 if (bits > MANT_DIG)
931 int i;
932 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
933 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
934 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
935 : least_idx;
936 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
937 : least_bit - 1;
939 if (least_bit == 0)
940 memcpy (retval, &num[least_idx],
941 RETURN_LIMB_SIZE * sizeof (mp_limb_t));
942 else
944 for (i = least_idx; i < numsize - 1; ++i)
945 retval[i - least_idx] = (num[i] >> least_bit)
946 | (num[i + 1]
947 << (BITS_PER_MP_LIMB - least_bit));
948 if (i - least_idx < RETURN_LIMB_SIZE)
949 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
952 /* Check whether any limb beside the ones in RETVAL are non-zero. */
953 for (i = 0; num[i] == 0; ++i)
956 return round_and_return (retval, bits - 1, negative,
957 num[round_idx], round_bit,
958 int_no < dig_no || i < round_idx);
959 /* NOTREACHED */
961 else if (dig_no == int_no)
963 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
964 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
966 if (target_bit == is_bit)
968 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
969 numsize * sizeof (mp_limb_t));
970 /* FIXME: the following loop can be avoided if we assume a
971 maximal MANT_DIG value. */
972 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
974 else if (target_bit > is_bit)
976 (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
977 num, numsize, target_bit - is_bit);
978 /* FIXME: the following loop can be avoided if we assume a
979 maximal MANT_DIG value. */
980 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
982 else
984 mp_limb_t cy;
985 assert (numsize < RETURN_LIMB_SIZE);
987 cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
988 num, numsize, is_bit - target_bit);
989 retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
990 /* FIXME: the following loop can be avoided if we assume a
991 maximal MANT_DIG value. */
992 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
995 return round_and_return (retval, bits - 1, negative, 0, 0, 0);
996 /* NOTREACHED */
999 /* Store the bits we already have. */
1000 memcpy (retval, num, numsize * sizeof (mp_limb_t));
1001 #if RETURN_LIMB_SIZE > 1
1002 if (numsize < RETURN_LIMB_SIZE)
1003 retval[numsize] = 0;
1004 #endif
1007 /* We have to compute at least some of the fractional digits. */
1009 /* We construct a fraction and the result of the division gives us
1010 the needed digits. The denominator is 1.0 multiplied by the
1011 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1012 123e-6 gives 123 / 1000000. */
1014 int expbit;
1015 int cnt;
1016 int neg_exp;
1017 int more_bits;
1018 mp_limb_t cy;
1019 mp_limb_t *psrc = den;
1020 mp_limb_t *pdest = num;
1021 const struct mp_power *ttab = &_fpioconst_pow10[0];
1023 assert (dig_no > int_no && exponent <= 0);
1026 /* For the fractional part we need not process too many digits. One
1027 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1028 ceil(BITS / 3) =: N
1029 digits we should have enough bits for the result. The remaining
1030 decimal digits give us the information that more bits are following.
1031 This can be used while rounding. (One added as a safety margin.) */
1032 if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
1034 dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
1035 more_bits = 1;
1037 else
1038 more_bits = 0;
1040 neg_exp = dig_no - int_no - exponent;
1042 /* Construct the denominator. */
1043 densize = 0;
1044 expbit = 1;
1047 if ((neg_exp & expbit) != 0)
1049 mp_limb_t cy;
1050 neg_exp ^= expbit;
1052 if (densize == 0)
1054 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
1055 memcpy (psrc, &ttab->array[_FPIO_CONST_OFFSET],
1056 densize * sizeof (mp_limb_t));
1058 else
1060 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
1061 ttab->arraysize - _FPIO_CONST_OFFSET,
1062 psrc, densize);
1063 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
1064 if (cy == 0)
1065 --densize;
1066 SWAP (psrc, pdest);
1069 expbit <<= 1;
1070 ++ttab;
1072 while (neg_exp != 0);
1074 if (psrc == num)
1075 memcpy (den, num, densize * sizeof (mp_limb_t));
1077 /* Read the fractional digits from the string. */
1078 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent);
1081 /* We now have to shift both numbers so that the highest bit in the
1082 denominator is set. In the same process we copy the numerator to
1083 a high place in the array so that the division constructs the wanted
1084 digits. This is done by a "quasi fix point" number representation.
1086 num: ddddddddddd . 0000000000000000000000
1087 |--- m ---|
1088 den: ddddddddddd n >= m
1089 |--- n ---|
1092 count_leading_zeros (cnt, den[densize - 1]);
1094 if (cnt > 0)
1096 /* Don't call `mpn_shift' with a count of zero since the specification
1097 does not allow this. */
1098 (void) __mpn_lshift (den, den, densize, cnt);
1099 cy = __mpn_lshift (num, num, numsize, cnt);
1100 if (cy != 0)
1101 num[numsize++] = cy;
1104 /* Now we are ready for the division. But it is not necessary to
1105 do a full multi-precision division because we only need a small
1106 number of bits for the result. So we do not use __mpn_divmod
1107 here but instead do the division here by hand and stop whenever
1108 the needed number of bits is reached. The code itself comes
1109 from the GNU MP Library by Torbj\"orn Granlund. */
1111 exponent = bits;
1113 switch (densize)
1115 case 1:
1117 mp_limb_t d, n, quot;
1118 int used = 0;
1120 n = num[0];
1121 d = den[0];
1122 assert (numsize == 1 && n < d);
1126 udiv_qrnnd (quot, n, n, 0, d);
1128 #define got_limb \
1129 if (bits == 0) \
1131 register int cnt; \
1132 if (quot == 0) \
1133 cnt = BITS_PER_MP_LIMB; \
1134 else \
1135 count_leading_zeros (cnt, quot); \
1136 exponent -= cnt; \
1137 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1139 used = MANT_DIG + cnt; \
1140 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1141 bits = MANT_DIG + 1; \
1143 else \
1145 /* Note that we only clear the second element. */ \
1146 /* The conditional is determined at compile time. */ \
1147 if (RETURN_LIMB_SIZE > 1) \
1148 retval[1] = 0; \
1149 retval[0] = quot; \
1150 bits = -cnt; \
1153 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1154 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1155 quot); \
1156 else \
1158 used = MANT_DIG - bits; \
1159 if (used > 0) \
1160 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1162 bits += BITS_PER_MP_LIMB
1164 got_limb;
1166 while (bits <= MANT_DIG);
1168 return round_and_return (retval, exponent - 1, negative,
1169 quot, BITS_PER_MP_LIMB - 1 - used,
1170 more_bits || n != 0);
1172 case 2:
1174 mp_limb_t d0, d1, n0, n1;
1175 mp_limb_t quot = 0;
1176 int used = 0;
1178 d0 = den[0];
1179 d1 = den[1];
1181 if (numsize < densize)
1183 if (num[0] >= d1)
1185 /* The numerator of the number occupies fewer bits than
1186 the denominator but the one limb is bigger than the
1187 high limb of the numerator. */
1188 n1 = 0;
1189 n0 = num[0];
1191 else
1193 if (bits <= 0)
1194 exponent -= BITS_PER_MP_LIMB;
1195 else
1197 if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
1198 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1199 BITS_PER_MP_LIMB, 0);
1200 else
1202 used = MANT_DIG - bits;
1203 if (used > 0)
1204 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1206 bits += BITS_PER_MP_LIMB;
1208 n1 = num[0];
1209 n0 = 0;
1212 else
1214 n1 = num[1];
1215 n0 = num[0];
1218 while (bits <= MANT_DIG)
1220 mp_limb_t r;
1222 if (n1 == d1)
1224 /* QUOT should be either 111..111 or 111..110. We need
1225 special treatment of this rare case as normal division
1226 would give overflow. */
1227 quot = ~(mp_limb_t) 0;
1229 r = n0 + d1;
1230 if (r < d1) /* Carry in the addition? */
1232 add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1233 goto have_quot;
1235 n1 = d0 - (d0 != 0);
1236 n0 = -d0;
1238 else
1240 udiv_qrnnd (quot, r, n1, n0, d1);
1241 umul_ppmm (n1, n0, d0, quot);
1244 q_test:
1245 if (n1 > r || (n1 == r && n0 > 0))
1247 /* The estimated QUOT was too large. */
1248 --quot;
1250 sub_ddmmss (n1, n0, n1, n0, 0, d0);
1251 r += d1;
1252 if (r >= d1) /* If not carry, test QUOT again. */
1253 goto q_test;
1255 sub_ddmmss (n1, n0, r, 0, n1, n0);
1257 have_quot:
1258 got_limb;
1261 return round_and_return (retval, exponent - 1, negative,
1262 quot, BITS_PER_MP_LIMB - 1 - used,
1263 more_bits || n1 != 0 || n0 != 0);
1265 default:
1267 int i;
1268 mp_limb_t cy, dX, d1, n0, n1;
1269 mp_limb_t quot = 0;
1270 int used = 0;
1272 dX = den[densize - 1];
1273 d1 = den[densize - 2];
1275 /* The division does not work if the upper limb of the two-limb
1276 numerator is greater than the denominator. */
1277 if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1278 num[numsize++] = 0;
1280 if (numsize < densize)
1282 mp_size_t empty = densize - numsize;
1284 if (bits <= 0)
1286 register int i;
1287 for (i = numsize; i > 0; --i)
1288 num[i + empty] = num[i - 1];
1289 MPN_ZERO (num, empty + 1);
1290 exponent -= empty * BITS_PER_MP_LIMB;
1292 else
1294 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1296 /* We make a difference here because the compiler
1297 cannot optimize the `else' case that good and
1298 this reflects all currently used FLOAT types
1299 and GMP implementations. */
1300 register int i;
1301 #if RETURN_LIMB_SIZE <= 2
1302 assert (empty == 1);
1303 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1304 BITS_PER_MP_LIMB, 0);
1305 #else
1306 for (i = RETURN_LIMB_SIZE; i > empty; --i)
1307 retval[i] = retval[i - empty];
1308 #endif
1309 #if RETURN_LIMB_SIZE > 1
1310 retval[1] = 0;
1311 #endif
1312 for (i = numsize; i > 0; --i)
1313 num[i + empty] = num[i - 1];
1314 MPN_ZERO (num, empty + 1);
1316 else
1318 used = MANT_DIG - bits;
1319 if (used >= BITS_PER_MP_LIMB)
1321 register int i;
1322 (void) __mpn_lshift (&retval[used
1323 / BITS_PER_MP_LIMB],
1324 retval, RETURN_LIMB_SIZE,
1325 used % BITS_PER_MP_LIMB);
1326 for (i = used / BITS_PER_MP_LIMB; i >= 0; --i)
1327 retval[i] = 0;
1329 else if (used > 0)
1330 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1332 bits += empty * BITS_PER_MP_LIMB;
1335 else
1337 int i;
1338 assert (numsize == densize);
1339 for (i = numsize; i > 0; --i)
1340 num[i] = num[i - 1];
1343 den[densize] = 0;
1344 n0 = num[densize];
1346 while (bits <= MANT_DIG)
1348 if (n0 == dX)
1349 /* This might over-estimate QUOT, but it's probably not
1350 worth the extra code here to find out. */
1351 quot = ~(mp_limb_t) 0;
1352 else
1354 mp_limb_t r;
1356 udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1357 umul_ppmm (n1, n0, d1, quot);
1359 while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1361 --quot;
1362 r += dX;
1363 if (r < dX) /* I.e. "carry in previous addition?" */
1364 break;
1365 n1 -= n0 < d1;
1366 n0 -= d1;
1370 /* Possible optimization: We already have (q * n0) and (1 * n1)
1371 after the calculation of QUOT. Taking advantage of this, we
1372 could make this loop make two iterations less. */
1374 cy = __mpn_submul_1 (num, den, densize + 1, quot);
1376 if (num[densize] != cy)
1378 cy = __mpn_add_n (num, num, den, densize);
1379 assert (cy != 0);
1380 --quot;
1382 n0 = num[densize] = num[densize - 1];
1383 for (i = densize - 1; i > 0; --i)
1384 num[i] = num[i - 1];
1386 got_limb;
1389 for (i = densize; num[i] == 0 && i >= 0; --i)
1391 return round_and_return (retval, exponent - 1, negative,
1392 quot, BITS_PER_MP_LIMB - 1 - used,
1393 more_bits || i >= 0);
1398 /* NOTREACHED */
1401 /* External user entry point. */
1403 FLOAT
1404 #ifdef weak_function
1405 weak_function
1406 #endif
1407 STRTOF (nptr, endptr LOCALE_PARAM)
1408 const STRING_TYPE *nptr;
1409 STRING_TYPE **endptr;
1410 LOCALE_PARAM_DECL
1412 return INTERNAL (STRTOF) (nptr, endptr, 0 LOCALE_PARAM);