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,2001 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'))
577 /* Check for `INF' or `INFINITY'. */
578 if (TOLOWER (c
) == L_('i')
579 && ((STRNCASECMP (cp
, L_("inf"), 3) == 0 && (matched
= 3))
580 || (STRNCASECMP (cp
, L_("infinity"), 8) == 0 && (matched
= 8))))
582 /* Return +/- infinity. */
584 *endptr
= (STRING_TYPE
*) (cp
+ matched
);
586 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
589 if (TOLOWER (c
) == L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
596 /* Match `(n-char-sequence-digit)'. */
599 const STRING_TYPE
*startp
= cp
;
602 while ((*cp
>= L_('0') && *cp
<= L_('9'))
603 || (TOLOWER (*cp
) >= L_('a') && TOLOWER (*cp
) <= L_('z'))
607 /* The closing brace is missing. Only match the NAN
612 /* This is a system-dependent way to specify the
613 bitmask used for the NaN. We expect it to be
614 a number which is put in the mantissa of the
617 unsigned long long int mant
;
619 mant
= STRTOULL (startp
+ 1, &endp
, 0);
621 SET_MANTISSA (retval
, mant
);
626 *endptr
= (STRING_TYPE
*) cp
;
631 /* It is really a text we do not recognize. */
635 /* First look whether we are faced with a hexadecimal number. */
636 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
638 /* Okay, it is a hexa-decimal number. Remember this and skip
639 the characters. BTW: hexadecimal numbers must not be
647 /* Record the start of the digits, in case we will check their grouping. */
648 start_of_digits
= startp
= cp
;
650 /* Ignore leading zeroes. This helps us to avoid useless computations. */
652 while (c
== L
'0' || (thousands
!= L
'\0' && c
== thousands
))
655 if (thousands
== NULL
)
660 /* We also have the multibyte thousands string. */
665 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
666 if (c
!= thousands
[cnt
])
668 if (thousands
[cnt
] != '\0')
676 /* If no other digit but a '0' is found the result is 0.0.
677 Return current read pointer. */
678 if ((c
< L_('0') || c
> L_('9'))
679 && (base
== 16 && (c
< TOLOWER (L_('a')) || c
> TOLOWER (L_('f'))))
683 && ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
684 if (decimal
[cnt
] != cp
[cnt
])
686 decimal
[cnt
] != '\0'; })
688 && (base
== 16 && (cp
== start_of_digits
|| TOLOWER (c
) != L_('p')))
689 && (base
!= 16 && TOLOWER (c
) != L_('e')))
691 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
692 /* If TP is at the start of the digits, there was no correctly
693 grouped prefix of the string; so no number found. */
694 RETURN (0.0, tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
697 /* Remember first significant digit and read following characters until the
698 decimal point, exponent character or any non-FP number character. */
703 if ((c
>= L_('0') && c
<= L_('9'))
704 || (base
== 16 && TOLOWER (c
) >= L_('a') && TOLOWER (c
) <= L_('f')))
709 if (thousands
== L
'\0' || c
!= thousands
)
710 /* Not a digit or separator: end of the integer part. */
713 if (thousands
== NULL
)
717 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
718 if (thousands
[cnt
] != cp
[cnt
])
720 if (thousands
[cnt
] != '\0')
728 if (grouping
&& dig_no
> 0)
730 /* Check the grouping of the digits. */
731 tp
= correctly_grouped_prefix (start_of_digits
, cp
, thousands
, grouping
);
734 /* Less than the entire string was correctly grouped. */
736 if (tp
== start_of_digits
)
737 /* No valid group of numbers at all: no valid number. */
741 /* The number is validly grouped, but consists
742 only of zeroes. The whole value is zero. */
745 /* Recompute DIG_NO so we won't read more digits than
746 are properly grouped. */
749 for (tp
= startp
; tp
< cp
; ++tp
)
750 if (*tp
>= L_('0') && *tp
<= L_('9'))
760 /* We have the number digits in the integer part. Whether these are all or
761 any is really a fractional digit will be decided later. */
763 lead_zero
= int_no
== 0 ? -1 : 0;
765 /* Read the fractional digits. A special case are the 'american style'
766 numbers like `16.' i.e. with decimal but without trailing digits. */
771 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
772 if (decimal
[cnt
] != cp
[cnt
])
774 decimal
[cnt
] == '\0'; })
780 while ((c
>= L_('0') && c
<= L_('9')) ||
781 (base
== 16 && TOLOWER (c
) >= L_('a') && TOLOWER (c
) <= L_('f')))
783 if (c
!= L_('0') && lead_zero
== -1)
784 lead_zero
= dig_no
- int_no
;
790 /* Remember start of exponent (if any). */
794 if ((base
== 16 && TOLOWER (c
) == L_('p'))
795 || (base
!= 16 && TOLOWER (c
) == L_('e')))
797 int exp_negative
= 0;
805 else if (c
== L_('+'))
808 if (c
>= L_('0') && c
<= L_('9'))
812 /* Get the exponent limit. */
814 exp_limit
= (exp_negative
?
815 -MIN_EXP
+ MANT_DIG
+ 4 * int_no
:
816 MAX_EXP
- 4 * int_no
+ lead_zero
);
818 exp_limit
= (exp_negative
?
819 -MIN_10_EXP
+ MANT_DIG
+ int_no
:
820 MAX_10_EXP
- int_no
+ lead_zero
);
826 if (exponent
> exp_limit
)
827 /* The exponent is too large/small to represent a valid
832 /* We have to take care for special situation: a joker
833 might have written "0.0e100000" which is in fact
836 result
= negative
? -0.0 : 0.0;
839 /* Overflow or underflow. */
840 __set_errno (ERANGE
);
841 result
= (exp_negative
? 0.0 :
842 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
845 /* Accept all following digits as part of the exponent. */
848 while (*cp
>= L_('0') && *cp
<= L_('9'));
854 exponent
+= c
- L_('0');
857 while (c
>= L_('0') && c
<= L_('9'));
860 exponent
= -exponent
;
866 /* We don't want to have to work with trailing zeroes after the radix. */
869 while (expp
[-1] == L_('0'))
874 assert (dig_no
>= int_no
);
877 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
880 while (expp
[-1] < L_('0') || expp
[-1] > L_('9'))
883 if (expp
[-1] != L_('0'))
891 while (dig_no
> 0 && exponent
< 0);
895 /* The whole string is parsed. Store the address of the next character. */
897 *endptr
= (STRING_TYPE
*) cp
;
900 return negative
? -0.0 : 0.0;
904 /* Find the decimal point */
906 while (*startp
!= decimal
)
911 if (*startp
== decimal
[0])
913 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
914 if (decimal
[cnt
] != startp
[cnt
])
916 if (decimal
[cnt
] == '\0')
922 startp
+= lead_zero
+ decimal_len
;
923 exponent
-= base
== 16 ? 4 * lead_zero
: lead_zero
;
927 /* If the BASE is 16 we can use a simpler algorithm. */
930 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
931 4, 4, 4, 4, 4, 4, 4, 4 };
932 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
933 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
936 while (!ISXDIGIT (*startp
))
938 while (*startp
== L_('0'))
940 if (ISDIGIT (*startp
))
941 val
= *startp
++ - L_('0');
943 val
= 10 + TOLOWER (*startp
++) - L_('a');
945 /* We cannot have a leading zero. */
948 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
950 /* We don't have to care for wrapping. This is the normal
951 case so we add the first clause in the `if' expression as
952 an optimization. It is a compile-time constant and so does
953 not cost anything. */
954 retval
[idx
] = val
<< (pos
- bits
+ 1);
959 retval
[idx
--] = val
>> (bits
- pos
- 1);
960 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
961 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
964 /* Adjust the exponent for the bits we are shifting in. */
965 exponent
+= bits
- 1 + (int_no
- 1) * 4;
967 while (--dig_no
> 0 && idx
>= 0)
969 if (!ISXDIGIT (*startp
))
970 startp
+= decimal_len
;
971 if (ISDIGIT (*startp
))
972 val
= *startp
++ - L_('0');
974 val
= 10 + TOLOWER (*startp
++) - L_('a');
978 retval
[idx
] |= val
<< (pos
- 4 + 1);
983 retval
[idx
--] |= val
>> (4 - pos
- 1);
984 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
986 return round_and_return (retval
, exponent
, negative
, val
,
987 BITS_PER_MP_LIMB
- 1, dig_no
> 0);
990 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
994 /* We ran out of digits. */
995 MPN_ZERO (retval
, idx
);
997 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
1000 /* Now we have the number of digits in total and the integer digits as well
1001 as the exponent and its sign. We can decide whether the read digits are
1002 really integer digits or belong to the fractional part; i.e. we normalize
1005 register int incr
= (exponent
< 0 ? MAX (-int_no
, exponent
)
1006 : MIN (dig_no
- int_no
, exponent
));
1011 if (int_no
+ exponent
> MAX_10_EXP
+ 1)
1013 __set_errno (ERANGE
);
1014 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1017 if (exponent
< MIN_10_EXP
- (DIG
+ 1))
1019 __set_errno (ERANGE
);
1025 /* Read the integer part as a multi-precision number to NUM. */
1026 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1027 #ifndef USE_WIDE_CHAR
1028 , decimal
, decimal_len
, thousands
1034 /* We now multiply the gained number by the given power of ten. */
1035 mp_limb_t
*psrc
= num
;
1036 mp_limb_t
*pdest
= den
;
1038 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1042 if ((exponent
& expbit
) != 0)
1044 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1048 /* FIXME: not the whole multiplication has to be
1049 done. If we have the needed number of bits we
1050 only need the information whether more non-zero
1052 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1053 cy
= __mpn_mul (pdest
, psrc
, numsize
,
1054 &__tens
[ttab
->arrayoff
1055 + _FPIO_CONST_OFFSET
],
1058 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1059 + _FPIO_CONST_OFFSET
],
1060 size
, psrc
, numsize
);
1064 (void) SWAP (psrc
, pdest
);
1069 while (exponent
!= 0);
1072 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1075 /* Determine how many bits of the result we already have. */
1076 count_leading_zeros (bits
, num
[numsize
- 1]);
1077 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1079 /* Now we know the exponent of the number in base two.
1080 Check it against the maximum possible exponent. */
1083 __set_errno (ERANGE
);
1084 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
1087 /* We have already the first BITS bits of the result. Together with
1088 the information whether more non-zero bits follow this is enough
1089 to determine the result. */
1090 if (bits
> MANT_DIG
)
1093 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1094 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1095 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1097 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1101 memcpy (retval
, &num
[least_idx
],
1102 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1105 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1106 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1108 << (BITS_PER_MP_LIMB
- least_bit
));
1109 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1110 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1113 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1114 for (i
= 0; num
[i
] == 0; ++i
)
1117 return round_and_return (retval
, bits
- 1, negative
,
1118 num
[round_idx
], round_bit
,
1119 int_no
< dig_no
|| i
< round_idx
);
1122 else if (dig_no
== int_no
)
1124 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1125 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1127 if (target_bit
== is_bit
)
1129 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1130 numsize
* sizeof (mp_limb_t
));
1131 /* FIXME: the following loop can be avoided if we assume a
1132 maximal MANT_DIG value. */
1133 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1135 else if (target_bit
> is_bit
)
1137 (void) __mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1138 num
, numsize
, target_bit
- is_bit
);
1139 /* FIXME: the following loop can be avoided if we assume a
1140 maximal MANT_DIG value. */
1141 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1146 assert (numsize
< RETURN_LIMB_SIZE
);
1148 cy
= __mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1149 num
, numsize
, is_bit
- target_bit
);
1150 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1151 /* FIXME: the following loop can be avoided if we assume a
1152 maximal MANT_DIG value. */
1153 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1156 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1160 /* Store the bits we already have. */
1161 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1162 #if RETURN_LIMB_SIZE > 1
1163 if (numsize
< RETURN_LIMB_SIZE
)
1164 retval
[numsize
] = 0;
1168 /* We have to compute at least some of the fractional digits. */
1170 /* We construct a fraction and the result of the division gives us
1171 the needed digits. The denominator is 1.0 multiplied by the
1172 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1173 123e-6 gives 123 / 1000000. */
1179 mp_limb_t
*psrc
= den
;
1180 mp_limb_t
*pdest
= num
;
1181 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1183 assert (dig_no
> int_no
&& exponent
<= 0);
1186 /* For the fractional part we need not process too many digits. One
1187 decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
1189 digits we should have enough bits for the result. The remaining
1190 decimal digits give us the information that more bits are following.
1191 This can be used while rounding. (One added as a safety margin.) */
1192 if (dig_no
- int_no
> (MANT_DIG
- bits
+ 2) / 3 + 1)
1194 dig_no
= int_no
+ (MANT_DIG
- bits
+ 2) / 3 + 1;
1200 neg_exp
= dig_no
- int_no
- exponent
;
1202 /* Construct the denominator. */
1207 if ((neg_exp
& expbit
) != 0)
1214 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1215 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1216 densize
* sizeof (mp_limb_t
));
1220 cy
= __mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1221 + _FPIO_CONST_OFFSET
],
1222 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1224 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1227 (void) SWAP (psrc
, pdest
);
1233 while (neg_exp
!= 0);
1236 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1238 /* Read the fractional digits from the string. */
1239 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1240 #ifndef USE_WIDE_CHAR
1241 , decimal
, decimal_len
, thousands
1245 /* We now have to shift both numbers so that the highest bit in the
1246 denominator is set. In the same process we copy the numerator to
1247 a high place in the array so that the division constructs the wanted
1248 digits. This is done by a "quasi fix point" number representation.
1250 num: ddddddddddd . 0000000000000000000000
1252 den: ddddddddddd n >= m
1256 count_leading_zeros (cnt
, den
[densize
- 1]);
1260 /* Don't call `mpn_shift' with a count of zero since the specification
1261 does not allow this. */
1262 (void) __mpn_lshift (den
, den
, densize
, cnt
);
1263 cy
= __mpn_lshift (num
, num
, numsize
, cnt
);
1265 num
[numsize
++] = cy
;
1268 /* Now we are ready for the division. But it is not necessary to
1269 do a full multi-precision division because we only need a small
1270 number of bits for the result. So we do not use __mpn_divmod
1271 here but instead do the division here by hand and stop whenever
1272 the needed number of bits is reached. The code itself comes
1273 from the GNU MP Library by Torbj\"orn Granlund. */
1281 mp_limb_t d
, n
, quot
;
1286 assert (numsize
== 1 && n
< d
);
1290 udiv_qrnnd (quot
, n
, n
, 0, d
);
1297 cnt = BITS_PER_MP_LIMB; \
1299 count_leading_zeros (cnt, quot); \
1301 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1303 used = MANT_DIG + cnt; \
1304 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1305 bits = MANT_DIG + 1; \
1309 /* Note that we only clear the second element. */ \
1310 /* The conditional is determined at compile time. */ \
1311 if (RETURN_LIMB_SIZE > 1) \
1317 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1318 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1322 used = MANT_DIG - bits; \
1324 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1326 bits += BITS_PER_MP_LIMB
1330 while (bits
<= MANT_DIG
);
1332 return round_and_return (retval
, exponent
- 1, negative
,
1333 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1334 more_bits
|| n
!= 0);
1338 mp_limb_t d0
, d1
, n0
, n1
;
1345 if (numsize
< densize
)
1349 /* The numerator of the number occupies fewer bits than
1350 the denominator but the one limb is bigger than the
1351 high limb of the numerator. */
1358 exponent
-= BITS_PER_MP_LIMB
;
1361 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1362 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1363 BITS_PER_MP_LIMB
, 0);
1366 used
= MANT_DIG
- bits
;
1368 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1370 bits
+= BITS_PER_MP_LIMB
;
1382 while (bits
<= MANT_DIG
)
1388 /* QUOT should be either 111..111 or 111..110. We need
1389 special treatment of this rare case as normal division
1390 would give overflow. */
1391 quot
= ~(mp_limb_t
) 0;
1394 if (r
< d1
) /* Carry in the addition? */
1396 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1399 n1
= d0
- (d0
!= 0);
1404 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1405 umul_ppmm (n1
, n0
, d0
, quot
);
1409 if (n1
> r
|| (n1
== r
&& n0
> 0))
1411 /* The estimated QUOT was too large. */
1414 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1416 if (r
>= d1
) /* If not carry, test QUOT again. */
1419 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1425 return round_and_return (retval
, exponent
- 1, negative
,
1426 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1427 more_bits
|| n1
!= 0 || n0
!= 0);
1432 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1436 dX
= den
[densize
- 1];
1437 d1
= den
[densize
- 2];
1439 /* The division does not work if the upper limb of the two-limb
1440 numerator is greater than the denominator. */
1441 if (__mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1444 if (numsize
< densize
)
1446 mp_size_t empty
= densize
- numsize
;
1451 for (i
= numsize
; i
> 0; --i
)
1452 num
[i
+ empty
] = num
[i
- 1];
1453 MPN_ZERO (num
, empty
+ 1);
1454 exponent
-= empty
* BITS_PER_MP_LIMB
;
1458 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1460 /* We make a difference here because the compiler
1461 cannot optimize the `else' case that good and
1462 this reflects all currently used FLOAT types
1463 and GMP implementations. */
1465 #if RETURN_LIMB_SIZE <= 2
1466 assert (empty
== 1);
1467 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1468 BITS_PER_MP_LIMB
, 0);
1470 for (i
= RETURN_LIMB_SIZE
; i
> empty
; --i
)
1471 retval
[i
] = retval
[i
- empty
];
1473 #if RETURN_LIMB_SIZE > 1
1476 for (i
= numsize
; i
> 0; --i
)
1477 num
[i
+ empty
] = num
[i
- 1];
1478 MPN_ZERO (num
, empty
+ 1);
1482 used
= MANT_DIG
- bits
;
1483 if (used
>= BITS_PER_MP_LIMB
)
1486 (void) __mpn_lshift (&retval
[used
1487 / BITS_PER_MP_LIMB
],
1488 retval
, RETURN_LIMB_SIZE
,
1489 used
% BITS_PER_MP_LIMB
);
1490 for (i
= used
/ BITS_PER_MP_LIMB
; i
>= 0; --i
)
1494 __mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1496 bits
+= empty
* BITS_PER_MP_LIMB
;
1502 assert (numsize
== densize
);
1503 for (i
= numsize
; i
> 0; --i
)
1504 num
[i
] = num
[i
- 1];
1510 while (bits
<= MANT_DIG
)
1513 /* This might over-estimate QUOT, but it's probably not
1514 worth the extra code here to find out. */
1515 quot
= ~(mp_limb_t
) 0;
1520 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1521 umul_ppmm (n1
, n0
, d1
, quot
);
1523 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1527 if (r
< dX
) /* I.e. "carry in previous addition?" */
1534 /* Possible optimization: We already have (q * n0) and (1 * n1)
1535 after the calculation of QUOT. Taking advantage of this, we
1536 could make this loop make two iterations less. */
1538 cy
= __mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1540 if (num
[densize
] != cy
)
1542 cy
= __mpn_add_n (num
, num
, den
, densize
);
1546 n0
= num
[densize
] = num
[densize
- 1];
1547 for (i
= densize
- 1; i
> 0; --i
)
1548 num
[i
] = num
[i
- 1];
1553 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1555 return round_and_return (retval
, exponent
- 1, negative
,
1556 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1557 more_bits
|| i
>= 0);
1565 /* External user entry point. */
1568 #ifdef weak_function
1571 STRTOF (nptr
, endptr LOCALE_PARAM
)
1572 const STRING_TYPE
*nptr
;
1573 STRING_TYPE
**endptr
;
1576 return INTERNAL (STRTOF
) (nptr
, endptr
, 0 LOCALE_PARAM
);