Update.
[glibc.git] / stdlib / strtod.c
blob154e2049d9f93dd894d4339811951383ecbdcb5b
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
91 #if defined _LIBC || defined HAVE_WCHAR_H
92 # include <wchar.h>
93 #endif
95 #ifdef USE_WIDE_CHAR
96 # include <wctype.h>
97 # define STRING_TYPE wchar_t
98 # define CHAR_TYPE wint_t
99 # define L_(Ch) L##Ch
100 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
101 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
102 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
103 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
104 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
105 # define STRNCASECMP(S1, S2, N) __wcsncasecmp_l ((S1), (S2), (N), loc)
106 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
107 # else
108 # define ISSPACE(Ch) iswspace (Ch)
109 # define ISDIGIT(Ch) iswdigit (Ch)
110 # define ISXDIGIT(Ch) iswxdigit (Ch)
111 # define TOLOWER(Ch) towlower (Ch)
112 # define STRNCASECMP(S1, S2, N) __wcsncasecmp ((S1), (S2), (N))
113 # define STRTOULL(S, E, B) __wcstoull_internal ((S), (E), (B), 0)
114 # endif
115 #else
116 # define STRING_TYPE char
117 # define CHAR_TYPE char
118 # define L_(Ch) Ch
119 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
120 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
121 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
122 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
123 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
124 # define STRNCASECMP(S1, S2, N) __strncasecmp_l ((S1), (S2), (N), loc)
125 # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc)
126 # else
127 # define ISSPACE(Ch) isspace (Ch)
128 # define ISDIGIT(Ch) isdigit (Ch)
129 # define ISXDIGIT(Ch) isxdigit (Ch)
130 # define TOLOWER(Ch) tolower (Ch)
131 # define STRNCASECMP(S1, S2, N) __strncasecmp ((S1), (S2), (N))
132 # define STRTOULL(S, E, B) __strtoull_internal ((S), (E), 0, (B))
133 # endif
134 #endif
137 /* Constants we need from float.h; select the set for the FLOAT precision. */
138 #define MANT_DIG PASTE(FLT,_MANT_DIG)
139 #define DIG PASTE(FLT,_DIG)
140 #define MAX_EXP PASTE(FLT,_MAX_EXP)
141 #define MIN_EXP PASTE(FLT,_MIN_EXP)
142 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
143 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
145 /* Extra macros required to get FLT expanded before the pasting. */
146 #define PASTE(a,b) PASTE1(a,b)
147 #define PASTE1(a,b) a##b
149 /* Function to construct a floating point number from an MP integer
150 containing the fraction bits, a base 2 exponent, and a sign flag. */
151 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
153 /* Definitions according to limb size used. */
154 #if BITS_PER_MP_LIMB == 32
155 # define MAX_DIG_PER_LIMB 9
156 # define MAX_FAC_PER_LIMB 1000000000UL
157 #elif BITS_PER_MP_LIMB == 64
158 # define MAX_DIG_PER_LIMB 19
159 # define MAX_FAC_PER_LIMB 10000000000000000000UL
160 #else
161 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
162 #endif
165 /* Local data structure. */
166 static const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] =
167 { 0, 10, 100,
168 1000, 10000, 100000,
169 1000000, 10000000, 100000000,
170 1000000000
171 #if BITS_PER_MP_LIMB > 32
172 , 10000000000U, 100000000000U,
173 1000000000000U, 10000000000000U, 100000000000000U,
174 1000000000000000U, 10000000000000000U, 100000000000000000U,
175 1000000000000000000U, 10000000000000000000U
176 #endif
177 #if BITS_PER_MP_LIMB > 64
178 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
179 #endif
182 #ifndef howmany
183 #define howmany(x,y) (((x)+((y)-1))/(y))
184 #endif
185 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
187 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
188 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
189 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
191 #define RETURN(val,end) \
192 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
193 return val; } while (0)
195 /* Maximum size necessary for mpn integers to hold floating point numbers. */
196 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
197 + 2)
198 /* Declare an mpn integer variable that big. */
199 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
200 /* Copy an mpn integer value. */
201 #define MPN_ASSIGN(dst, src) \
202 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
205 /* Return a floating point number of the needed type according to the given
206 multi-precision number after possible rounding. */
207 static inline FLOAT
208 round_and_return (mp_limb_t *retval, int exponent, int negative,
209 mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
211 if (exponent < MIN_EXP - 1)
213 mp_size_t shift = MIN_EXP - 1 - exponent;
215 if (shift > MANT_DIG)
217 __set_errno (EDOM);
218 return 0.0;
221 more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
222 if (shift == MANT_DIG)
223 /* This is a special case to handle the very seldom case where
224 the mantissa will be empty after the shift. */
226 int i;
228 round_limb = retval[RETURN_LIMB_SIZE - 1];
229 round_bit = BITS_PER_MP_LIMB - 1;
230 for (i = 0; i < RETURN_LIMB_SIZE; ++i)
231 more_bits |= retval[i] != 0;
232 MPN_ZERO (retval, RETURN_LIMB_SIZE);
234 else if (shift >= BITS_PER_MP_LIMB)
236 int i;
238 round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
239 round_bit = (shift - 1) % BITS_PER_MP_LIMB;
240 for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
241 more_bits |= retval[i] != 0;
242 more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
243 != 0);
245 (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
246 RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
247 shift % BITS_PER_MP_LIMB);
248 MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
249 shift / BITS_PER_MP_LIMB);
251 else if (shift > 0)
253 round_limb = retval[0];
254 round_bit = shift - 1;
255 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
257 exponent = MIN_EXP - 2;
260 if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
261 && (more_bits || (retval[0] & 1) != 0
262 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
264 mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
266 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
267 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
268 (retval[RETURN_LIMB_SIZE - 1]
269 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
271 ++exponent;
272 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
273 retval[RETURN_LIMB_SIZE - 1]
274 |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
276 else if (exponent == MIN_EXP - 2
277 && (retval[RETURN_LIMB_SIZE - 1]
278 & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
279 != 0)
280 /* The number was denormalized but now normalized. */
281 exponent = MIN_EXP - 1;
284 if (exponent > MAX_EXP)
285 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
287 return MPN2FLOAT (retval, exponent, negative);
291 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
292 into N. Return the size of the number limbs in NSIZE at the first
293 character od the string that is not part of the integer as the function
294 value. If the EXPONENT is small enough to be taken as an additional
295 factor for the resulting number (see code) multiply by it. */
296 static inline const STRING_TYPE *
297 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
298 int *exponent)
300 /* Number of digits for actual limb. */
301 int cnt = 0;
302 mp_limb_t low = 0;
303 mp_limb_t start;
305 *nsize = 0;
306 assert (digcnt > 0);
309 if (cnt == MAX_DIG_PER_LIMB)
311 if (*nsize == 0)
313 n[0] = low;
314 *nsize = 1;
316 else
318 mp_limb_t cy;
319 cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
320 cy += __mpn_add_1 (n, n, *nsize, low);
321 if (cy != 0)
323 n[*nsize] = cy;
324 ++(*nsize);
327 cnt = 0;
328 low = 0;
331 /* There might be thousands separators or radix characters in
332 the string. But these all can be ignored because we know the
333 format of the number is correct and we have an exact number
334 of characters to read. */
335 while (*str < L_('0') || *str > L_('9'))
336 ++str;
337 low = low * 10 + *str++ - L_('0');
338 ++cnt;
340 while (--digcnt > 0);
342 if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
344 low *= _tens_in_limb[*exponent];
345 start = _tens_in_limb[cnt + *exponent];
346 *exponent = 0;
348 else
349 start = _tens_in_limb[cnt];
351 if (*nsize == 0)
353 n[0] = low;
354 *nsize = 1;
356 else
358 mp_limb_t cy;
359 cy = __mpn_mul_1 (n, n, *nsize, start);
360 cy += __mpn_add_1 (n, n, *nsize, low);
361 if (cy != 0)
362 n[(*nsize)++] = cy;
365 return str;
369 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
370 with the COUNT most significant bits of LIMB.
372 Tege doesn't like this function so I have to write it here myself. :)
373 --drepper */
374 static inline void
375 __mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
376 mp_limb_t limb)
378 if (count == BITS_PER_MP_LIMB)
380 /* Optimize the case of shifting by exactly a word:
381 just copy words, with no actual bit-shifting. */
382 mp_size_t i;
383 for (i = size - 1; i > 0; --i)
384 ptr[i] = ptr[i - 1];
385 ptr[0] = limb;
387 else
389 (void) __mpn_lshift (ptr, ptr, size, count);
390 ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
395 #define INTERNAL(x) INTERNAL1(x)
396 #define INTERNAL1(x) __##x##_internal
398 /* This file defines a function to check for correct grouping. */
399 #include "grouping.h"
402 /* Return a floating point number with the value of the given string NPTR.
403 Set *ENDPTR to the character after the last used one. If the number is
404 smaller than the smallest representable number, set `errno' to ERANGE and
405 return 0.0. If the number is too big to be represented, set `errno' to
406 ERANGE and return HUGE_VAL with the appropriate sign. */
407 FLOAT
408 INTERNAL (STRTOF) (nptr, endptr, group LOCALE_PARAM)
409 const STRING_TYPE *nptr;
410 STRING_TYPE **endptr;
411 int group;
412 LOCALE_PARAM_DECL
414 int negative; /* The sign of the number. */
415 MPN_VAR (num); /* MP representation of the number. */
416 int exponent; /* Exponent of the number. */
418 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
419 int base = 10;
421 /* When we have to compute fractional digits we form a fraction with a
422 second multi-precision number (and we sometimes need a second for
423 temporary results). */
424 MPN_VAR (den);
426 /* Representation for the return value. */
427 mp_limb_t retval[RETURN_LIMB_SIZE];
428 /* Number of bits currently in result value. */
429 int bits;
431 /* Running pointer after the last character processed in the string. */
432 const STRING_TYPE *cp, *tp;
433 /* Start of significant part of the number. */
434 const STRING_TYPE *startp, *start_of_digits;
435 /* Points at the character following the integer and fractional digits. */
436 const STRING_TYPE *expp;
437 /* Total number of digit and number of digits in integer part. */
438 int dig_no, int_no, lead_zero;
439 /* Contains the last character read. */
440 CHAR_TYPE c;
442 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
443 there. So define it ourselves if it remains undefined. */
444 #ifndef _WINT_T
445 typedef unsigned int wint_t;
446 #endif
447 /* The radix character of the current locale. */
448 wchar_t decimal;
449 /* The thousands character of the current locale. */
450 wchar_t thousands = L'\0';
451 /* The numeric grouping specification of the current locale,
452 in the format described in <locale.h>. */
453 const char *grouping;
455 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
456 struct locale_data *current = loc->__locales[LC_NUMERIC];
457 #endif
459 if (group)
461 grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
462 if (*grouping <= 0 || *grouping == CHAR_MAX)
463 grouping = NULL;
464 else
466 /* Figure out the thousands separator character. */
467 thousands = __btowc (*_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP));
468 if (thousands == WEOF)
469 thousands = L'\0';
470 if (thousands == L'\0')
471 grouping = NULL;
474 else
475 grouping = NULL;
477 /* Find the locale's decimal point character. */
478 decimal = __btowc (*_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT));
479 if (decimal == WEOF)
480 decimal = L'.';
481 assert (decimal != L'\0');
483 /* Prepare number representation. */
484 exponent = 0;
485 negative = 0;
486 bits = 0;
488 /* Parse string to get maximal legal prefix. We need the number of
489 characters of the integer part, the fractional part and the exponent. */
490 cp = nptr - 1;
491 /* Ignore leading white space. */
493 c = *++cp;
494 while (ISSPACE (c));
496 /* Get sign of the result. */
497 if (c == L_('-'))
499 negative = 1;
500 c = *++cp;
502 else if (c == L_('+'))
503 c = *++cp;
505 /* Return 0.0 if no legal string is found.
506 No character is used even if a sign was found. */
507 if ((c < L_('0') || c > L_('9'))
508 && ((wchar_t) c != decimal || cp[1] < L_('0') || cp[1] > L_('9')))
510 int matched = 0;
511 /* Check for `INF' or `INFINITY'. */
512 if (TOLOWER (c) == L_('i')
513 && ((STRNCASECMP (cp, L_("inf"), 3) == 0 && (matched = 3))
514 || (STRNCASECMP (cp, L_("infinity"), 8) == 0 && (matched = 8))))
516 /* Return +/- infinity. */
517 if (endptr != NULL)
518 *endptr = (STRING_TYPE *) (cp + matched);
520 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
523 if (TOLOWER (c) == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
525 /* Return NaN. */
526 FLOAT retval = NAN;
528 cp += 3;
530 /* Match `(n-char-sequence-digit)'. */
531 if (*cp == L_('('))
533 const STRING_TYPE *startp = cp;
535 ++cp;
536 while ((*cp >= L_('0') && *cp <= L_('9'))
537 || (TOLOWER (*cp) >= L_('a') && TOLOWER (*cp) <= L_('z'))
538 || *cp == L_('_'));
540 if (*cp != L_(')'))
541 /* The closing brace is missing. Only match the NAN
542 part. */
543 cp = startp;
544 else
546 /* This is a system-dependent way to specify the
547 bitmask used for the NaN. We expect it to be
548 a number which is put in the mantissa of the
549 number. */
550 STRING_TYPE *endp;
551 unsigned long long int mant;
553 mant = STRTOULL (startp + 1, &endp, 0);
554 if (endp == cp)
555 SET_MANTISSA (retval, mant);
559 if (endptr != NULL)
560 *endptr = (STRING_TYPE *) cp;
562 return retval;
565 /* It is really a text we do not recognize. */
566 RETURN (0.0, nptr);
569 /* First look whether we are faced with a hexadecimal number. */
570 if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
572 /* Okay, it is a hexa-decimal number. Remember this and skip
573 the characters. BTW: hexadecimal numbers must not be
574 grouped. */
575 base = 16;
576 cp += 2;
577 c = *cp;
578 grouping = NULL;
581 /* Record the start of the digits, in case we will check their grouping. */
582 start_of_digits = startp = cp;
584 /* Ignore leading zeroes. This helps us to avoid useless computations. */
585 while (c == L_('0') || (thousands != L'\0' && (wchar_t) c == thousands))
586 c = *++cp;
588 /* If no other digit but a '0' is found the result is 0.0.
589 Return current read pointer. */
590 if ((c < L_('0') || c > L_('9')) &&
591 (base == 16 && (c < TOLOWER (L_('a')) || c > TOLOWER (L_('f')))) &&
592 (wchar_t) c != decimal &&
593 (base == 16 && (cp == start_of_digits || TOLOWER (c) != L_('p'))) &&
594 (base != 16 && TOLOWER (c) != L_('e')))
596 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
597 /* If TP is at the start of the digits, there was no correctly
598 grouped prefix of the string; so no number found. */
599 RETURN (0.0, tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
602 /* Remember first significant digit and read following characters until the
603 decimal point, exponent character or any non-FP number character. */
604 startp = cp;
605 dig_no = 0;
606 while (dig_no < (base == 16 ? HEXNDIG : NDIG) ||
607 /* If parsing grouping info, keep going past useful digits
608 so we can check all the grouping separators. */
609 grouping)
611 if ((c >= L_('0') && c <= L_('9'))
612 || (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f')))
613 ++dig_no;
614 else if (thousands == L'\0' || (wchar_t) c != thousands)
615 /* Not a digit or separator: end of the integer part. */
616 break;
617 c = *++cp;
620 if (grouping && dig_no > 0)
622 /* Check the grouping of the digits. */
623 tp = correctly_grouped_prefix (start_of_digits, cp, thousands, grouping);
624 if (cp != tp)
626 /* Less than the entire string was correctly grouped. */
628 if (tp == start_of_digits)
629 /* No valid group of numbers at all: no valid number. */
630 RETURN (0.0, nptr);
632 if (tp < startp)
633 /* The number is validly grouped, but consists
634 only of zeroes. The whole value is zero. */
635 RETURN (0.0, tp);
637 /* Recompute DIG_NO so we won't read more digits than
638 are properly grouped. */
639 cp = tp;
640 dig_no = 0;
641 for (tp = startp; tp < cp; ++tp)
642 if (*tp >= L_('0') && *tp <= L_('9'))
643 ++dig_no;
645 int_no = dig_no;
646 lead_zero = 0;
648 goto number_parsed;
652 if (dig_no >= (base == 16 ? HEXNDIG : NDIG))
653 /* Too many digits to be representable. Assigning this to EXPONENT
654 allows us to read the full number but return HUGE_VAL after parsing. */
655 exponent = MAX_10_EXP;
657 /* We have the number digits in the integer part. Whether these are all or
658 any is really a fractional digit will be decided later. */
659 int_no = dig_no;
660 lead_zero = int_no == 0 ? -1 : 0;
662 /* Read the fractional digits. A special case are the 'american style'
663 numbers like `16.' i.e. with decimal but without trailing digits. */
664 if ((wchar_t) c == decimal)
666 c = *++cp;
667 while (c >= L_('0') && c <= L_('9') ||
668 (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f')))
670 if (c != L_('0') && lead_zero == -1)
671 lead_zero = dig_no - int_no;
672 ++dig_no;
673 c = *++cp;
677 /* Remember start of exponent (if any). */
678 expp = cp;
680 /* Read exponent. */
681 if ((base == 16 && TOLOWER (c) == L_('p'))
682 || (base != 16 && TOLOWER (c) == L_('e')))
684 int exp_negative = 0;
686 c = *++cp;
687 if (c == L_('-'))
689 exp_negative = 1;
690 c = *++cp;
692 else if (c == L_('+'))
693 c = *++cp;
695 if (c >= L_('0') && c <= L_('9'))
697 int exp_limit;
699 /* Get the exponent limit. */
700 if (base == 16)
701 exp_limit = (exp_negative ?
702 -MIN_EXP + MANT_DIG - 4 * int_no :
703 MAX_EXP - 4 * int_no + lead_zero);
704 else
705 exp_limit = (exp_negative ?
706 -MIN_10_EXP + MANT_DIG - int_no :
707 MAX_10_EXP - int_no + lead_zero);
711 exponent *= 10;
713 if (exponent > exp_limit)
714 /* The exponent is too large/small to represent a valid
715 number. */
717 FLOAT result;
719 /* We have to take care for special situation: a joker
720 might have written "0.0e100000" which is in fact
721 zero. */
722 if (lead_zero == -1)
723 result = negative ? -0.0 : 0.0;
724 else
726 /* Overflow or underflow. */
727 __set_errno (ERANGE);
728 result = (exp_negative ? 0.0 :
729 negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
732 /* Accept all following digits as part of the exponent. */
734 ++cp;
735 while (*cp >= L_('0') && *cp <= L_('9'));
737 RETURN (result, cp);
738 /* NOTREACHED */
741 exponent += c - L_('0');
742 c = *++cp;
744 while (c >= L_('0') && c <= L_('9'));
746 if (exp_negative)
747 exponent = -exponent;
749 else
750 cp = expp;
753 /* We don't want to have to work with trailing zeroes after the radix. */
754 if (dig_no > int_no)
756 while (expp[-1] == L_('0'))
758 --expp;
759 --dig_no;
761 assert (dig_no >= int_no);
764 number_parsed:
766 /* The whole string is parsed. Store the address of the next character. */
767 if (endptr)
768 *endptr = (STRING_TYPE *) cp;
770 if (dig_no == 0)
771 return negative ? -0.0 : 0.0;
773 if (lead_zero)
775 /* Find the decimal point */
776 while ((wchar_t) *startp != decimal)
777 ++startp;
778 startp += lead_zero + 1;
779 exponent -= base == 16 ? 4 * lead_zero : lead_zero;
780 dig_no -= lead_zero;
783 /* If the BASE is 16 we can use a simpler algorithm. */
784 if (base == 16)
786 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
787 4, 4, 4, 4, 4, 4, 4, 4 };
788 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
789 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
790 mp_limb_t val;
792 while (!ISXDIGIT (*startp))
793 ++startp;
794 if (ISDIGIT (*startp))
795 val = *startp++ - L_('0');
796 else
797 val = 10 + TOLOWER (*startp++) - L_('a');
798 bits = nbits[val];
800 if (pos + 1 >= 4 || pos + 1 >= bits)
802 /* We don't have to care for wrapping. This is the normal
803 case so we add the first clause in the `if' expression as
804 an optimization. It is a compile-time constant and so does
805 not cost anything. */
806 retval[idx] = val << (pos - bits + 1);
807 pos -= bits;
809 else
811 retval[idx--] = val >> (bits - pos - 1);
812 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
813 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
816 while (--dig_no > 0 && idx >= 0)
818 while (!ISXDIGIT (*startp))
819 ++startp;
820 if (ISDIGIT (*startp))
821 val = *startp++ - L_('0');
822 else
823 val = 10 + TOLOWER (*startp++) - L_('a');
825 if (pos + 1 >= 4)
827 retval[idx] |= val << (pos - 4 + 1);
828 pos -= 4;
830 else
832 retval[idx--] |= val >> (4 - pos - 1);
833 val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
834 if (idx < 0)
835 return round_and_return (retval, exponent, negative, val,
836 BITS_PER_MP_LIMB - 1, dig_no > 0);
838 retval[idx] = val;
839 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
843 /* We ran out of digits. */
844 MPN_ZERO (retval, idx);
846 return round_and_return (retval, exponent, negative, 0, 0, 0);
849 /* Now we have the number of digits in total and the integer digits as well
850 as the exponent and its sign. We can decide whether the read digits are
851 really integer digits or belong to the fractional part; i.e. we normalize
852 123e-2 to 1.23. */
854 register int incr = (exponent < 0 ? MAX (-int_no, exponent)
855 : MIN (dig_no - int_no, exponent));
856 int_no += incr;
857 exponent -= incr;
860 if (int_no + exponent > MAX_10_EXP + 1)
862 __set_errno (ERANGE);
863 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
866 if (exponent < MIN_10_EXP - (DIG + 1))
868 __set_errno (ERANGE);
869 return 0.0;
872 if (int_no > 0)
874 /* Read the integer part as a multi-precision number to NUM. */
875 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent);
877 if (exponent > 0)
879 /* We now multiply the gained number by the given power of ten. */
880 mp_limb_t *psrc = num;
881 mp_limb_t *pdest = den;
882 int expbit = 1;
883 const struct mp_power *ttab = &_fpioconst_pow10[0];
887 if ((exponent & expbit) != 0)
889 mp_limb_t cy;
890 exponent ^= expbit;
892 /* FIXME: not the whole multiplication has to be
893 done. If we have the needed number of bits we
894 only need the information whether more non-zero
895 bits follow. */
896 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
897 cy = __mpn_mul (pdest, psrc, numsize,
898 &ttab->array[_FPIO_CONST_OFFSET],
899 ttab->arraysize - _FPIO_CONST_OFFSET);
900 else
901 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
902 ttab->arraysize - _FPIO_CONST_OFFSET,
903 psrc, numsize);
904 numsize += ttab->arraysize - _FPIO_CONST_OFFSET;
905 if (cy == 0)
906 --numsize;
907 SWAP (psrc, pdest);
909 expbit <<= 1;
910 ++ttab;
912 while (exponent != 0);
914 if (psrc == den)
915 memcpy (num, den, numsize * sizeof (mp_limb_t));
918 /* Determine how many bits of the result we already have. */
919 count_leading_zeros (bits, num[numsize - 1]);
920 bits = numsize * BITS_PER_MP_LIMB - bits;
922 /* Now we know the exponent of the number in base two.
923 Check it against the maximum possible exponent. */
924 if (bits > MAX_EXP)
926 __set_errno (ERANGE);
927 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
930 /* We have already the first BITS bits of the result. Together with
931 the information whether more non-zero bits follow this is enough
932 to determine the result. */
933 if (bits > MANT_DIG)
935 int i;
936 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
937 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
938 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
939 : least_idx;
940 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
941 : least_bit - 1;
943 if (least_bit == 0)
944 memcpy (retval, &num[least_idx],
945 RETURN_LIMB_SIZE * sizeof (mp_limb_t));
946 else
948 for (i = least_idx; i < numsize - 1; ++i)
949 retval[i - least_idx] = (num[i] >> least_bit)
950 | (num[i + 1]
951 << (BITS_PER_MP_LIMB - least_bit));
952 if (i - least_idx < RETURN_LIMB_SIZE)
953 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
956 /* Check whether any limb beside the ones in RETVAL are non-zero. */
957 for (i = 0; num[i] == 0; ++i)
960 return round_and_return (retval, bits - 1, negative,
961 num[round_idx], round_bit,
962 int_no < dig_no || i < round_idx);
963 /* NOTREACHED */
965 else if (dig_no == int_no)
967 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
968 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
970 if (target_bit == is_bit)
972 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
973 numsize * sizeof (mp_limb_t));
974 /* FIXME: the following loop can be avoided if we assume a
975 maximal MANT_DIG value. */
976 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
978 else if (target_bit > is_bit)
980 (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
981 num, numsize, target_bit - is_bit);
982 /* FIXME: the following loop can be avoided if we assume a
983 maximal MANT_DIG value. */
984 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
986 else
988 mp_limb_t cy;
989 assert (numsize < RETURN_LIMB_SIZE);
991 cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
992 num, numsize, is_bit - target_bit);
993 retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
994 /* FIXME: the following loop can be avoided if we assume a
995 maximal MANT_DIG value. */
996 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
999 return round_and_return (retval, bits - 1, negative, 0, 0, 0);
1000 /* NOTREACHED */
1003 /* Store the bits we already have. */
1004 memcpy (retval, num, numsize * sizeof (mp_limb_t));
1005 #if RETURN_LIMB_SIZE > 1
1006 if (numsize < RETURN_LIMB_SIZE)
1007 retval[numsize] = 0;
1008 #endif
1011 /* We have to compute at least some of the fractional digits. */
1013 /* We construct a fraction and the result of the division gives us
1014 the needed digits. The denominator is 1.0 multiplied by the
1015 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1016 123e-6 gives 123 / 1000000. */
1018 int expbit;
1019 int cnt;
1020 int neg_exp;
1021 int more_bits;
1022 mp_limb_t cy;
1023 mp_limb_t *psrc = den;
1024 mp_limb_t *pdest = num;
1025 const struct mp_power *ttab = &_fpioconst_pow10[0];
1027 assert (dig_no > int_no && exponent <= 0);
1030 /* For the fractional part we need not process too many digits. One
1031 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1032 ceil(BITS / 3) =: N
1033 digits we should have enough bits for the result. The remaining
1034 decimal digits give us the information that more bits are following.
1035 This can be used while rounding. (One added as a safety margin.) */
1036 if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1)
1038 dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1;
1039 more_bits = 1;
1041 else
1042 more_bits = 0;
1044 neg_exp = dig_no - int_no - exponent;
1046 /* Construct the denominator. */
1047 densize = 0;
1048 expbit = 1;
1051 if ((neg_exp & expbit) != 0)
1053 mp_limb_t cy;
1054 neg_exp ^= expbit;
1056 if (densize == 0)
1058 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
1059 memcpy (psrc, &ttab->array[_FPIO_CONST_OFFSET],
1060 densize * sizeof (mp_limb_t));
1062 else
1064 cy = __mpn_mul (pdest, &ttab->array[_FPIO_CONST_OFFSET],
1065 ttab->arraysize - _FPIO_CONST_OFFSET,
1066 psrc, densize);
1067 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
1068 if (cy == 0)
1069 --densize;
1070 SWAP (psrc, pdest);
1073 expbit <<= 1;
1074 ++ttab;
1076 while (neg_exp != 0);
1078 if (psrc == num)
1079 memcpy (den, num, densize * sizeof (mp_limb_t));
1081 /* Read the fractional digits from the string. */
1082 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent);
1085 /* We now have to shift both numbers so that the highest bit in the
1086 denominator is set. In the same process we copy the numerator to
1087 a high place in the array so that the division constructs the wanted
1088 digits. This is done by a "quasi fix point" number representation.
1090 num: ddddddddddd . 0000000000000000000000
1091 |--- m ---|
1092 den: ddddddddddd n >= m
1093 |--- n ---|
1096 count_leading_zeros (cnt, den[densize - 1]);
1098 if (cnt > 0)
1100 /* Don't call `mpn_shift' with a count of zero since the specification
1101 does not allow this. */
1102 (void) __mpn_lshift (den, den, densize, cnt);
1103 cy = __mpn_lshift (num, num, numsize, cnt);
1104 if (cy != 0)
1105 num[numsize++] = cy;
1108 /* Now we are ready for the division. But it is not necessary to
1109 do a full multi-precision division because we only need a small
1110 number of bits for the result. So we do not use __mpn_divmod
1111 here but instead do the division here by hand and stop whenever
1112 the needed number of bits is reached. The code itself comes
1113 from the GNU MP Library by Torbj\"orn Granlund. */
1115 exponent = bits;
1117 switch (densize)
1119 case 1:
1121 mp_limb_t d, n, quot;
1122 int used = 0;
1124 n = num[0];
1125 d = den[0];
1126 assert (numsize == 1 && n < d);
1130 udiv_qrnnd (quot, n, n, 0, d);
1132 #define got_limb \
1133 if (bits == 0) \
1135 register int cnt; \
1136 if (quot == 0) \
1137 cnt = BITS_PER_MP_LIMB; \
1138 else \
1139 count_leading_zeros (cnt, quot); \
1140 exponent -= cnt; \
1141 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1143 used = MANT_DIG + cnt; \
1144 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1145 bits = MANT_DIG + 1; \
1147 else \
1149 /* Note that we only clear the second element. */ \
1150 /* The conditional is determined at compile time. */ \
1151 if (RETURN_LIMB_SIZE > 1) \
1152 retval[1] = 0; \
1153 retval[0] = quot; \
1154 bits = -cnt; \
1157 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1158 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1159 quot); \
1160 else \
1162 used = MANT_DIG - bits; \
1163 if (used > 0) \
1164 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1166 bits += BITS_PER_MP_LIMB
1168 got_limb;
1170 while (bits <= MANT_DIG);
1172 return round_and_return (retval, exponent - 1, negative,
1173 quot, BITS_PER_MP_LIMB - 1 - used,
1174 more_bits || n != 0);
1176 case 2:
1178 mp_limb_t d0, d1, n0, n1;
1179 mp_limb_t quot = 0;
1180 int used = 0;
1182 d0 = den[0];
1183 d1 = den[1];
1185 if (numsize < densize)
1187 if (num[0] >= d1)
1189 /* The numerator of the number occupies fewer bits than
1190 the denominator but the one limb is bigger than the
1191 high limb of the numerator. */
1192 n1 = 0;
1193 n0 = num[0];
1195 else
1197 if (bits <= 0)
1198 exponent -= BITS_PER_MP_LIMB;
1199 else
1201 if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
1202 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1203 BITS_PER_MP_LIMB, 0);
1204 else
1206 used = MANT_DIG - bits;
1207 if (used > 0)
1208 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1210 bits += BITS_PER_MP_LIMB;
1212 n1 = num[0];
1213 n0 = 0;
1216 else
1218 n1 = num[1];
1219 n0 = num[0];
1222 while (bits <= MANT_DIG)
1224 mp_limb_t r;
1226 if (n1 == d1)
1228 /* QUOT should be either 111..111 or 111..110. We need
1229 special treatment of this rare case as normal division
1230 would give overflow. */
1231 quot = ~(mp_limb_t) 0;
1233 r = n0 + d1;
1234 if (r < d1) /* Carry in the addition? */
1236 add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1237 goto have_quot;
1239 n1 = d0 - (d0 != 0);
1240 n0 = -d0;
1242 else
1244 udiv_qrnnd (quot, r, n1, n0, d1);
1245 umul_ppmm (n1, n0, d0, quot);
1248 q_test:
1249 if (n1 > r || (n1 == r && n0 > 0))
1251 /* The estimated QUOT was too large. */
1252 --quot;
1254 sub_ddmmss (n1, n0, n1, n0, 0, d0);
1255 r += d1;
1256 if (r >= d1) /* If not carry, test QUOT again. */
1257 goto q_test;
1259 sub_ddmmss (n1, n0, r, 0, n1, n0);
1261 have_quot:
1262 got_limb;
1265 return round_and_return (retval, exponent - 1, negative,
1266 quot, BITS_PER_MP_LIMB - 1 - used,
1267 more_bits || n1 != 0 || n0 != 0);
1269 default:
1271 int i;
1272 mp_limb_t cy, dX, d1, n0, n1;
1273 mp_limb_t quot = 0;
1274 int used = 0;
1276 dX = den[densize - 1];
1277 d1 = den[densize - 2];
1279 /* The division does not work if the upper limb of the two-limb
1280 numerator is greater than the denominator. */
1281 if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1282 num[numsize++] = 0;
1284 if (numsize < densize)
1286 mp_size_t empty = densize - numsize;
1288 if (bits <= 0)
1290 register int i;
1291 for (i = numsize; i > 0; --i)
1292 num[i + empty] = num[i - 1];
1293 MPN_ZERO (num, empty + 1);
1294 exponent -= empty * BITS_PER_MP_LIMB;
1296 else
1298 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1300 /* We make a difference here because the compiler
1301 cannot optimize the `else' case that good and
1302 this reflects all currently used FLOAT types
1303 and GMP implementations. */
1304 register int i;
1305 #if RETURN_LIMB_SIZE <= 2
1306 assert (empty == 1);
1307 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1308 BITS_PER_MP_LIMB, 0);
1309 #else
1310 for (i = RETURN_LIMB_SIZE; i > empty; --i)
1311 retval[i] = retval[i - empty];
1312 #endif
1313 #if RETURN_LIMB_SIZE > 1
1314 retval[1] = 0;
1315 #endif
1316 for (i = numsize; i > 0; --i)
1317 num[i + empty] = num[i - 1];
1318 MPN_ZERO (num, empty + 1);
1320 else
1322 used = MANT_DIG - bits;
1323 if (used >= BITS_PER_MP_LIMB)
1325 register int i;
1326 (void) __mpn_lshift (&retval[used
1327 / BITS_PER_MP_LIMB],
1328 retval, RETURN_LIMB_SIZE,
1329 used % BITS_PER_MP_LIMB);
1330 for (i = used / BITS_PER_MP_LIMB; i >= 0; --i)
1331 retval[i] = 0;
1333 else if (used > 0)
1334 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1336 bits += empty * BITS_PER_MP_LIMB;
1339 else
1341 int i;
1342 assert (numsize == densize);
1343 for (i = numsize; i > 0; --i)
1344 num[i] = num[i - 1];
1347 den[densize] = 0;
1348 n0 = num[densize];
1350 while (bits <= MANT_DIG)
1352 if (n0 == dX)
1353 /* This might over-estimate QUOT, but it's probably not
1354 worth the extra code here to find out. */
1355 quot = ~(mp_limb_t) 0;
1356 else
1358 mp_limb_t r;
1360 udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1361 umul_ppmm (n1, n0, d1, quot);
1363 while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1365 --quot;
1366 r += dX;
1367 if (r < dX) /* I.e. "carry in previous addition?" */
1368 break;
1369 n1 -= n0 < d1;
1370 n0 -= d1;
1374 /* Possible optimization: We already have (q * n0) and (1 * n1)
1375 after the calculation of QUOT. Taking advantage of this, we
1376 could make this loop make two iterations less. */
1378 cy = __mpn_submul_1 (num, den, densize + 1, quot);
1380 if (num[densize] != cy)
1382 cy = __mpn_add_n (num, num, den, densize);
1383 assert (cy != 0);
1384 --quot;
1386 n0 = num[densize] = num[densize - 1];
1387 for (i = densize - 1; i > 0; --i)
1388 num[i] = num[i - 1];
1390 got_limb;
1393 for (i = densize; num[i] == 0 && i >= 0; --i)
1395 return round_and_return (retval, exponent - 1, negative,
1396 quot, BITS_PER_MP_LIMB - 1 - used,
1397 more_bits || i >= 0);
1402 /* NOTREACHED */
1405 /* External user entry point. */
1407 FLOAT
1408 #ifdef weak_function
1409 weak_function
1410 #endif
1411 STRTOF (nptr, endptr LOCALE_PARAM)
1412 const STRING_TYPE *nptr;
1413 STRING_TYPE **endptr;
1414 LOCALE_PARAM_DECL
1416 return INTERNAL (STRTOF) (nptr, endptr, 0 LOCALE_PARAM);