1 /* Convert string representing a number to float value, using given locale.
2 Copyright (C) 1997,1998,2002,2004,2005,2006,2007,2008,2009,2010
3 Free Software Foundation, Inc.
4 This file is part of the GNU C Library.
5 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
7 The GNU C Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU Lesser General Public
9 License as published by the Free Software Foundation; either
10 version 2.1 of the License, or (at your option) any later version.
12 The GNU C Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 Lesser General Public License for more details.
17 You should have received a copy of the GNU Lesser General Public
18 License along with the GNU C Library; if not, write to the Free
19 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
32 #include "../printf/quadmath-printf.h"
33 #include "../printf/fpioconst.h"
38 # define STRING_TYPE wchar_t
39 # define CHAR_TYPE wint_t
41 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
42 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
43 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
44 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
45 # define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
46 # define STRNCASECMP(S1, S2, N) \
47 __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
48 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
50 # define STRING_TYPE char
51 # define CHAR_TYPE char
53 # define ISSPACE(Ch) isspace (Ch)
54 # define ISDIGIT(Ch) isdigit (Ch)
55 # define ISXDIGIT(Ch) isxdigit (Ch)
56 # define TOLOWER(Ch) tolower (Ch)
57 # define TOLOWER_C(Ch) \
58 ({__typeof(Ch) __tlc = (Ch); \
59 (__tlc >= 'A' && __tlc <= 'Z') ? __tlc - 'A' + 'a' : __tlc; })
60 # define STRNCASECMP(S1, S2, N) \
61 __quadmath_strncasecmp_c (S1, S2, N)
63 # define STRTOULL(S, E, B) strtoull (S, E, B)
65 # define STRTOULL(S, E, B) strtoul (S, E, B)
69 __quadmath_strncasecmp_c (const char *s1
, const char *s2
, size_t n
)
71 const unsigned char *p1
= (const unsigned char *) s1
;
72 const unsigned char *p2
= (const unsigned char *) s2
;
74 if (p1
== p2
|| n
== 0)
76 while ((result
= TOLOWER_C (*p1
) - TOLOWER_C (*p2
++)) == 0)
77 if (*p1
++ == '\0' || --n
== 0)
85 /* Constants we need from float.h; select the set for the FLOAT precision. */
86 #define MANT_DIG PASTE(FLT,_MANT_DIG)
87 #define DIG PASTE(FLT,_DIG)
88 #define MAX_EXP PASTE(FLT,_MAX_EXP)
89 #define MIN_EXP PASTE(FLT,_MIN_EXP)
90 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
91 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
93 /* Extra macros required to get FLT expanded before the pasting. */
94 #define PASTE(a,b) PASTE1(a,b)
95 #define PASTE1(a,b) a##b
97 /* Function to construct a floating point number from an MP integer
98 containing the fraction bits, a base 2 exponent, and a sign flag. */
99 extern FLOAT
MPN2FLOAT (mp_srcptr mpn
, int exponent
, int negative
);
101 /* Definitions according to limb size used. */
102 #if BITS_PER_MP_LIMB == 32
103 # define MAX_DIG_PER_LIMB 9
104 # define MAX_FAC_PER_LIMB 1000000000UL
105 #elif BITS_PER_MP_LIMB == 64
106 # define MAX_DIG_PER_LIMB 19
107 # define MAX_FAC_PER_LIMB 10000000000000000000ULL
109 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
112 #define _tens_in_limb __quadmath_tens_in_limb
113 extern const mp_limb_t _tens_in_limb
[MAX_DIG_PER_LIMB
+ 1] attribute_hidden
;
116 #define howmany(x,y) (((x)+((y)-1))/(y))
118 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
120 #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
121 #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
122 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
124 #define RETURN(val,end) \
125 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
126 return val; } while (0)
128 /* Maximum size necessary for mpn integers to hold floating point numbers. */
129 #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
131 /* Declare an mpn integer variable that big. */
132 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
133 /* Copy an mpn integer value. */
134 #define MPN_ASSIGN(dst, src) \
135 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
138 /* Return a floating point number of the needed type according to the given
139 multi-precision number after possible rounding. */
141 round_and_return (mp_limb_t
*retval
, int exponent
, int negative
,
142 mp_limb_t round_limb
, mp_size_t round_bit
, int more_bits
)
144 if (exponent
< MIN_EXP
- 1)
146 mp_size_t shift
= MIN_EXP
- 1 - exponent
;
148 if (shift
> MANT_DIG
)
150 #if defined HAVE_ERRNO_H && defined EDOM
156 more_bits
|= (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0;
157 if (shift
== MANT_DIG
)
158 /* This is a special case to handle the very seldom case where
159 the mantissa will be empty after the shift. */
163 round_limb
= retval
[RETURN_LIMB_SIZE
- 1];
164 round_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
165 for (i
= 0; i
< RETURN_LIMB_SIZE
; ++i
)
166 more_bits
|= retval
[i
] != 0;
167 MPN_ZERO (retval
, RETURN_LIMB_SIZE
);
169 else if (shift
>= BITS_PER_MP_LIMB
)
173 round_limb
= retval
[(shift
- 1) / BITS_PER_MP_LIMB
];
174 round_bit
= (shift
- 1) % BITS_PER_MP_LIMB
;
175 for (i
= 0; i
< (shift
- 1) / BITS_PER_MP_LIMB
; ++i
)
176 more_bits
|= retval
[i
] != 0;
177 more_bits
|= ((round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1))
180 (void) mpn_rshift (retval
, &retval
[shift
/ BITS_PER_MP_LIMB
],
181 RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
),
182 shift
% BITS_PER_MP_LIMB
);
183 MPN_ZERO (&retval
[RETURN_LIMB_SIZE
- (shift
/ BITS_PER_MP_LIMB
)],
184 shift
/ BITS_PER_MP_LIMB
);
188 round_limb
= retval
[0];
189 round_bit
= shift
- 1;
190 (void) mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
192 /* This is a hook for the m68k long double format, where the
193 exponent bias is the same for normalized and denormalized
196 # define DENORM_EXP (MIN_EXP - 2)
198 exponent
= DENORM_EXP
;
199 #if defined HAVE_ERRNO_H && defined ERANGE
204 if ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
205 && (more_bits
|| (retval
[0] & 1) != 0
206 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
208 mp_limb_t cy
= mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
210 if (((MANT_DIG
% BITS_PER_MP_LIMB
) == 0 && cy
) ||
211 ((MANT_DIG
% BITS_PER_MP_LIMB
) != 0 &&
212 (retval
[RETURN_LIMB_SIZE
- 1]
213 & (((mp_limb_t
) 1) << (MANT_DIG
% BITS_PER_MP_LIMB
))) != 0))
216 (void) mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, 1);
217 retval
[RETURN_LIMB_SIZE
- 1]
218 |= ((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
);
220 else if (exponent
== DENORM_EXP
221 && (retval
[RETURN_LIMB_SIZE
- 1]
222 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
224 /* The number was denormalized but now normalized. */
225 exponent
= MIN_EXP
- 1;
228 if (exponent
> MAX_EXP
)
229 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
231 return MPN2FLOAT (retval
, exponent
, negative
);
235 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
236 into N. Return the size of the number limbs in NSIZE at the first
237 character od the string that is not part of the integer as the function
238 value. If the EXPONENT is small enough to be taken as an additional
239 factor for the resulting number (see code) multiply by it. */
240 static const STRING_TYPE
*
241 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
243 #ifndef USE_WIDE_CHAR
244 , const char *decimal
, size_t decimal_len
, const char *thousands
249 /* Number of digits for actual limb. */
258 if (cnt
== MAX_DIG_PER_LIMB
)
268 cy
= mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
269 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
280 /* There might be thousands separators or radix characters in
281 the string. But these all can be ignored because we know the
282 format of the number is correct and we have an exact number
283 of characters to read. */
285 if (*str
< L
'0' || *str
> L
'9')
288 if (*str
< '0' || *str
> '9')
291 if (thousands
!= NULL
&& *str
== *thousands
292 && ({ for (inner
= 1; thousands
[inner
] != '\0'; ++inner
)
293 if (thousands
[inner
] != str
[inner
])
295 thousands
[inner
] == '\0'; }))
301 low
= low
* 10 + *str
++ - L_('0');
304 while (--digcnt
> 0);
306 if (*exponent
> 0 && cnt
+ *exponent
<= MAX_DIG_PER_LIMB
)
308 low
*= _tens_in_limb
[*exponent
];
309 start
= _tens_in_limb
[cnt
+ *exponent
];
313 start
= _tens_in_limb
[cnt
];
323 cy
= mpn_mul_1 (n
, n
, *nsize
, start
);
324 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
333 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
334 with the COUNT most significant bits of LIMB.
336 Tege doesn't like this function so I have to write it here myself. :)
339 __attribute ((always_inline
))
340 mpn_lshift_1 (mp_limb_t
*ptr
, mp_size_t size
, unsigned int count
,
343 if (__builtin_constant_p (count
) && count
== BITS_PER_MP_LIMB
)
345 /* Optimize the case of shifting by exactly a word:
346 just copy words, with no actual bit-shifting. */
348 for (i
= size
- 1; i
> 0; --i
)
354 (void) mpn_lshift (ptr
, ptr
, size
, count
);
355 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
360 #define INTERNAL(x) INTERNAL1(x)
361 #define INTERNAL1(x) __##x##_internal
362 #ifndef ____STRTOF_INTERNAL
363 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
366 /* This file defines a function to check for correct grouping. */
367 #include "grouping.h"
370 /* Return a floating point number with the value of the given string NPTR.
371 Set *ENDPTR to the character after the last used one. If the number is
372 smaller than the smallest representable number, set `errno' to ERANGE and
373 return 0.0. If the number is too big to be represented, set `errno' to
374 ERANGE and return HUGE_VAL with the appropriate sign. */
376 ____STRTOF_INTERNAL (nptr
, endptr
, group
)
377 const STRING_TYPE
*nptr
;
378 STRING_TYPE
**endptr
;
381 int negative
; /* The sign of the number. */
382 MPN_VAR (num
); /* MP representation of the number. */
383 int exponent
; /* Exponent of the number. */
385 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
388 /* When we have to compute fractional digits we form a fraction with a
389 second multi-precision number (and we sometimes need a second for
390 temporary results). */
393 /* Representation for the return value. */
394 mp_limb_t retval
[RETURN_LIMB_SIZE
];
395 /* Number of bits currently in result value. */
398 /* Running pointer after the last character processed in the string. */
399 const STRING_TYPE
*cp
, *tp
;
400 /* Start of significant part of the number. */
401 const STRING_TYPE
*startp
, *start_of_digits
;
402 /* Points at the character following the integer and fractional digits. */
403 const STRING_TYPE
*expp
;
404 /* Total number of digit and number of digits in integer part. */
405 int dig_no
, int_no
, lead_zero
;
406 /* Contains the last character read. */
409 /* The radix character of the current locale. */
416 /* The thousands character of the current locale. */
418 wchar_t thousands
= L
'\0';
420 const char *thousands
= NULL
;
422 /* The numeric grouping specification of the current locale,
423 in the format described in <locale.h>. */
424 const char *grouping
;
425 /* Used in several places. */
428 #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
429 const struct lconv
*lc
= localeconv ();
432 if (__builtin_expect (group
, 0))
434 #ifdef USE_NL_LANGINFO
435 grouping
= nl_langinfo (GROUPING
);
436 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
440 /* Figure out the thousands separator character. */
442 thousands
= nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC
);
443 if (thousands
== L
'\0')
446 thousands
= nl_langinfo (THOUSANDS_SEP
);
447 if (*thousands
== '\0')
454 #elif defined USE_LOCALECONV
455 grouping
= lc
->grouping
;
456 if (grouping
== NULL
|| *grouping
<= 0 || *grouping
== CHAR_MAX
)
460 /* Figure out the thousands separator character. */
461 thousands
= lc
->thousands_sep
;
462 if (thousands
== NULL
|| *thousands
== '\0')
475 /* Find the locale's decimal point character. */
477 decimal
= nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC
);
478 assert (decimal
!= L
'\0');
479 # define decimal_len 1
481 #ifdef USE_NL_LANGINFO
482 decimal
= nl_langinfo (DECIMAL_POINT
);
483 decimal_len
= strlen (decimal
);
484 assert (decimal_len
> 0);
485 #elif defined USE_LOCALECONV
486 decimal
= lc
->decimal_point
;
487 if (decimal
== NULL
|| *decimal
== '\0')
489 decimal_len
= strlen (decimal
);
496 /* Prepare number representation. */
501 /* Parse string to get maximal legal prefix. We need the number of
502 characters of the integer part, the fractional part and the exponent. */
504 /* Ignore leading white space. */
509 /* Get sign of the result. */
515 else if (c
== L_('+'))
518 /* Return 0.0 if no legal string is found.
519 No character is used even if a sign was found. */
521 if (c
== (wint_t) decimal
522 && (wint_t) cp
[1] >= L
'0' && (wint_t) cp
[1] <= L
'9')
524 /* We accept it. This funny construct is here only to indent
525 the code correctly. */
528 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
529 if (cp
[cnt
] != decimal
[cnt
])
531 if (decimal
[cnt
] == '\0' && cp
[cnt
] >= '0' && cp
[cnt
] <= '9')
533 /* We accept it. This funny construct is here only to indent
534 the code correctly. */
537 else if (c
< L_('0') || c
> L_('9'))
539 /* Check for `INF' or `INFINITY'. */
540 CHAR_TYPE lowc
= TOLOWER_C (c
);
542 if (lowc
== L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
544 /* Return +/- infinity. */
546 *endptr
= (STRING_TYPE
*)
547 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
550 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
553 if (lowc
== L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
560 /* Match `(n-char-sequence-digit)'. */
563 const STRING_TYPE
*startp
= cp
;
566 while ((*cp
>= L_('0') && *cp
<= L_('9'))
567 || ({ CHAR_TYPE lo
= TOLOWER (*cp
);
568 lo
>= L_('a') && lo
<= L_('z'); })
572 /* The closing brace is missing. Only match the NAN
577 /* This is a system-dependent way to specify the
578 bitmask used for the NaN. We expect it to be
579 a number which is put in the mantissa of the
582 unsigned long long int mant
;
584 mant
= STRTOULL (startp
+ 1, &endp
, 0);
586 SET_MANTISSA (retval
, mant
);
588 /* Consume the closing brace. */
594 *endptr
= (STRING_TYPE
*) cp
;
599 /* It is really a text we do not recognize. */
603 /* First look whether we are faced with a hexadecimal number. */
604 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
606 /* Okay, it is a hexa-decimal number. Remember this and skip
607 the characters. BTW: hexadecimal numbers must not be
615 /* Record the start of the digits, in case we will check their grouping. */
616 start_of_digits
= startp
= cp
;
618 /* Ignore leading zeroes. This helps us to avoid useless computations. */
620 while (c
== L
'0' || ((wint_t) thousands
!= L
'\0' && c
== (wint_t) thousands
))
623 if (__builtin_expect (thousands
== NULL
, 1))
628 /* We also have the multibyte thousands string. */
633 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
634 if (thousands
[cnt
] != cp
[cnt
])
636 if (thousands
[cnt
] != '\0')
645 /* If no other digit but a '0' is found the result is 0.0.
646 Return current read pointer. */
647 CHAR_TYPE lowc
= TOLOWER (c
);
648 if (!((c
>= L_('0') && c
<= L_('9'))
649 || (base
== 16 && lowc
>= L_('a') && lowc
<= L_('f'))
652 c
== (wint_t) decimal
654 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
655 if (decimal
[cnt
] != cp
[cnt
])
657 decimal
[cnt
] == '\0'; })
659 /* '0x.' alone is not a valid hexadecimal number.
660 '.' alone is not valid either, but that has been checked
663 || cp
!= start_of_digits
664 || (cp
[decimal_len
] >= L_('0') && cp
[decimal_len
] <= L_('9'))
665 || ({ CHAR_TYPE lo
= TOLOWER (cp
[decimal_len
]);
666 lo
>= L_('a') && lo
<= L_('f'); })))
667 || (base
== 16 && (cp
!= start_of_digits
669 || (base
!= 16 && lowc
== L_('e'))))
672 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
675 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
678 /* If TP is at the start of the digits, there was no correctly
679 grouped prefix of the string; so no number found. */
680 RETURN (negative
? -0.0 : 0.0,
681 tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
684 /* Remember first significant digit and read following characters until the
685 decimal point, exponent character or any non-FP number character. */
690 if ((c
>= L_('0') && c
<= L_('9'))
692 && ({ CHAR_TYPE lo
= TOLOWER (c
);
693 lo
>= L_('a') && lo
<= L_('f'); })))
698 if (__builtin_expect ((wint_t) thousands
== L
'\0', 1)
699 || c
!= (wint_t) thousands
)
700 /* Not a digit or separator: end of the integer part. */
703 if (__builtin_expect (thousands
== NULL
, 1))
707 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
708 if (thousands
[cnt
] != cp
[cnt
])
710 if (thousands
[cnt
] != '\0')
719 if (__builtin_expect (grouping
!= NULL
, 0) && cp
> start_of_digits
)
721 /* Check the grouping of the digits. */
723 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
726 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
731 /* Less than the entire string was correctly grouped. */
733 if (tp
== start_of_digits
)
734 /* No valid group of numbers at all: no valid number. */
738 /* The number is validly grouped, but consists
739 only of zeroes. The whole value is zero. */
740 RETURN (negative
? -0.0 : 0.0, tp
);
742 /* Recompute DIG_NO so we won't read more digits than
743 are properly grouped. */
746 for (tp
= startp
; tp
< cp
; ++tp
)
747 if (*tp
>= L_('0') && *tp
<= L_('9'))
757 /* We have the number of digits in the integer part. Whether these
758 are all or any is really a fractional digit will be decided
761 lead_zero
= int_no
== 0 ? -1 : 0;
763 /* Read the fractional digits. A special case are the 'american
764 style' numbers like `16.' i.e. with decimal point but without
768 c
== (wint_t) decimal
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 && ({ CHAR_TYPE lo
= TOLOWER (c
);
781 lo
>= L_('a') && lo
<= L_('f'); })))
783 if (c
!= L_('0') && lead_zero
== -1)
784 lead_zero
= dig_no
- int_no
;
790 /* Remember start of exponent (if any). */
795 if ((base
== 16 && lowc
== L_('p'))
796 || (base
!= 16 && lowc
== L_('e')))
798 int exp_negative
= 0;
806 else if (c
== L_('+'))
809 if (c
>= L_('0') && c
<= L_('9'))
813 /* Get the exponent limit. */
815 exp_limit
= (exp_negative
?
816 -MIN_EXP
+ MANT_DIG
+ 4 * int_no
:
817 MAX_EXP
- 4 * int_no
+ 4 * lead_zero
+ 3);
819 exp_limit
= (exp_negative
?
820 -MIN_10_EXP
+ MANT_DIG
+ int_no
:
821 MAX_10_EXP
- int_no
+ lead_zero
+ 1);
826 exponent
+= c
- L_('0');
828 if (__builtin_expect (exponent
> exp_limit
, 0))
829 /* The exponent is too large/small to represent a valid
834 /* We have to take care for special situation: a joker
835 might have written "0.0e100000" which is in fact
838 result
= negative
? -0.0 : 0.0;
841 /* Overflow or underflow. */
842 #if defined HAVE_ERRNO_H && defined ERANGE
845 result
= (exp_negative
? (negative
? -0.0 : 0.0) :
846 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
849 /* Accept all following digits as part of the exponent. */
852 while (*cp
>= L_('0') && *cp
<= L_('9'));
860 while (c
>= L_('0') && c
<= L_('9'));
863 exponent
= -exponent
;
869 /* We don't want to have to work with trailing zeroes after the radix. */
872 while (expp
[-1] == L_('0'))
877 assert (dig_no
>= int_no
);
880 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
883 while (! (base
== 16 ? ISXDIGIT (expp
[-1]) : ISDIGIT (expp
[-1])))
886 if (expp
[-1] != L_('0'))
892 exponent
+= base
== 16 ? 4 : 1;
894 while (dig_no
> 0 && exponent
< 0);
898 /* The whole string is parsed. Store the address of the next character. */
900 *endptr
= (STRING_TYPE
*) cp
;
903 return negative
? -0.0 : 0.0;
907 /* Find the decimal point */
909 while (*startp
!= decimal
)
914 if (*startp
== decimal
[0])
916 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
917 if (decimal
[cnt
] != startp
[cnt
])
919 if (decimal
[cnt
] == '\0')
925 startp
+= lead_zero
+ decimal_len
;
926 exponent
-= base
== 16 ? 4 * lead_zero
: lead_zero
;
930 /* If the BASE is 16 we can use a simpler algorithm. */
933 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
934 4, 4, 4, 4, 4, 4, 4, 4 };
935 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
936 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
939 while (!ISXDIGIT (*startp
))
941 while (*startp
== L_('0'))
943 if (ISDIGIT (*startp
))
944 val
= *startp
++ - L_('0');
946 val
= 10 + TOLOWER (*startp
++) - L_('a');
948 /* We cannot have a leading zero. */
951 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
953 /* We don't have to care for wrapping. This is the normal
954 case so we add the first clause in the `if' expression as
955 an optimization. It is a compile-time constant and so does
956 not cost anything. */
957 retval
[idx
] = val
<< (pos
- bits
+ 1);
962 retval
[idx
--] = val
>> (bits
- pos
- 1);
963 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
964 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
967 /* Adjust the exponent for the bits we are shifting in. */
968 exponent
+= bits
- 1 + (int_no
- 1) * 4;
970 while (--dig_no
> 0 && idx
>= 0)
972 if (!ISXDIGIT (*startp
))
973 startp
+= decimal_len
;
974 if (ISDIGIT (*startp
))
975 val
= *startp
++ - L_('0');
977 val
= 10 + TOLOWER (*startp
++) - L_('a');
981 retval
[idx
] |= val
<< (pos
- 4 + 1);
986 retval
[idx
--] |= val
>> (4 - pos
- 1);
987 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
989 return round_and_return (retval
, exponent
, negative
, val
,
990 BITS_PER_MP_LIMB
- 1, dig_no
> 0);
993 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
997 /* We ran out of digits. */
998 MPN_ZERO (retval
, idx
);
1000 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
1003 /* Now we have the number of digits in total and the integer digits as well
1004 as the exponent and its sign. We can decide whether the read digits are
1005 really integer digits or belong to the fractional part; i.e. we normalize
1008 register int incr
= (exponent
< 0 ? MAX (-int_no
, exponent
)
1009 : MIN (dig_no
- int_no
, exponent
));
1014 if (__builtin_expect (int_no
+ exponent
> MAX_10_EXP
+ 1, 0))
1016 #if defined HAVE_ERRNO_H && defined ERANGE
1019 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1022 if (__builtin_expect (exponent
< MIN_10_EXP
- (DIG
+ 1), 0))
1024 #if defined HAVE_ERRNO_H && defined ERANGE
1027 return negative
? -0.0 : 0.0;
1032 /* Read the integer part as a multi-precision number to NUM. */
1033 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1034 #ifndef USE_WIDE_CHAR
1035 , decimal
, decimal_len
, thousands
1041 /* We now multiply the gained number by the given power of ten. */
1042 mp_limb_t
*psrc
= num
;
1043 mp_limb_t
*pdest
= den
;
1045 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1049 if ((exponent
& expbit
) != 0)
1051 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1055 /* FIXME: not the whole multiplication has to be
1056 done. If we have the needed number of bits we
1057 only need the information whether more non-zero
1059 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1060 cy
= mpn_mul (pdest
, psrc
, numsize
,
1061 &__tens
[ttab
->arrayoff
1062 + _FPIO_CONST_OFFSET
],
1065 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1066 + _FPIO_CONST_OFFSET
],
1067 size
, psrc
, numsize
);
1071 (void) SWAP (psrc
, pdest
);
1076 while (exponent
!= 0);
1079 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1082 /* Determine how many bits of the result we already have. */
1083 count_leading_zeros (bits
, num
[numsize
- 1]);
1084 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1086 /* Now we know the exponent of the number in base two.
1087 Check it against the maximum possible exponent. */
1088 if (__builtin_expect (bits
> MAX_EXP
, 0))
1090 #if defined HAVE_ERRNO_H && defined ERANGE
1093 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1096 /* We have already the first BITS bits of the result. Together with
1097 the information whether more non-zero bits follow this is enough
1098 to determine the result. */
1099 if (bits
> MANT_DIG
)
1102 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1103 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1104 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1106 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1110 memcpy (retval
, &num
[least_idx
],
1111 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1114 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1115 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1117 << (BITS_PER_MP_LIMB
- least_bit
));
1118 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1119 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1122 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1123 for (i
= 0; num
[i
] == 0; ++i
)
1126 return round_and_return (retval
, bits
- 1, negative
,
1127 num
[round_idx
], round_bit
,
1128 int_no
< dig_no
|| i
< round_idx
);
1131 else if (dig_no
== int_no
)
1133 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1134 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1136 if (target_bit
== is_bit
)
1138 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1139 numsize
* sizeof (mp_limb_t
));
1140 /* FIXME: the following loop can be avoided if we assume a
1141 maximal MANT_DIG value. */
1142 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1144 else if (target_bit
> is_bit
)
1146 (void) mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1147 num
, numsize
, target_bit
- is_bit
);
1148 /* FIXME: the following loop can be avoided if we assume a
1149 maximal MANT_DIG value. */
1150 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1155 assert (numsize
< RETURN_LIMB_SIZE
);
1157 cy
= mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1158 num
, numsize
, is_bit
- target_bit
);
1159 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1160 /* FIXME: the following loop can be avoided if we assume a
1161 maximal MANT_DIG value. */
1162 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1165 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1169 /* Store the bits we already have. */
1170 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1171 #if RETURN_LIMB_SIZE > 1
1172 if (numsize
< RETURN_LIMB_SIZE
)
1173 # if RETURN_LIMB_SIZE == 2
1174 retval
[numsize
] = 0;
1176 MPN_ZERO (retval
+ numsize
, RETURN_LIMB_SIZE
- numsize
);
1181 /* We have to compute at least some of the fractional digits. */
1183 /* We construct a fraction and the result of the division gives us
1184 the needed digits. The denominator is 1.0 multiplied by the
1185 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1186 123e-6 gives 123 / 1000000. */
1192 mp_limb_t
*psrc
= den
;
1193 mp_limb_t
*pdest
= num
;
1194 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1196 assert (dig_no
> int_no
&& exponent
<= 0);
1199 /* For the fractional part we need not process too many digits. One
1200 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1202 digits we should have enough bits for the result. The remaining
1203 decimal digits give us the information that more bits are following.
1204 This can be used while rounding. (Two added as a safety margin.) */
1205 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 2)
1207 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 2;
1213 neg_exp
= dig_no
- int_no
- exponent
;
1215 /* Construct the denominator. */
1220 if ((neg_exp
& expbit
) != 0)
1227 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1228 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1229 densize
* sizeof (mp_limb_t
));
1233 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1234 + _FPIO_CONST_OFFSET
],
1235 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1237 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1240 (void) SWAP (psrc
, pdest
);
1246 while (neg_exp
!= 0);
1249 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1251 /* Read the fractional digits from the string. */
1252 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1253 #ifndef USE_WIDE_CHAR
1254 , decimal
, decimal_len
, thousands
1258 /* We now have to shift both numbers so that the highest bit in the
1259 denominator is set. In the same process we copy the numerator to
1260 a high place in the array so that the division constructs the wanted
1261 digits. This is done by a "quasi fix point" number representation.
1263 num: ddddddddddd . 0000000000000000000000
1265 den: ddddddddddd n >= m
1269 count_leading_zeros (cnt
, den
[densize
- 1]);
1273 /* Don't call `mpn_shift' with a count of zero since the specification
1274 does not allow this. */
1275 (void) mpn_lshift (den
, den
, densize
, cnt
);
1276 cy
= mpn_lshift (num
, num
, numsize
, cnt
);
1278 num
[numsize
++] = cy
;
1281 /* Now we are ready for the division. But it is not necessary to
1282 do a full multi-precision division because we only need a small
1283 number of bits for the result. So we do not use mpn_divmod
1284 here but instead do the division here by hand and stop whenever
1285 the needed number of bits is reached. The code itself comes
1286 from the GNU MP Library by Torbj\"orn Granlund. */
1294 mp_limb_t d
, n
, quot
;
1299 assert (numsize
== 1 && n
< d
);
1303 udiv_qrnnd (quot
, n
, n
, 0, d
);
1310 cnt = BITS_PER_MP_LIMB; \
1312 count_leading_zeros (cnt, quot); \
1314 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1316 used = MANT_DIG + cnt; \
1317 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1318 bits = MANT_DIG + 1; \
1322 /* Note that we only clear the second element. */ \
1323 /* The conditional is determined at compile time. */ \
1324 if (RETURN_LIMB_SIZE > 1) \
1330 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1331 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1335 used = MANT_DIG - bits; \
1337 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1339 bits += BITS_PER_MP_LIMB
1343 while (bits
<= MANT_DIG
);
1345 return round_and_return (retval
, exponent
- 1, negative
,
1346 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1347 more_bits
|| n
!= 0);
1351 mp_limb_t d0
, d1
, n0
, n1
;
1358 if (numsize
< densize
)
1362 /* The numerator of the number occupies fewer bits than
1363 the denominator but the one limb is bigger than the
1364 high limb of the numerator. */
1371 exponent
-= BITS_PER_MP_LIMB
;
1374 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1375 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1376 BITS_PER_MP_LIMB
, 0);
1379 used
= MANT_DIG
- bits
;
1381 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1383 bits
+= BITS_PER_MP_LIMB
;
1395 while (bits
<= MANT_DIG
)
1401 /* QUOT should be either 111..111 or 111..110. We need
1402 special treatment of this rare case as normal division
1403 would give overflow. */
1404 quot
= ~(mp_limb_t
) 0;
1407 if (r
< d1
) /* Carry in the addition? */
1409 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1412 n1
= d0
- (d0
!= 0);
1417 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1418 umul_ppmm (n1
, n0
, d0
, quot
);
1422 if (n1
> r
|| (n1
== r
&& n0
> 0))
1424 /* The estimated QUOT was too large. */
1427 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1429 if (r
>= d1
) /* If not carry, test QUOT again. */
1432 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1438 return round_and_return (retval
, exponent
- 1, negative
,
1439 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1440 more_bits
|| n1
!= 0 || n0
!= 0);
1445 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1449 dX
= den
[densize
- 1];
1450 d1
= den
[densize
- 2];
1452 /* The division does not work if the upper limb of the two-limb
1453 numerator is greater than the denominator. */
1454 if (mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1457 if (numsize
< densize
)
1459 mp_size_t empty
= densize
- numsize
;
1463 exponent
-= empty
* BITS_PER_MP_LIMB
;
1466 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1468 /* We make a difference here because the compiler
1469 cannot optimize the `else' case that good and
1470 this reflects all currently used FLOAT types
1471 and GMP implementations. */
1472 #if RETURN_LIMB_SIZE <= 2
1473 assert (empty
== 1);
1474 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1475 BITS_PER_MP_LIMB
, 0);
1477 for (i
= RETURN_LIMB_SIZE
- 1; i
>= empty
; --i
)
1478 retval
[i
] = retval
[i
- empty
];
1485 used
= MANT_DIG
- bits
;
1486 if (used
>= BITS_PER_MP_LIMB
)
1489 (void) mpn_lshift (&retval
[used
1490 / BITS_PER_MP_LIMB
],
1493 - used
/ BITS_PER_MP_LIMB
),
1494 used
% BITS_PER_MP_LIMB
);
1495 for (i
= used
/ BITS_PER_MP_LIMB
- 1; i
>= 0; --i
)
1499 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1501 bits
+= empty
* BITS_PER_MP_LIMB
;
1503 for (i
= numsize
; i
> 0; --i
)
1504 num
[i
+ empty
] = num
[i
- 1];
1505 MPN_ZERO (num
, empty
+ 1);
1510 assert (numsize
== densize
);
1511 for (i
= numsize
; i
> 0; --i
)
1512 num
[i
] = num
[i
- 1];
1518 while (bits
<= MANT_DIG
)
1521 /* This might over-estimate QUOT, but it's probably not
1522 worth the extra code here to find out. */
1523 quot
= ~(mp_limb_t
) 0;
1528 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1529 umul_ppmm (n1
, n0
, d1
, quot
);
1531 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1535 if (r
< dX
) /* I.e. "carry in previous addition?" */
1542 /* Possible optimization: We already have (q * n0) and (1 * n1)
1543 after the calculation of QUOT. Taking advantage of this, we
1544 could make this loop make two iterations less. */
1546 cy
= mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1548 if (num
[densize
] != cy
)
1550 cy
= mpn_add_n (num
, num
, den
, densize
);
1554 n0
= num
[densize
] = num
[densize
- 1];
1555 for (i
= densize
- 1; i
> 0; --i
)
1556 num
[i
] = num
[i
- 1];
1561 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1563 return round_and_return (retval
, exponent
- 1, negative
,
1564 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1565 more_bits
|| i
>= 0);