1 /* Read decimal floating point numbers.
2 This file is part of the GNU C Library.
3 Copyright (C) 1995,96,97,98,99,2000,01,02 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. */
69 #include <gmp-mparam.h>
71 #include "fpioconst.h"
77 /* We use this code also for the extended locale handling where the
78 function gets as an additional argument the locale which has to be
79 used. To access the values we have to redefine the _NL_CURRENT
81 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
83 # define _NL_CURRENT(category, item) \
84 (current->values[_NL_ITEM_INDEX (item)].string)
85 # define LOCALE_PARAM , loc
86 # define LOCALE_PARAM_DECL __locale_t loc;
89 # define LOCALE_PARAM_DECL
92 #if defined _LIBC || defined HAVE_WCHAR_H
98 # define STRING_TYPE wchar_t
99 # define CHAR_TYPE wint_t
100 # define L_(Ch) L##Ch
101 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
102 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
103 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
104 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
105 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
106 # define STRNCASECMP(S1, S2, N) __wcsncasecmp_l ((S1), (S2), (N), loc)
107 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
109 # define ISSPACE(Ch) iswspace (Ch)
110 # define ISDIGIT(Ch) iswdigit (Ch)
111 # define ISXDIGIT(Ch) iswxdigit (Ch)
112 # define TOLOWER(Ch) towlower (Ch)
113 # define STRNCASECMP(S1, S2, N) __wcsncasecmp ((S1), (S2), (N))
114 # define STRTOULL(S, E, B) __wcstoull_internal ((S), (E), (B), 0)
117 # define STRING_TYPE char
118 # define CHAR_TYPE char
120 # ifdef USE_IN_EXTENDED_LOCALE_MODEL
121 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
122 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
123 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
124 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
125 # define STRNCASECMP(S1, S2, N) __strncasecmp_l ((S1), (S2), (N), loc)
126 # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc)
128 # define ISSPACE(Ch) isspace (Ch)
129 # define ISDIGIT(Ch) isdigit (Ch)
130 # define ISXDIGIT(Ch) isxdigit (Ch)
131 # define TOLOWER(Ch) tolower (Ch)
132 # define STRNCASECMP(S1, S2, N) __strncasecmp ((S1), (S2), (N))
133 # define STRTOULL(S, E, B) __strtoull_internal ((S), (E), 0, (B))
138 /* Constants we need from float.h; select the set for the FLOAT precision. */
139 #define MANT_DIG PASTE(FLT,_MANT_DIG)
140 #define DIG PASTE(FLT,_DIG)
141 #define MAX_EXP PASTE(FLT,_MAX_EXP)
142 #define MIN_EXP PASTE(FLT,_MIN_EXP)
143 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
144 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
146 /* Extra macros required to get FLT expanded before the pasting. */
147 #define PASTE(a,b) PASTE1(a,b)
148 #define PASTE1(a,b) a##b
150 /* Function to construct a floating point number from an MP integer
151 containing the fraction bits, a base 2 exponent, and a sign flag. */
152 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
154 /* Definitions according to limb size used. */
155 #if BITS_PER_MP_LIMB == 32
156 # define MAX_DIG_PER_LIMB 9
157 # define MAX_FAC_PER_LIMB 1000000000UL
158 #elif BITS_PER_MP_LIMB == 64
159 # define MAX_DIG_PER_LIMB 19
160 # define MAX_FAC_PER_LIMB 10000000000000000000UL
162 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
166 /* Local data structure. */
167 static const mp_limb_t _tens_in_limb
[MAX_DIG_PER_LIMB
+ 1] =
170 1000000, 10000000, 100000000,
172 #if BITS_PER_MP_LIMB > 32
173 , 10000000000U, 100000000000U,
174 1000000000000U, 10000000000000U, 100000000000000U,
175 1000000000000000U, 10000000000000000U, 100000000000000000U,
176 1000000000000000000U, 10000000000000000000U
178 #if BITS_PER_MP_LIMB > 64
179 #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
184 #define howmany(x,y) (((x)+((y)-1))/(y))
186 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
188 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
189 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
190 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
192 #define RETURN(val,end) \
193 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
194 return val; } while (0)
196 /* Maximum size necessary for mpn integers to hold floating point numbers. */
197 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
199 /* Declare an mpn integer variable that big. */
200 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
201 /* Copy an mpn integer value. */
202 #define MPN_ASSIGN(dst, src) \
203 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
206 /* Return a floating point number of the needed type according to the given
207 multi-precision number after possible rounding. */
209 round_and_return (mp_limb_t
*retval
, int exponent
, int negative
,
210 mp_limb_t round_limb
, mp_size_t round_bit
, int more_bits
)
212 if (exponent
< MIN_EXP
- 1)
214 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
216 if (shift
> MANT_DIG
)
222 more_bits
|= (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0;
223 if (shift
== MANT_DIG
)
224 /* This is a special case to handle the very seldom case where
225 the mantissa will be empty after the shift. */
229 round_limb
= retval
[RETURN_LIMB_SIZE
- 1];
230 round_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
231 for (i
= 0; i
< RETURN_LIMB_SIZE
; ++i
)
232 more_bits
|= retval
[i
] != 0;
233 MPN_ZERO (retval
, RETURN_LIMB_SIZE
);
235 else if (shift
>= BITS_PER_MP_LIMB
)
239 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
240 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
241 for (i
= 0; i
< (shift
- 1) / BITS_PER_MP_LIMB
; ++i
)
242 more_bits
|= retval
[i
] != 0;
243 more_bits
|= ((round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1))
246 (void) __mpn_rshift (retval
, &retval
[shift
/ BITS_PER_MP_LIMB
],
247 RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
),
248 shift
% BITS_PER_MP_LIMB
);
249 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
)],
250 shift
/ BITS_PER_MP_LIMB
);
254 round_limb
= retval
[0];
255 round_bit
= shift
- 1;
256 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
258 /* This is a hook for the m68k long double format, where the
259 exponent bias is the same for normalized and denormalized
262 # define DENORM_EXP (MIN_EXP - 2)
264 exponent
= DENORM_EXP
;
267 if ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
268 && (more_bits
|| (retval
[0] & 1) != 0
269 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
271 mp_limb_t cy
= __mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
273 if (((MANT_DIG
% BITS_PER_MP_LIMB
) == 0 && cy
) ||
274 ((MANT_DIG
% BITS_PER_MP_LIMB
) != 0 &&
275 (retval
[RETURN_LIMB_SIZE
- 1]
276 & (((mp_limb_t
) 1) << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0))
279 (void) __mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
280 retval
[RETURN_LIMB_SIZE
- 1]
281 |= ((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
);
283 else if (exponent
== DENORM_EXP
284 && (retval
[RETURN_LIMB_SIZE
- 1]
285 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
287 /* The number was denormalized but now normalized. */
288 exponent
= MIN_EXP
- 1;
291 if (exponent
> MAX_EXP
)
292 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
294 return MPN2FLOAT (retval
, exponent
, negative
);
298 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
299 into N. Return the size of the number limbs in NSIZE at the first
300 character od the string that is not part of the integer as the function
301 value. If the EXPONENT is small enough to be taken as an additional
302 factor for the resulting number (see code) multiply by it. */
303 static inline const STRING_TYPE
*
304 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
306 #ifndef USE_WIDE_CHAR
307 , const char *decimal
, size_t decimal_len
, const char *thousands
312 /* Number of digits for actual limb. */
321 if (cnt
== MAX_DIG_PER_LIMB
)
331 cy
= __mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
332 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
343 /* There might be thousands separators or radix characters in
344 the string. But these all can be ignored because we know the
345 format of the number is correct and we have an exact number
346 of characters to read. */
348 if (*str
< L
'0' || *str
> L
'9')
351 if (*str
< '0' || *str
> '9')
354 if (thousands
!= NULL
&& *str
== *thousands
355 && ({ for (inner
= 1; thousands
[inner
] != '\0'; ++inner
)
356 if (thousands
[inner
] != str
[inner
])
358 thousands
[inner
] == '\0'; }))
364 low
= low
* 10 + *str
++ - L_('0');
367 while (--digcnt
> 0);
369 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
371 low
*= _tens_in_limb
[*exponent
];
372 start
= _tens_in_limb
[cnt
+ *exponent
];
376 start
= _tens_in_limb
[cnt
];
386 cy
= __mpn_mul_1 (n
, n
, *nsize
, start
);
387 cy
+= __mpn_add_1 (n
, n
, *nsize
, low
);
396 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
397 with the COUNT most significant bits of LIMB.
399 Tege doesn't like this function so I have to write it here myself. :)
402 __mpn_lshift_1 (mp_limb_t
*ptr
, mp_size_t size
, unsigned int count
,
405 if (count
== BITS_PER_MP_LIMB
)
407 /* Optimize the case of shifting by exactly a word:
408 just copy words, with no actual bit-shifting. */
410 for (i
= size
- 1; i
> 0; --i
)
416 (void) __mpn_lshift (ptr
, ptr
, size
, count
);
417 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
422 #define INTERNAL(x) INTERNAL1(x)
423 #define INTERNAL1(x) __##x##_internal
425 /* This file defines a function to check for correct grouping. */
426 #include "grouping.h"
429 /* Return a floating point number with the value of the given string NPTR.
430 Set *ENDPTR to the character after the last used one. If the number is
431 smaller than the smallest representable number, set `errno' to ERANGE and
432 return 0.0. If the number is too big to be represented, set `errno' to
433 ERANGE and return HUGE_VAL with the appropriate sign. */
435 INTERNAL (STRTOF
) (nptr
, endptr
, group LOCALE_PARAM
)
436 const STRING_TYPE
*nptr
;
437 STRING_TYPE
**endptr
;
441 int negative
; /* The sign of the number. */
442 MPN_VAR (num
); /* MP representation of the number. */
443 int exponent
; /* Exponent of the number. */
445 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
448 /* When we have to compute fractional digits we form a fraction with a
449 second multi-precision number (and we sometimes need a second for
450 temporary results). */
453 /* Representation for the return value. */
454 mp_limb_t retval
[RETURN_LIMB_SIZE
];
455 /* Number of bits currently in result value. */
458 /* Running pointer after the last character processed in the string. */
459 const STRING_TYPE
*cp
, *tp
;
460 /* Start of significant part of the number. */
461 const STRING_TYPE
*startp
, *start_of_digits
;
462 /* Points at the character following the integer and fractional digits. */
463 const STRING_TYPE
*expp
;
464 /* Total number of digit and number of digits in integer part. */
465 int dig_no
, int_no
, lead_zero
;
466 /* Contains the last character read. */
469 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
470 there. So define it ourselves if it remains undefined. */
472 typedef unsigned int wint_t;
474 /* The radix character of the current locale. */
481 /* The thousands character of the current locale. */
483 wchar_t thousands
= L
'\0';
485 const char *thousands
= NULL
;
487 /* The numeric grouping specification of the current locale,
488 in the format described in <locale.h>. */
489 const char *grouping
;
490 /* Used in several places. */
493 #ifdef USE_IN_EXTENDED_LOCALE_MODEL
494 struct locale_data
*current
= loc
->__locales
[LC_NUMERIC
];
499 grouping
= _NL_CURRENT (LC_NUMERIC
, GROUPING
);
500 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
504 /* Figure out the thousands separator character. */
506 thousands
= _NL_CURRENT_WORD (LC_NUMERIC
,
507 _NL_NUMERIC_THOUSANDS_SEP_WC
);
508 if (thousands
== L
'\0')
511 thousands
= _NL_CURRENT (LC_NUMERIC
, THOUSANDS_SEP
);
512 if (*thousands
== '\0')
523 /* Find the locale's decimal point character. */
525 decimal
= _NL_CURRENT_WORD (LC_NUMERIC
, _NL_NUMERIC_DECIMAL_POINT_WC
);
526 assert (decimal
!= L
'\0');
527 # define decimal_len 1
529 decimal
= _NL_CURRENT (LC_NUMERIC
, DECIMAL_POINT
);
530 decimal_len
= strlen (decimal
);
531 assert (decimal_len
> 0);
534 /* Prepare number representation. */
539 /* Parse string to get maximal legal prefix. We need the number of
540 characters of the integer part, the fractional part and the exponent. */
542 /* Ignore leading white space. */
547 /* Get sign of the result. */
553 else if (c
== L_('+'))
556 /* Return 0.0 if no legal string is found.
557 No character is used even if a sign was found. */
559 if (c
== decimal
&& cp
[1] >= L
'0' && cp
[1] <= L
'9')
561 /* We accept it. This funny construct is here only to indent
562 the code directly. */
565 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
566 if (cp
[cnt
] != decimal
[cnt
])
568 if (decimal
[cnt
] == '\0' && cp
[1] >= '0' && cp
[1] <= '9')
570 /* We accept it. This funny construct is here only to indent
571 the code directly. */
574 else if (c
< L_('0') || c
> L_('9'))
576 /* Check for `INF' or `INFINITY'. */
577 if (TOLOWER (c
) == L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
579 /* Return +/- infinity. */
581 *endptr
= (STRING_TYPE
*)
582 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
585 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
588 if (TOLOWER (c
) == L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
595 /* Match `(n-char-sequence-digit)'. */
598 const STRING_TYPE
*startp
= cp
;
601 while ((*cp
>= L_('0') && *cp
<= L_('9'))
602 || (TOLOWER (*cp
) >= L_('a') && TOLOWER (*cp
) <= L_('z'))
606 /* The closing brace is missing. Only match the NAN
611 /* This is a system-dependent way to specify the
612 bitmask used for the NaN. We expect it to be
613 a number which is put in the mantissa of the
616 unsigned long long int mant
;
618 mant
= STRTOULL (startp
+ 1, &endp
, 0);
620 SET_MANTISSA (retval
, mant
);
625 *endptr
= (STRING_TYPE
*) cp
;
630 /* It is really a text we do not recognize. */
634 /* First look whether we are faced with a hexadecimal number. */
635 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
637 /* Okay, it is a hexa-decimal number. Remember this and skip
638 the characters. BTW: hexadecimal numbers must not be
646 /* Record the start of the digits, in case we will check their grouping. */
647 start_of_digits
= startp
= cp
;
649 /* Ignore leading zeroes. This helps us to avoid useless computations. */
651 while (c
== L
'0' || (thousands
!= L
'\0' && c
== thousands
))
654 if (thousands
== NULL
)
659 /* We also have the multibyte thousands string. */
664 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
665 if (c
!= thousands
[cnt
])
667 if (thousands
[cnt
] != '\0')
675 /* If no other digit but a '0' is found the result is 0.0.
676 Return current read pointer. */
677 if ((c
< L_('0') || c
> L_('9'))
678 && (base
== 16 && (c
< TOLOWER (L_('a')) || c
> TOLOWER (L_('f'))))
682 && ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
683 if (decimal
[cnt
] != cp
[cnt
])
685 decimal
[cnt
] != '\0'; })
687 && (base
== 16 && (cp
== start_of_digits
|| TOLOWER (c
) != L_('p')))
688 && (base
!= 16 && TOLOWER (c
) != L_('e')))
690 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
691 /* If TP is at the start of the digits, there was no correctly
692 grouped prefix of the string; so no number found. */
693 RETURN (0.0, tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
696 /* Remember first significant digit and read following characters until the
697 decimal point, exponent character or any non-FP number character. */
702 if ((c
>= L_('0') && c
<= L_('9'))
703 || (base
== 16 && TOLOWER (c
) >= L_('a') && TOLOWER (c
) <= L_('f')))
708 if (thousands
== L
'\0' || c
!= thousands
)
709 /* Not a digit or separator: end of the integer part. */
712 if (thousands
== NULL
)
716 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
717 if (thousands
[cnt
] != cp
[cnt
])
719 if (thousands
[cnt
] != '\0')
727 if (grouping
&& dig_no
> 0)
729 /* Check the grouping of the digits. */
730 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
733 /* Less than the entire string was correctly grouped. */
735 if (tp
== start_of_digits
)
736 /* No valid group of numbers at all: no valid number. */
740 /* The number is validly grouped, but consists
741 only of zeroes. The whole value is zero. */
744 /* Recompute DIG_NO so we won't read more digits than
745 are properly grouped. */
748 for (tp
= startp
; tp
< cp
; ++tp
)
749 if (*tp
>= L_('0') && *tp
<= L_('9'))
759 /* We have the number digits in the integer part. Whether these are all or
760 any is really a fractional digit will be decided later. */
762 lead_zero
= int_no
== 0 ? -1 : 0;
764 /* Read the fractional digits. A special case are the 'american style'
765 numbers like `16.' i.e. with decimal but without trailing digits. */
770 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
771 if (decimal
[cnt
] != cp
[cnt
])
773 decimal
[cnt
] == '\0'; })
779 while ((c
>= L_('0') && c
<= L_('9')) ||
780 (base
== 16 && TOLOWER (c
) >= L_('a') && TOLOWER (c
) <= L_('f')))
782 if (c
!= L_('0') && lead_zero
== -1)
783 lead_zero
= dig_no
- int_no
;
789 /* Remember start of exponent (if any). */
793 if ((base
== 16 && TOLOWER (c
) == L_('p'))
794 || (base
!= 16 && TOLOWER (c
) == L_('e')))
796 int exp_negative
= 0;
804 else if (c
== L_('+'))
807 if (c
>= L_('0') && c
<= L_('9'))
811 /* Get the exponent limit. */
813 exp_limit
= (exp_negative
?
814 -MIN_EXP
+ MANT_DIG
+ 4 * int_no
:
815 MAX_EXP
- 4 * int_no
+ lead_zero
);
817 exp_limit
= (exp_negative
?
818 -MIN_10_EXP
+ MANT_DIG
+ int_no
:
819 MAX_10_EXP
- int_no
+ lead_zero
);
825 if (exponent
> exp_limit
)
826 /* The exponent is too large/small to represent a valid
831 /* We have to take care for special situation: a joker
832 might have written "0.0e100000" which is in fact
835 result
= negative
? -0.0 : 0.0;
838 /* Overflow or underflow. */
839 __set_errno (ERANGE
);
840 result
= (exp_negative
? 0.0 :
841 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
844 /* Accept all following digits as part of the exponent. */
847 while (*cp
>= L_('0') && *cp
<= L_('9'));
853 exponent
+= c
- L_('0');
856 while (c
>= L_('0') && c
<= L_('9'));
859 exponent
= -exponent
;
865 /* We don't want to have to work with trailing zeroes after the radix. */
868 while (expp
[-1] == L_('0'))
873 assert (dig_no
>= int_no
);
876 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
879 while (expp
[-1] < L_('0') || expp
[-1] > L_('9'))
882 if (expp
[-1] != L_('0'))
890 while (dig_no
> 0 && exponent
< 0);
894 /* The whole string is parsed. Store the address of the next character. */
896 *endptr
= (STRING_TYPE
*) cp
;
899 return negative
? -0.0 : 0.0;
903 /* Find the decimal point */
905 while (*startp
!= decimal
)
910 if (*startp
== decimal
[0])
912 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
913 if (decimal
[cnt
] != startp
[cnt
])
915 if (decimal
[cnt
] == '\0')
921 startp
+= lead_zero
+ decimal_len
;
922 exponent
-= base
== 16 ? 4 * lead_zero
: lead_zero
;
926 /* If the BASE is 16 we can use a simpler algorithm. */
929 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
930 4, 4, 4, 4, 4, 4, 4, 4 };
931 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
932 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
935 while (!ISXDIGIT (*startp
))
937 while (*startp
== L_('0'))
939 if (ISDIGIT (*startp
))
940 val
= *startp
++ - L_('0');
942 val
= 10 + TOLOWER (*startp
++) - L_('a');
944 /* We cannot have a leading zero. */
947 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
949 /* We don't have to care for wrapping. This is the normal
950 case so we add the first clause in the `if' expression as
951 an optimization. It is a compile-time constant and so does
952 not cost anything. */
953 retval
[idx
] = val
<< (pos
- bits
+ 1);
958 retval
[idx
--] = val
>> (bits
- pos
- 1);
959 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
960 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
963 /* Adjust the exponent for the bits we are shifting in. */
964 exponent
+= bits
- 1 + (int_no
- 1) * 4;
966 while (--dig_no
> 0 && idx
>= 0)
968 if (!ISXDIGIT (*startp
))
969 startp
+= decimal_len
;
970 if (ISDIGIT (*startp
))
971 val
= *startp
++ - L_('0');
973 val
= 10 + TOLOWER (*startp
++) - L_('a');
977 retval
[idx
] |= val
<< (pos
- 4 + 1);
982 retval
[idx
--] |= val
>> (4 - pos
- 1);
983 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
985 return round_and_return (retval
, exponent
, negative
, val
,
986 BITS_PER_MP_LIMB
- 1, dig_no
> 0);
989 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
993 /* We ran out of digits. */
994 MPN_ZERO (retval
, idx
);
996 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
999 /* Now we have the number of digits in total and the integer digits as well
1000 as the exponent and its sign. We can decide whether the read digits are
1001 really integer digits or belong to the fractional part; i.e. we normalize
1004 register int incr
= (exponent
< 0 ? MAX (-int_no
, exponent
)
1005 : MIN (dig_no
- int_no
, exponent
));
1010 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
1012 __set_errno (ERANGE
);
1013 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1016 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
1018 __set_errno (ERANGE
);
1024 /* Read the integer part as a multi-precision number to NUM. */
1025 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1026 #ifndef USE_WIDE_CHAR
1027 , decimal
, decimal_len
, thousands
1033 /* We now multiply the gained number by the given power of ten. */
1034 mp_limb_t
*psrc
= num
;
1035 mp_limb_t
*pdest
= den
;
1037 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1041 if ((exponent
& expbit
) != 0)
1043 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1047 /* FIXME: not the whole multiplication has to be
1048 done. If we have the needed number of bits we
1049 only need the information whether more non-zero
1051 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1052 cy
= __mpn_mul (pdest
, psrc
, numsize
,
1053 &__tens
[ttab
->arrayoff
1054 + _FPIO_CONST_OFFSET
],
1057 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1058 + _FPIO_CONST_OFFSET
],
1059 size
, psrc
, numsize
);
1063 (void) SWAP (psrc
, pdest
);
1068 while (exponent
!= 0);
1071 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1074 /* Determine how many bits of the result we already have. */
1075 count_leading_zeros (bits
, num
[numsize
- 1]);
1076 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1078 /* Now we know the exponent of the number in base two.
1079 Check it against the maximum possible exponent. */
1082 __set_errno (ERANGE
);
1083 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1086 /* We have already the first BITS bits of the result. Together with
1087 the information whether more non-zero bits follow this is enough
1088 to determine the result. */
1089 if (bits
> MANT_DIG
)
1092 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1093 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1094 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1096 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1100 memcpy (retval
, &num
[least_idx
],
1101 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1104 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1105 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1107 << (BITS_PER_MP_LIMB
- least_bit
));
1108 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1109 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1112 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1113 for (i
= 0; num
[i
] == 0; ++i
)
1116 return round_and_return (retval
, bits
- 1, negative
,
1117 num
[round_idx
], round_bit
,
1118 int_no
< dig_no
|| i
< round_idx
);
1121 else if (dig_no
== int_no
)
1123 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1124 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1126 if (target_bit
== is_bit
)
1128 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1129 numsize
* sizeof (mp_limb_t
));
1130 /* FIXME: the following loop can be avoided if we assume a
1131 maximal MANT_DIG value. */
1132 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1134 else if (target_bit
> is_bit
)
1136 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1137 num
, numsize
, target_bit
- is_bit
);
1138 /* FIXME: the following loop can be avoided if we assume a
1139 maximal MANT_DIG value. */
1140 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1145 assert (numsize
< RETURN_LIMB_SIZE
);
1147 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1148 num
, numsize
, is_bit
- target_bit
);
1149 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1150 /* FIXME: the following loop can be avoided if we assume a
1151 maximal MANT_DIG value. */
1152 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1155 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1159 /* Store the bits we already have. */
1160 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1161 #if RETURN_LIMB_SIZE > 1
1162 if (numsize
< RETURN_LIMB_SIZE
)
1163 retval
[numsize
] = 0;
1167 /* We have to compute at least some of the fractional digits. */
1169 /* We construct a fraction and the result of the division gives us
1170 the needed digits. The denominator is 1.0 multiplied by the
1171 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1172 123e-6 gives 123 / 1000000. */
1178 mp_limb_t
*psrc
= den
;
1179 mp_limb_t
*pdest
= num
;
1180 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1182 assert (dig_no
> int_no
&& exponent
<= 0);
1185 /* For the fractional part we need not process too many digits. One
1186 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1188 digits we should have enough bits for the result. The remaining
1189 decimal digits give us the information that more bits are following.
1190 This can be used while rounding. (One added as a safety margin.) */
1191 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
1193 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
1199 neg_exp
= dig_no
- int_no
- exponent
;
1201 /* Construct the denominator. */
1206 if ((neg_exp
& expbit
) != 0)
1213 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1214 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1215 densize
* sizeof (mp_limb_t
));
1219 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1220 + _FPIO_CONST_OFFSET
],
1221 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1223 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1226 (void) SWAP (psrc
, pdest
);
1232 while (neg_exp
!= 0);
1235 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1237 /* Read the fractional digits from the string. */
1238 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1239 #ifndef USE_WIDE_CHAR
1240 , decimal
, decimal_len
, thousands
1244 /* We now have to shift both numbers so that the highest bit in the
1245 denominator is set. In the same process we copy the numerator to
1246 a high place in the array so that the division constructs the wanted
1247 digits. This is done by a "quasi fix point" number representation.
1249 num: ddddddddddd . 0000000000000000000000
1251 den: ddddddddddd n >= m
1255 count_leading_zeros (cnt
, den
[densize
- 1]);
1259 /* Don't call `mpn_shift' with a count of zero since the specification
1260 does not allow this. */
1261 (void) __mpn_lshift (den
, den
, densize
, cnt
);
1262 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
1264 num
[numsize
++] = cy
;
1267 /* Now we are ready for the division. But it is not necessary to
1268 do a full multi-precision division because we only need a small
1269 number of bits for the result. So we do not use __mpn_divmod
1270 here but instead do the division here by hand and stop whenever
1271 the needed number of bits is reached. The code itself comes
1272 from the GNU MP Library by Torbj\"orn Granlund. */
1280 mp_limb_t d
, n
, quot
;
1285 assert (numsize
== 1 && n
< d
);
1289 udiv_qrnnd (quot
, n
, n
, 0, d
);
1296 cnt = BITS_PER_MP_LIMB; \
1298 count_leading_zeros (cnt, quot); \
1300 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1302 used = MANT_DIG + cnt; \
1303 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1304 bits = MANT_DIG + 1; \
1308 /* Note that we only clear the second element. */ \
1309 /* The conditional is determined at compile time. */ \
1310 if (RETURN_LIMB_SIZE > 1) \
1316 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1317 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1321 used = MANT_DIG - bits; \
1323 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1325 bits += BITS_PER_MP_LIMB
1329 while (bits
<= MANT_DIG
);
1331 return round_and_return (retval
, exponent
- 1, negative
,
1332 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1333 more_bits
|| n
!= 0);
1337 mp_limb_t d0
, d1
, n0
, n1
;
1344 if (numsize
< densize
)
1348 /* The numerator of the number occupies fewer bits than
1349 the denominator but the one limb is bigger than the
1350 high limb of the numerator. */
1357 exponent
-= BITS_PER_MP_LIMB
;
1360 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1361 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1362 BITS_PER_MP_LIMB
, 0);
1365 used
= MANT_DIG
- bits
;
1367 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1369 bits
+= BITS_PER_MP_LIMB
;
1381 while (bits
<= MANT_DIG
)
1387 /* QUOT should be either 111..111 or 111..110. We need
1388 special treatment of this rare case as normal division
1389 would give overflow. */
1390 quot
= ~(mp_limb_t
) 0;
1393 if (r
< d1
) /* Carry in the addition? */
1395 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1398 n1
= d0
- (d0
!= 0);
1403 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1404 umul_ppmm (n1
, n0
, d0
, quot
);
1408 if (n1
> r
|| (n1
== r
&& n0
> 0))
1410 /* The estimated QUOT was too large. */
1413 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1415 if (r
>= d1
) /* If not carry, test QUOT again. */
1418 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1424 return round_and_return (retval
, exponent
- 1, negative
,
1425 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1426 more_bits
|| n1
!= 0 || n0
!= 0);
1431 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1435 dX
= den
[densize
- 1];
1436 d1
= den
[densize
- 2];
1438 /* The division does not work if the upper limb of the two-limb
1439 numerator is greater than the denominator. */
1440 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1443 if (numsize
< densize
)
1445 mp_size_t empty
= densize
- numsize
;
1450 for (i
= numsize
; i
> 0; --i
)
1451 num
[i
+ empty
] = num
[i
- 1];
1452 MPN_ZERO (num
, empty
+ 1);
1453 exponent
-= empty
* BITS_PER_MP_LIMB
;
1457 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1459 /* We make a difference here because the compiler
1460 cannot optimize the `else' case that good and
1461 this reflects all currently used FLOAT types
1462 and GMP implementations. */
1464 #if RETURN_LIMB_SIZE <= 2
1465 assert (empty
== 1);
1466 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1467 BITS_PER_MP_LIMB
, 0);
1469 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1470 retval
[i
] = retval
[i
- empty
];
1472 #if RETURN_LIMB_SIZE > 1
1475 for (i
= numsize
; i
> 0; --i
)
1476 num
[i
+ empty
] = num
[i
- 1];
1477 MPN_ZERO (num
, empty
+ 1);
1481 used
= MANT_DIG
- bits
;
1482 if (used
>= BITS_PER_MP_LIMB
)
1485 (void) __mpn_lshift (&retval
[used
1486 / BITS_PER_MP_LIMB
],
1487 retval
, RETURN_LIMB_SIZE
,
1488 used
% BITS_PER_MP_LIMB
);
1489 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1493 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1495 bits
+= empty
* BITS_PER_MP_LIMB
;
1501 assert (numsize
== densize
);
1502 for (i
= numsize
; i
> 0; --i
)
1503 num
[i
] = num
[i
- 1];
1509 while (bits
<= MANT_DIG
)
1512 /* This might over-estimate QUOT, but it's probably not
1513 worth the extra code here to find out. */
1514 quot
= ~(mp_limb_t
) 0;
1519 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1520 umul_ppmm (n1
, n0
, d1
, quot
);
1522 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1526 if (r
< dX
) /* I.e. "carry in previous addition?" */
1533 /* Possible optimization: We already have (q * n0) and (1 * n1)
1534 after the calculation of QUOT. Taking advantage of this, we
1535 could make this loop make two iterations less. */
1537 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1539 if (num
[densize
] != cy
)
1541 cy
= __mpn_add_n (num
, num
, den
, densize
);
1545 n0
= num
[densize
] = num
[densize
- 1];
1546 for (i
= densize
- 1; i
> 0; --i
)
1547 num
[i
] = num
[i
- 1];
1552 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1554 return round_and_return (retval
, exponent
- 1, negative
,
1555 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1556 more_bits
|| i
>= 0);
1564 && !(defined USE_IN_EXTENDED_LOCALE_MODEL && defined USE_WIDE_CHAR)
1565 libc_hidden_def (INTERNAL (STRTOF
))
1568 /* External user entry point. */
1571 #ifdef weak_function
1574 STRTOF (nptr
, endptr LOCALE_PARAM
)
1575 const STRING_TYPE
*nptr
;
1576 STRING_TYPE
**endptr
;
1579 return INTERNAL (STRTOF
) (nptr
, endptr
, 0 LOCALE_PARAM
);