1 /* Read decimal floating point numbers.
2 This file is part of the GNU C Library.
3 Copyright (C) 1995-2002, 2003 Free Software Foundation, Inc.
4 Contributed by Ulrich Drepper <drepper@gnu.org>, 1995.
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the 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 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, write to the Free
18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
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. */
28 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
29 # define STRTOF __wcstod_l
31 # define STRTOF wcstod
34 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
35 # define STRTOF __strtod_l
37 # define STRTOF strtod
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; \
45 if ((mant & 0xfffffffffffffULL) == 0) \
46 mant = 0x8000000000000ULL; \
47 u.ieee.mantissa0 = ((mant) >> 32) & 0xfffff; \
48 u.ieee.mantissa1 = (mant) & 0xffffffff; \
52 /* End of configuration part. */
58 #include "../locale/localeinfo.h"
64 /* The gmp headers need some configuration frobs. */
67 /* Include gmp-mparam.h first, such that definitions of _SHORT_LIMB
68 and _LONG_LONG_LIMB in it can take effect into gmp.h. */
69 #include <gmp-mparam.h>
73 #include "fpioconst.h"
79 /* We use this code also for the extended locale handling where the
80 function gets as an additional argument the locale which has to be
81 used. To access the values we have to redefine the _NL_CURRENT
83 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
85 # define _NL_CURRENT(category, item) \
86 (current->values[_NL_ITEM_INDEX (item)].string)
87 # define LOCALE_PARAM , loc
88 # define LOCALE_PARAM_DECL __locale_t loc;
91 # define LOCALE_PARAM_DECL
94 #if defined _LIBC || defined HAVE_WCHAR_H
100 # define STRING_TYPE wchar_t
101 # define CHAR_TYPE wint_t
102 # define L_(Ch) L##Ch
103 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
104 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
105 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
106 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
107 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
108 # define STRNCASECMP(S1, S2, N) __wcsncasecmp_l ((S1), (S2), (N), loc)
109 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
111 # define ISSPACE(Ch) iswspace (Ch)
112 # define ISDIGIT(Ch) iswdigit (Ch)
113 # define ISXDIGIT(Ch) iswxdigit (Ch)
114 # define TOLOWER(Ch) towlower (Ch)
115 # define STRNCASECMP(S1, S2, N) __wcsncasecmp ((S1), (S2), (N))
116 # define STRTOULL(S, E, B) __wcstoull_internal ((S), (E), (B), 0)
119 # define STRING_TYPE char
120 # define CHAR_TYPE char
122 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
123 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
124 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
125 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
126 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
127 # define STRNCASECMP(S1, S2, N) __strncasecmp_l ((S1), (S2), (N), loc)
128 # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc)
130 # define ISSPACE(Ch) isspace (Ch)
131 # define ISDIGIT(Ch) isdigit (Ch)
132 # define ISXDIGIT(Ch) isxdigit (Ch)
133 # define TOLOWER(Ch) tolower (Ch)
134 # define STRNCASECMP(S1, S2, N) __strncasecmp ((S1), (S2), (N))
135 # define STRTOULL(S, E, B) __strtoull_internal ((S), (E), 0, (B))
140 /* Constants we need from float.h; select the set for the FLOAT precision. */
141 #define MANT_DIG PASTE(FLT,_MANT_DIG)
142 #define DIG PASTE(FLT,_DIG)
143 #define MAX_EXP PASTE(FLT,_MAX_EXP)
144 #define MIN_EXP PASTE(FLT,_MIN_EXP)
145 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
146 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
148 /* Extra macros required to get FLT expanded before the pasting. */
149 #define PASTE(a,b) PASTE1(a,b)
150 #define PASTE1(a,b) a##b
152 /* Function to construct a floating point number from an MP integer
153 containing the fraction bits, a base 2 exponent, and a sign flag. */
154 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
156 /* Definitions according to limb size used. */
157 #if BITS_PER_MP_LIMB == 32
158 # define MAX_DIG_PER_LIMB 9
159 # define MAX_FAC_PER_LIMB 1000000000UL
160 #elif BITS_PER_MP_LIMB == 64
161 # define MAX_DIG_PER_LIMB 19
162 # define MAX_FAC_PER_LIMB 10000000000000000000ULL
164 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
168 /* Local data structure. */
169 static const mp_limb_t _tens_in_limb
[MAX_DIG_PER_LIMB
+ 1] =
171 1000, 10000, 100000L,
172 1000000L, 10000000L, 100000000L,
174 #if BITS_PER_MP_LIMB > 32
175 , 10000000000ULL, 100000000000ULL,
176 1000000000000ULL, 10000000000000ULL, 100000000000000ULL,
177 1000000000000000ULL, 10000000000000000ULL, 100000000000000000ULL,
178 1000000000000000000ULL, 10000000000000000000ULL
180 #if BITS_PER_MP_LIMB > 64
181 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
186 #define howmany(x,y) (((x)+((y)-1))/(y))
188 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
190 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
191 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
192 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
194 #define RETURN(val,end) \
195 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
196 return val; } while (0)
198 /* Maximum size necessary for mpn integers to hold floating point numbers. */
199 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
201 /* Declare an mpn integer variable that big. */
202 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
203 /* Copy an mpn integer value. */
204 #define MPN_ASSIGN(dst, src) \
205 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
208 /* Return a floating point number of the needed type according to the given
209 multi-precision number after possible rounding. */
211 round_and_return (mp_limb_t
*retval
, int exponent
, int negative
,
212 mp_limb_t round_limb
, mp_size_t round_bit
, int more_bits
)
214 if (exponent
< MIN_EXP
- 1)
216 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
218 if (shift
> MANT_DIG
)
224 more_bits
|= (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0;
225 if (shift
== MANT_DIG
)
226 /* This is a special case to handle the very seldom case where
227 the mantissa will be empty after the shift. */
231 round_limb
= retval
[RETURN_LIMB_SIZE
- 1];
232 round_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
233 for (i
= 0; i
< RETURN_LIMB_SIZE
; ++i
)
234 more_bits
|= retval
[i
] != 0;
235 MPN_ZERO (retval
, RETURN_LIMB_SIZE
);
237 else if (shift
>= BITS_PER_MP_LIMB
)
241 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
242 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
243 for (i
= 0; i
< (shift
- 1) / BITS_PER_MP_LIMB
; ++i
)
244 more_bits
|= retval
[i
] != 0;
245 more_bits
|= ((round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1))
248 (void) __mpn_rshift (retval
, &retval
[shift
/ BITS_PER_MP_LIMB
],
249 RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
),
250 shift
% BITS_PER_MP_LIMB
);
251 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
)],
252 shift
/ BITS_PER_MP_LIMB
);
256 round_limb
= retval
[0];
257 round_bit
= shift
- 1;
258 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
260 /* This is a hook for the m68k long double format, where the
261 exponent bias is the same for normalized and denormalized
264 # define DENORM_EXP (MIN_EXP - 2)
266 exponent
= DENORM_EXP
;
269 if ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
270 && (more_bits
|| (retval
[0] & 1) != 0
271 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
273 mp_limb_t cy
= __mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
275 if (((MANT_DIG
% BITS_PER_MP_LIMB
) == 0 && cy
) ||
276 ((MANT_DIG
% BITS_PER_MP_LIMB
) != 0 &&
277 (retval
[RETURN_LIMB_SIZE
- 1]
278 & (((mp_limb_t
) 1) << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0))
281 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
282 retval
[RETURN_LIMB_SIZE
- 1]
283 |= ((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
);
285 else if (exponent
== DENORM_EXP
286 && (retval
[RETURN_LIMB_SIZE
- 1]
287 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
289 /* The number was denormalized but now normalized. */
290 exponent
= MIN_EXP
- 1;
293 if (exponent
> MAX_EXP
)
294 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
296 return MPN2FLOAT (retval
, exponent
, negative
);
300 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
301 into N. Return the size of the number limbs in NSIZE at the first
302 character od the string that is not part of the integer as the function
303 value. If the EXPONENT is small enough to be taken as an additional
304 factor for the resulting number (see code) multiply by it. */
305 static inline const STRING_TYPE
*
306 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
308 #ifndef USE_WIDE_CHAR
309 , const char *decimal
, size_t decimal_len
, const char *thousands
314 /* Number of digits for actual limb. */
323 if (cnt
== MAX_DIG_PER_LIMB
)
333 cy
= __mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
334 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
345 /* There might be thousands separators or radix characters in
346 the string. But these all can be ignored because we know the
347 format of the number is correct and we have an exact number
348 of characters to read. */
350 if (*str
< L
'0' || *str
> L
'9')
353 if (*str
< '0' || *str
> '9')
356 if (thousands
!= NULL
&& *str
== *thousands
357 && ({ for (inner
= 1; thousands
[inner
] != '\0'; ++inner
)
358 if (thousands
[inner
] != str
[inner
])
360 thousands
[inner
] == '\0'; }))
366 low
= low
* 10 + *str
++ - L_('0');
369 while (--digcnt
> 0);
371 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
373 low
*= _tens_in_limb
[*exponent
];
374 start
= _tens_in_limb
[cnt
+ *exponent
];
378 start
= _tens_in_limb
[cnt
];
388 cy
= __mpn_mul_1 (n
, n
, *nsize
, start
);
389 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
398 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
399 with the COUNT most significant bits of LIMB.
401 Tege doesn't like this function so I have to write it here myself. :)
404 __mpn_lshift_1 (mp_limb_t
*ptr
, mp_size_t size
, unsigned int count
,
407 if (count
== BITS_PER_MP_LIMB
)
409 /* Optimize the case of shifting by exactly a word:
410 just copy words, with no actual bit-shifting. */
412 for (i
= size
- 1; i
> 0; --i
)
418 (void) __mpn_lshift (ptr
, ptr
, size
, count
);
419 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
424 #define INTERNAL(x) INTERNAL1(x)
425 #define INTERNAL1(x) __##x##_internal
427 /* This file defines a function to check for correct grouping. */
428 #include "grouping.h"
431 /* Return a floating point number with the value of the given string NPTR.
432 Set *ENDPTR to the character after the last used one. If the number is
433 smaller than the smallest representable number, set `errno' to ERANGE and
434 return 0.0. If the number is too big to be represented, set `errno' to
435 ERANGE and return HUGE_VAL with the appropriate sign. */
437 INTERNAL (STRTOF
) (nptr
, endptr
, group LOCALE_PARAM
)
438 const STRING_TYPE
*nptr
;
439 STRING_TYPE
**endptr
;
443 int negative
; /* The sign of the number. */
444 MPN_VAR (num
); /* MP representation of the number. */
445 int exponent
; /* Exponent of the number. */
447 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
450 /* When we have to compute fractional digits we form a fraction with a
451 second multi-precision number (and we sometimes need a second for
452 temporary results). */
455 /* Representation for the return value. */
456 mp_limb_t retval
[RETURN_LIMB_SIZE
];
457 /* Number of bits currently in result value. */
460 /* Running pointer after the last character processed in the string. */
461 const STRING_TYPE
*cp
, *tp
;
462 /* Start of significant part of the number. */
463 const STRING_TYPE
*startp
, *start_of_digits
;
464 /* Points at the character following the integer and fractional digits. */
465 const STRING_TYPE
*expp
;
466 /* Total number of digit and number of digits in integer part. */
467 int dig_no
, int_no
, lead_zero
;
468 /* Contains the last character read. */
471 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
472 there. So define it ourselves if it remains undefined. */
474 typedef unsigned int wint_t;
476 /* The radix character of the current locale. */
483 /* The thousands character of the current locale. */
485 wchar_t thousands
= L
'\0';
487 const char *thousands
= NULL
;
489 /* The numeric grouping specification of the current locale,
490 in the format described in <locale.h>. */
491 const char *grouping
;
492 /* Used in several places. */
495 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
496 struct locale_data
*current
= loc
->__locales
[LC_NUMERIC
];
501 grouping
= _NL_CURRENT (LC_NUMERIC
, GROUPING
);
502 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
506 /* Figure out the thousands separator character. */
508 thousands
= _NL_CURRENT_WORD (LC_NUMERIC
,
509 _NL_NUMERIC_THOUSANDS_SEP_WC
);
510 if (thousands
== L
'\0')
513 thousands
= _NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
);
514 if (*thousands
== '\0')
525 /* Find the locale's decimal point character. */
527 decimal
= _NL_CURRENT_WORD (LC_NUMERIC
, _NL_NUMERIC_DECIMAL_POINT_WC
);
528 assert (decimal
!= L
'\0');
529 # define decimal_len 1
531 decimal
= _NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
);
532 decimal_len
= strlen (decimal
);
533 assert (decimal_len
> 0);
536 /* Prepare number representation. */
541 /* Parse string to get maximal legal prefix. We need the number of
542 characters of the integer part, the fractional part and the exponent. */
544 /* Ignore leading white space. */
549 /* Get sign of the result. */
555 else if (c
== L_('+'))
558 /* Return 0.0 if no legal string is found.
559 No character is used even if a sign was found. */
561 if (c
== (wint_t) decimal
562 && (wint_t) cp
[1] >= L
'0' && (wint_t) cp
[1] <= L
'9')
564 /* We accept it. This funny construct is here only to indent
565 the code directly. */
568 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
569 if (cp
[cnt
] != decimal
[cnt
])
571 if (decimal
[cnt
] == '\0' && cp
[cnt
] >= '0' && cp
[cnt
] <= '9')
573 /* We accept it. This funny construct is here only to indent
574 the code directly. */
577 else if (c
< L_('0') || c
> L_('9'))
579 /* Check for `INF' or `INFINITY'. */
580 if (TOLOWER (c
) == L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
582 /* Return +/- infinity. */
584 *endptr
= (STRING_TYPE
*)
585 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
588 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
591 if (TOLOWER (c
) == L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
598 /* Match `(n-char-sequence-digit)'. */
601 const STRING_TYPE
*startp
= cp
;
604 while ((*cp
>= L_('0') && *cp
<= L_('9'))
605 || (TOLOWER (*cp
) >= L_('a') && TOLOWER (*cp
) <= L_('z'))
609 /* The closing brace is missing. Only match the NAN
614 /* This is a system-dependent way to specify the
615 bitmask used for the NaN. We expect it to be
616 a number which is put in the mantissa of the
619 unsigned long long int mant
;
621 mant
= STRTOULL (startp
+ 1, &endp
, 0);
623 SET_MANTISSA (retval
, mant
);
628 *endptr
= (STRING_TYPE
*) cp
;
633 /* It is really a text we do not recognize. */
637 /* First look whether we are faced with a hexadecimal number. */
638 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
640 /* Okay, it is a hexa-decimal number. Remember this and skip
641 the characters. BTW: hexadecimal numbers must not be
649 /* Record the start of the digits, in case we will check their grouping. */
650 start_of_digits
= startp
= cp
;
652 /* Ignore leading zeroes. This helps us to avoid useless computations. */
654 while (c
== L
'0' || ((wint_t) thousands
!= L
'\0' && c
== (wint_t) thousands
))
657 if (thousands
== NULL
)
662 /* We also have the multibyte thousands string. */
667 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
668 if (c
!= thousands
[cnt
])
670 if (thousands
[cnt
] != '\0')
678 /* If no other digit but a '0' is found the result is 0.0.
679 Return current read pointer. */
680 if ((c
< L_('0') || c
> L_('9'))
681 && (base
== 16 && (c
< (CHAR_TYPE
) TOLOWER (L_('a'))
682 || c
> (CHAR_TYPE
) TOLOWER (L_('f'))))
684 && c
!= (wint_t) decimal
686 && ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
687 if (decimal
[cnt
] != cp
[cnt
])
689 decimal
[cnt
] != '\0'; })
691 && (base
== 16 && (cp
== start_of_digits
692 || (CHAR_TYPE
) TOLOWER (c
) != L_('p')))
693 && (base
!= 16 && (CHAR_TYPE
) TOLOWER (c
) != L_('e')))
695 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
696 /* If TP is at the start of the digits, there was no correctly
697 grouped prefix of the string; so no number found. */
698 RETURN (0.0, tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
701 /* Remember first significant digit and read following characters until the
702 decimal point, exponent character or any non-FP number character. */
707 if ((c
>= L_('0') && c
<= L_('9'))
708 || (base
== 16 && (wint_t) TOLOWER (c
) >= L_('a')
709 && (wint_t) TOLOWER (c
) <= L_('f')))
714 if ((wint_t) thousands
== L
'\0' || c
!= (wint_t) thousands
)
715 /* Not a digit or separator: end of the integer part. */
718 if (thousands
== NULL
)
722 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
723 if (thousands
[cnt
] != cp
[cnt
])
725 if (thousands
[cnt
] != '\0')
733 if (grouping
&& dig_no
> 0)
735 /* Check the grouping of the digits. */
736 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
739 /* Less than the entire string was correctly grouped. */
741 if (tp
== start_of_digits
)
742 /* No valid group of numbers at all: no valid number. */
746 /* The number is validly grouped, but consists
747 only of zeroes. The whole value is zero. */
750 /* Recompute DIG_NO so we won't read more digits than
751 are properly grouped. */
754 for (tp
= startp
; tp
< cp
; ++tp
)
755 if (*tp
>= L_('0') && *tp
<= L_('9'))
765 /* We have the number digits in the integer part. Whether these are all or
766 any is really a fractional digit will be decided later. */
768 lead_zero
= int_no
== 0 ? -1 : 0;
770 /* Read the fractional digits. A special case are the 'american style'
771 numbers like `16.' i.e. with decimal but without trailing digits. */
774 c
== (wint_t) decimal
776 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
777 if (decimal
[cnt
] != cp
[cnt
])
779 decimal
[cnt
] == '\0'; })
785 while ((c
>= L_('0') && c
<= L_('9')) ||
786 (base
== 16 && TOLOWER (c
) >= L_('a') && TOLOWER (c
) <= L_('f')))
788 if (c
!= L_('0') && lead_zero
== -1)
789 lead_zero
= dig_no
- int_no
;
795 /* Remember start of exponent (if any). */
799 if ((base
== 16 && TOLOWER (c
) == L_('p'))
800 || (base
!= 16 && TOLOWER (c
) == L_('e')))
802 int exp_negative
= 0;
810 else if (c
== L_('+'))
813 if (c
>= L_('0') && c
<= L_('9'))
817 /* Get the exponent limit. */
819 exp_limit
= (exp_negative
?
820 -MIN_EXP
+ MANT_DIG
+ 4 * int_no
:
821 MAX_EXP
- 4 * int_no
+ lead_zero
);
823 exp_limit
= (exp_negative
?
824 -MIN_10_EXP
+ MANT_DIG
+ int_no
:
825 MAX_10_EXP
- int_no
+ lead_zero
);
831 if (exponent
> exp_limit
)
832 /* The exponent is too large/small to represent a valid
837 /* We have to take care for special situation: a joker
838 might have written "0.0e100000" which is in fact
841 result
= negative
? -0.0 : 0.0;
844 /* Overflow or underflow. */
845 __set_errno (ERANGE
);
846 result
= (exp_negative
? 0.0 :
847 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
850 /* Accept all following digits as part of the exponent. */
853 while (*cp
>= L_('0') && *cp
<= L_('9'));
859 exponent
+= c
- L_('0');
862 while (c
>= L_('0') && c
<= L_('9'));
865 exponent
= -exponent
;
871 /* We don't want to have to work with trailing zeroes after the radix. */
874 while (expp
[-1] == L_('0'))
879 assert (dig_no
>= int_no
);
882 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
885 while (! (base
== 16 ? ISXDIGIT (expp
[-1]) : ISDIGIT (expp
[-1])))
888 if (expp
[-1] != L_('0'))
896 while (dig_no
> 0 && exponent
< 0);
900 /* The whole string is parsed. Store the address of the next character. */
902 *endptr
= (STRING_TYPE
*) cp
;
905 return negative
? -0.0 : 0.0;
909 /* Find the decimal point */
911 while (*startp
!= decimal
)
916 if (*startp
== decimal
[0])
918 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
919 if (decimal
[cnt
] != startp
[cnt
])
921 if (decimal
[cnt
] == '\0')
927 startp
+= lead_zero
+ decimal_len
;
928 exponent
-= base
== 16 ? 4 * lead_zero
: lead_zero
;
932 /* If the BASE is 16 we can use a simpler algorithm. */
935 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
936 4, 4, 4, 4, 4, 4, 4, 4 };
937 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
938 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
941 while (!ISXDIGIT (*startp
))
943 while (*startp
== L_('0'))
945 if (ISDIGIT (*startp
))
946 val
= *startp
++ - L_('0');
948 val
= 10 + TOLOWER (*startp
++) - L_('a');
950 /* We cannot have a leading zero. */
953 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
955 /* We don't have to care for wrapping. This is the normal
956 case so we add the first clause in the `if' expression as
957 an optimization. It is a compile-time constant and so does
958 not cost anything. */
959 retval
[idx
] = val
<< (pos
- bits
+ 1);
964 retval
[idx
--] = val
>> (bits
- pos
- 1);
965 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
966 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
969 /* Adjust the exponent for the bits we are shifting in. */
970 exponent
+= bits
- 1 + (int_no
- 1) * 4;
972 while (--dig_no
> 0 && idx
>= 0)
974 if (!ISXDIGIT (*startp
))
975 startp
+= decimal_len
;
976 if (ISDIGIT (*startp
))
977 val
= *startp
++ - L_('0');
979 val
= 10 + TOLOWER (*startp
++) - L_('a');
983 retval
[idx
] |= val
<< (pos
- 4 + 1);
988 retval
[idx
--] |= val
>> (4 - pos
- 1);
989 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
991 return round_and_return (retval
, exponent
, negative
, val
,
992 BITS_PER_MP_LIMB
- 1, dig_no
> 0);
995 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
999 /* We ran out of digits. */
1000 MPN_ZERO (retval
, idx
);
1002 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
1005 /* Now we have the number of digits in total and the integer digits as well
1006 as the exponent and its sign. We can decide whether the read digits are
1007 really integer digits or belong to the fractional part; i.e. we normalize
1010 register int incr
= (exponent
< 0 ? MAX (-int_no
, exponent
)
1011 : MIN (dig_no
- int_no
, exponent
));
1016 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
1018 __set_errno (ERANGE
);
1019 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1022 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
1024 __set_errno (ERANGE
);
1030 /* Read the integer part as a multi-precision number to NUM. */
1031 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1032 #ifndef USE_WIDE_CHAR
1033 , decimal
, decimal_len
, thousands
1039 /* We now multiply the gained number by the given power of ten. */
1040 mp_limb_t
*psrc
= num
;
1041 mp_limb_t
*pdest
= den
;
1043 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1047 if ((exponent
& expbit
) != 0)
1049 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1053 /* FIXME: not the whole multiplication has to be
1054 done. If we have the needed number of bits we
1055 only need the information whether more non-zero
1057 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1058 cy
= __mpn_mul (pdest
, psrc
, numsize
,
1059 &__tens
[ttab
->arrayoff
1060 + _FPIO_CONST_OFFSET
],
1063 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1064 + _FPIO_CONST_OFFSET
],
1065 size
, psrc
, numsize
);
1069 (void) SWAP (psrc
, pdest
);
1074 while (exponent
!= 0);
1077 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1080 /* Determine how many bits of the result we already have. */
1081 count_leading_zeros (bits
, num
[numsize
- 1]);
1082 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1084 /* Now we know the exponent of the number in base two.
1085 Check it against the maximum possible exponent. */
1088 __set_errno (ERANGE
);
1089 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1092 /* We have already the first BITS bits of the result. Together with
1093 the information whether more non-zero bits follow this is enough
1094 to determine the result. */
1095 if (bits
> MANT_DIG
)
1098 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1099 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1100 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1102 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1106 memcpy (retval
, &num
[least_idx
],
1107 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1110 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1111 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1113 << (BITS_PER_MP_LIMB
- least_bit
));
1114 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1115 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1118 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1119 for (i
= 0; num
[i
] == 0; ++i
)
1122 return round_and_return (retval
, bits
- 1, negative
,
1123 num
[round_idx
], round_bit
,
1124 int_no
< dig_no
|| i
< round_idx
);
1127 else if (dig_no
== int_no
)
1129 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1130 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1132 if (target_bit
== is_bit
)
1134 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1135 numsize
* sizeof (mp_limb_t
));
1136 /* FIXME: the following loop can be avoided if we assume a
1137 maximal MANT_DIG value. */
1138 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1140 else if (target_bit
> is_bit
)
1142 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1143 num
, numsize
, target_bit
- is_bit
);
1144 /* FIXME: the following loop can be avoided if we assume a
1145 maximal MANT_DIG value. */
1146 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1151 assert (numsize
< RETURN_LIMB_SIZE
);
1153 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1154 num
, numsize
, is_bit
- target_bit
);
1155 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1156 /* FIXME: the following loop can be avoided if we assume a
1157 maximal MANT_DIG value. */
1158 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1161 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1165 /* Store the bits we already have. */
1166 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1167 #if RETURN_LIMB_SIZE > 1
1168 if (numsize
< RETURN_LIMB_SIZE
)
1169 retval
[numsize
] = 0;
1173 /* We have to compute at least some of the fractional digits. */
1175 /* We construct a fraction and the result of the division gives us
1176 the needed digits. The denominator is 1.0 multiplied by the
1177 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1178 123e-6 gives 123 / 1000000. */
1184 mp_limb_t
*psrc
= den
;
1185 mp_limb_t
*pdest
= num
;
1186 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1188 assert (dig_no
> int_no
&& exponent
<= 0);
1191 /* For the fractional part we need not process too many digits. One
1192 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1194 digits we should have enough bits for the result. The remaining
1195 decimal digits give us the information that more bits are following.
1196 This can be used while rounding. (Two added as a safety margin.) */
1197 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 2)
1199 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 2;
1205 neg_exp
= dig_no
- int_no
- exponent
;
1207 /* Construct the denominator. */
1212 if ((neg_exp
& expbit
) != 0)
1219 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1220 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1221 densize
* sizeof (mp_limb_t
));
1225 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1226 + _FPIO_CONST_OFFSET
],
1227 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1229 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1232 (void) SWAP (psrc
, pdest
);
1238 while (neg_exp
!= 0);
1241 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1243 /* Read the fractional digits from the string. */
1244 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1245 #ifndef USE_WIDE_CHAR
1246 , decimal
, decimal_len
, thousands
1250 /* We now have to shift both numbers so that the highest bit in the
1251 denominator is set. In the same process we copy the numerator to
1252 a high place in the array so that the division constructs the wanted
1253 digits. This is done by a "quasi fix point" number representation.
1255 num: ddddddddddd . 0000000000000000000000
1257 den: ddddddddddd n >= m
1261 count_leading_zeros (cnt
, den
[densize
- 1]);
1265 /* Don't call `mpn_shift' with a count of zero since the specification
1266 does not allow this. */
1267 (void) __mpn_lshift (den
, den
, densize
, cnt
);
1268 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
1270 num
[numsize
++] = cy
;
1273 /* Now we are ready for the division. But it is not necessary to
1274 do a full multi-precision division because we only need a small
1275 number of bits for the result. So we do not use __mpn_divmod
1276 here but instead do the division here by hand and stop whenever
1277 the needed number of bits is reached. The code itself comes
1278 from the GNU MP Library by Torbj\"orn Granlund. */
1286 mp_limb_t d
, n
, quot
;
1291 assert (numsize
== 1 && n
< d
);
1295 udiv_qrnnd (quot
, n
, n
, 0, d
);
1302 cnt = BITS_PER_MP_LIMB; \
1304 count_leading_zeros (cnt, quot); \
1306 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1308 used = MANT_DIG + cnt; \
1309 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1310 bits = MANT_DIG + 1; \
1314 /* Note that we only clear the second element. */ \
1315 /* The conditional is determined at compile time. */ \
1316 if (RETURN_LIMB_SIZE > 1) \
1322 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1323 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1327 used = MANT_DIG - bits; \
1329 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1331 bits += BITS_PER_MP_LIMB
1335 while (bits
<= MANT_DIG
);
1337 return round_and_return (retval
, exponent
- 1, negative
,
1338 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1339 more_bits
|| n
!= 0);
1343 mp_limb_t d0
, d1
, n0
, n1
;
1350 if (numsize
< densize
)
1354 /* The numerator of the number occupies fewer bits than
1355 the denominator but the one limb is bigger than the
1356 high limb of the numerator. */
1363 exponent
-= BITS_PER_MP_LIMB
;
1366 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1367 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1368 BITS_PER_MP_LIMB
, 0);
1371 used
= MANT_DIG
- bits
;
1373 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1375 bits
+= BITS_PER_MP_LIMB
;
1387 while (bits
<= MANT_DIG
)
1393 /* QUOT should be either 111..111 or 111..110. We need
1394 special treatment of this rare case as normal division
1395 would give overflow. */
1396 quot
= ~(mp_limb_t
) 0;
1399 if (r
< d1
) /* Carry in the addition? */
1401 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1404 n1
= d0
- (d0
!= 0);
1409 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1410 umul_ppmm (n1
, n0
, d0
, quot
);
1414 if (n1
> r
|| (n1
== r
&& n0
> 0))
1416 /* The estimated QUOT was too large. */
1419 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1421 if (r
>= d1
) /* If not carry, test QUOT again. */
1424 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1430 return round_and_return (retval
, exponent
- 1, negative
,
1431 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1432 more_bits
|| n1
!= 0 || n0
!= 0);
1437 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1441 dX
= den
[densize
- 1];
1442 d1
= den
[densize
- 2];
1444 /* The division does not work if the upper limb of the two-limb
1445 numerator is greater than the denominator. */
1446 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1449 if (numsize
< densize
)
1451 mp_size_t empty
= densize
- numsize
;
1456 for (i
= numsize
; i
> 0; --i
)
1457 num
[i
+ empty
] = num
[i
- 1];
1458 MPN_ZERO (num
, empty
+ 1);
1459 exponent
-= empty
* BITS_PER_MP_LIMB
;
1463 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1465 /* We make a difference here because the compiler
1466 cannot optimize the `else' case that good and
1467 this reflects all currently used FLOAT types
1468 and GMP implementations. */
1470 #if RETURN_LIMB_SIZE <= 2
1471 assert (empty
== 1);
1472 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1473 BITS_PER_MP_LIMB
, 0);
1475 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1476 retval
[i
] = retval
[i
- empty
];
1478 for (i
= numsize
; i
> 0; --i
)
1479 num
[i
+ empty
] = num
[i
- 1];
1480 MPN_ZERO (num
, empty
+ 1);
1484 used
= MANT_DIG
- bits
;
1485 if (used
>= BITS_PER_MP_LIMB
)
1488 (void) __mpn_lshift (&retval
[used
1489 / BITS_PER_MP_LIMB
],
1490 retval
, RETURN_LIMB_SIZE
,
1491 used
% BITS_PER_MP_LIMB
);
1492 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1496 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1498 bits
+= empty
* BITS_PER_MP_LIMB
;
1504 assert (numsize
== densize
);
1505 for (i
= numsize
; i
> 0; --i
)
1506 num
[i
] = num
[i
- 1];
1512 while (bits
<= MANT_DIG
)
1515 /* This might over-estimate QUOT, but it's probably not
1516 worth the extra code here to find out. */
1517 quot
= ~(mp_limb_t
) 0;
1522 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1523 umul_ppmm (n1
, n0
, d1
, quot
);
1525 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1529 if (r
< dX
) /* I.e. "carry in previous addition?" */
1536 /* Possible optimization: We already have (q * n0) and (1 * n1)
1537 after the calculation of QUOT. Taking advantage of this, we
1538 could make this loop make two iterations less. */
1540 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1542 if (num
[densize
] != cy
)
1544 cy
= __mpn_add_n (num
, num
, den
, densize
);
1548 n0
= num
[densize
] = num
[densize
- 1];
1549 for (i
= densize
- 1; i
> 0; --i
)
1550 num
[i
] = num
[i
- 1];
1555 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1557 return round_and_return (retval
, exponent
- 1, negative
,
1558 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1559 more_bits
|| i
>= 0);
1567 && !(defined USE_IN_EXTENDED_LOCALE_MODEL && defined USE_WIDE_CHAR)
1568 libc_hidden_def (INTERNAL (STRTOF
))
1571 /* External user entry point. */
1574 #ifdef weak_function
1577 STRTOF (nptr
, endptr LOCALE_PARAM
)
1578 const STRING_TYPE
*nptr
;
1579 STRING_TYPE
**endptr
;
1582 return INTERNAL (STRTOF
) (nptr
, endptr
, 0 LOCALE_PARAM
);