BZ#11261 fixed by previous commit.
[glibc.git] / stdlib / strtod_l.c
blob47247b548ccc160d6aac4366e07525bb8e0d47ce
1 /* Convert string representing a number to float value, using given locale.
2 Copyright (C) 1997-2013 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
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, see
18 <http://www.gnu.org/licenses/>. */
20 #include <xlocale.h>
22 extern double ____strtod_l_internal (const char *, char **, int, __locale_t);
23 extern unsigned long long int ____strtoull_l_internal (const char *, char **,
24 int, int, __locale_t);
26 /* Configuration part. These macros are defined by `strtold.c',
27 `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
28 `long double' and `float' versions of the reader. */
29 #ifndef FLOAT
30 # include <math_ldbl_opt.h>
31 # define FLOAT double
32 # define FLT DBL
33 # ifdef USE_WIDE_CHAR
34 # define STRTOF wcstod_l
35 # define __STRTOF __wcstod_l
36 # else
37 # define STRTOF strtod_l
38 # define __STRTOF __strtod_l
39 # endif
40 # define MPN2FLOAT __mpn_construct_double
41 # define FLOAT_HUGE_VAL HUGE_VAL
42 # define SET_MANTISSA(flt, mant) \
43 do { union ieee754_double u; \
44 u.d = (flt); \
45 if ((mant & 0xfffffffffffffULL) == 0) \
46 mant = 0x8000000000000ULL; \
47 u.ieee.mantissa0 = ((mant) >> 32) & 0xfffff; \
48 u.ieee.mantissa1 = (mant) & 0xffffffff; \
49 (flt) = u.d; \
50 } while (0)
51 #endif
52 /* End of configuration part. */
54 #include <ctype.h>
55 #include <errno.h>
56 #include <float.h>
57 #include <ieee754.h>
58 #include "../locale/localeinfo.h"
59 #include <locale.h>
60 #include <math.h>
61 #include <stdlib.h>
62 #include <string.h>
63 #include <stdint.h>
64 #include <rounding-mode.h>
65 #include <tininess.h>
67 /* The gmp headers need some configuration frobs. */
68 #define HAVE_ALLOCA 1
70 /* Include gmp-mparam.h first, such that definitions of _SHORT_LIMB
71 and _LONG_LONG_LIMB in it can take effect into gmp.h. */
72 #include <gmp-mparam.h>
73 #include <gmp.h>
74 #include "gmp-impl.h"
75 #include "longlong.h"
76 #include "fpioconst.h"
78 #include <assert.h>
81 /* We use this code for the extended locale handling where the
82 function gets as an additional argument the locale which has to be
83 used. To access the values we have to redefine the _NL_CURRENT and
84 _NL_CURRENT_WORD macros. */
85 #undef _NL_CURRENT
86 #define _NL_CURRENT(category, item) \
87 (current->values[_NL_ITEM_INDEX (item)].string)
88 #undef _NL_CURRENT_WORD
89 #define _NL_CURRENT_WORD(category, item) \
90 ((uint32_t) current->values[_NL_ITEM_INDEX (item)].word)
92 #if defined _LIBC || defined HAVE_WCHAR_H
93 # include <wchar.h>
94 #endif
96 #ifdef USE_WIDE_CHAR
97 # include <wctype.h>
98 # define STRING_TYPE wchar_t
99 # define CHAR_TYPE wint_t
100 # define L_(Ch) L##Ch
101 # define ISSPACE(Ch) __iswspace_l ((Ch), loc)
102 # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
103 # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
104 # define TOLOWER(Ch) __towlower_l ((Ch), loc)
105 # define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
106 # define STRNCASECMP(S1, S2, N) \
107 __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
108 # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
109 #else
110 # define STRING_TYPE char
111 # define CHAR_TYPE char
112 # define L_(Ch) Ch
113 # define ISSPACE(Ch) __isspace_l ((Ch), loc)
114 # define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
115 # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
116 # define TOLOWER(Ch) __tolower_l ((Ch), loc)
117 # define TOLOWER_C(Ch) __tolower_l ((Ch), _nl_C_locobj_ptr)
118 # define STRNCASECMP(S1, S2, N) \
119 __strncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
120 # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc)
121 #endif
124 /* Constants we need from float.h; select the set for the FLOAT precision. */
125 #define MANT_DIG PASTE(FLT,_MANT_DIG)
126 #define DIG PASTE(FLT,_DIG)
127 #define MAX_EXP PASTE(FLT,_MAX_EXP)
128 #define MIN_EXP PASTE(FLT,_MIN_EXP)
129 #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
130 #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
131 #define MAX_VALUE PASTE(FLT,_MAX)
132 #define MIN_VALUE PASTE(FLT,_MIN)
134 /* Extra macros required to get FLT expanded before the pasting. */
135 #define PASTE(a,b) PASTE1(a,b)
136 #define PASTE1(a,b) a##b
138 /* Function to construct a floating point number from an MP integer
139 containing the fraction bits, a base 2 exponent, and a sign flag. */
140 extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
142 /* Definitions according to limb size used. */
143 #if BITS_PER_MP_LIMB == 32
144 # define MAX_DIG_PER_LIMB 9
145 # define MAX_FAC_PER_LIMB 1000000000UL
146 #elif BITS_PER_MP_LIMB == 64
147 # define MAX_DIG_PER_LIMB 19
148 # define MAX_FAC_PER_LIMB 10000000000000000000ULL
149 #else
150 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
151 #endif
153 extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1];
155 #ifndef howmany
156 #define howmany(x,y) (((x)+((y)-1))/(y))
157 #endif
158 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
160 #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
162 #define RETURN(val,end) \
163 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
164 return val; } while (0)
166 /* Maximum size necessary for mpn integers to hold floating point
167 numbers. The largest number we need to hold is 10^n where 2^-n is
168 1/4 ulp of the smallest representable value (that is, n = MANT_DIG
169 - MIN_EXP + 2). Approximate using 10^3 < 2^10. */
170 #define MPNSIZE (howmany (1 + ((MANT_DIG - MIN_EXP + 2) * 10) / 3, \
171 BITS_PER_MP_LIMB) + 2)
172 /* Declare an mpn integer variable that big. */
173 #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
174 /* Copy an mpn integer value. */
175 #define MPN_ASSIGN(dst, src) \
176 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
179 /* Set errno and return an overflowing value with sign specified by
180 NEGATIVE. */
181 static FLOAT
182 overflow_value (int negative)
184 __set_errno (ERANGE);
185 #if FLT_EVAL_METHOD != 0
186 volatile
187 #endif
188 FLOAT result = (negative ? -MAX_VALUE : MAX_VALUE) * MAX_VALUE;
189 return result;
193 /* Set errno and return an underflowing value with sign specified by
194 NEGATIVE. */
195 static FLOAT
196 underflow_value (int negative)
198 __set_errno (ERANGE);
199 #if FLT_EVAL_METHOD != 0
200 volatile
201 #endif
202 FLOAT result = (negative ? -MIN_VALUE : MIN_VALUE) * MIN_VALUE;
203 return result;
207 /* Return a floating point number of the needed type according to the given
208 multi-precision number after possible rounding. */
209 static FLOAT
210 round_and_return (mp_limb_t *retval, intmax_t exponent, int negative,
211 mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
213 int mode = get_rounding_mode ();
215 if (exponent < MIN_EXP - 1)
217 if (exponent < MIN_EXP - 1 - MANT_DIG)
218 return underflow_value (negative);
220 mp_size_t shift = MIN_EXP - 1 - exponent;
221 bool is_tiny = true;
223 more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
224 if (shift == MANT_DIG)
225 /* This is a special case to handle the very seldom case where
226 the mantissa will be empty after the shift. */
228 int i;
230 round_limb = retval[RETURN_LIMB_SIZE - 1];
231 round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
232 for (i = 0; i < RETURN_LIMB_SIZE; ++i)
233 more_bits |= retval[i] != 0;
234 MPN_ZERO (retval, RETURN_LIMB_SIZE);
236 else if (shift >= BITS_PER_MP_LIMB)
238 int i;
240 round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
241 round_bit = (shift - 1) % BITS_PER_MP_LIMB;
242 for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
243 more_bits |= retval[i] != 0;
244 more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
245 != 0);
247 (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
248 RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
249 shift % BITS_PER_MP_LIMB);
250 MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
251 shift / BITS_PER_MP_LIMB);
253 else if (shift > 0)
255 if (TININESS_AFTER_ROUNDING && shift == 1)
257 /* Whether the result counts as tiny depends on whether,
258 after rounding to the normal precision, it still has
259 a subnormal exponent. */
260 mp_limb_t retval_normal[RETURN_LIMB_SIZE];
261 if (round_away (negative,
262 (retval[0] & 1) != 0,
263 (round_limb
264 & (((mp_limb_t) 1) << round_bit)) != 0,
265 (more_bits
266 || ((round_limb
267 & ((((mp_limb_t) 1) << round_bit) - 1))
268 != 0)),
269 mode))
271 mp_limb_t cy = __mpn_add_1 (retval_normal, retval,
272 RETURN_LIMB_SIZE, 1);
274 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
275 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
276 ((retval_normal[RETURN_LIMB_SIZE - 1]
277 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB)))
278 != 0)))
279 is_tiny = false;
282 round_limb = retval[0];
283 round_bit = shift - 1;
284 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
286 /* This is a hook for the m68k long double format, where the
287 exponent bias is the same for normalized and denormalized
288 numbers. */
289 #ifndef DENORM_EXP
290 # define DENORM_EXP (MIN_EXP - 2)
291 #endif
292 exponent = DENORM_EXP;
293 if (is_tiny
294 && ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
295 || more_bits
296 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
298 __set_errno (ERANGE);
299 volatile FLOAT force_underflow_exception = MIN_VALUE * MIN_VALUE;
300 (void) force_underflow_exception;
304 if (exponent > MAX_EXP)
305 goto overflow;
307 if (round_away (negative,
308 (retval[0] & 1) != 0,
309 (round_limb & (((mp_limb_t) 1) << round_bit)) != 0,
310 (more_bits
311 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0),
312 mode))
314 mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
316 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
317 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
318 (retval[RETURN_LIMB_SIZE - 1]
319 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
321 ++exponent;
322 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
323 retval[RETURN_LIMB_SIZE - 1]
324 |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
326 else if (exponent == DENORM_EXP
327 && (retval[RETURN_LIMB_SIZE - 1]
328 & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
329 != 0)
330 /* The number was denormalized but now normalized. */
331 exponent = MIN_EXP - 1;
334 if (exponent > MAX_EXP)
335 overflow:
336 return overflow_value (negative);
338 return MPN2FLOAT (retval, exponent, negative);
342 /* Read a multi-precision integer starting at STR with exactly DIGCNT digits
343 into N. Return the size of the number limbs in NSIZE at the first
344 character od the string that is not part of the integer as the function
345 value. If the EXPONENT is small enough to be taken as an additional
346 factor for the resulting number (see code) multiply by it. */
347 static const STRING_TYPE *
348 str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
349 intmax_t *exponent
350 #ifndef USE_WIDE_CHAR
351 , const char *decimal, size_t decimal_len, const char *thousands
352 #endif
356 /* Number of digits for actual limb. */
357 int cnt = 0;
358 mp_limb_t low = 0;
359 mp_limb_t start;
361 *nsize = 0;
362 assert (digcnt > 0);
365 if (cnt == MAX_DIG_PER_LIMB)
367 if (*nsize == 0)
369 n[0] = low;
370 *nsize = 1;
372 else
374 mp_limb_t cy;
375 cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
376 cy += __mpn_add_1 (n, n, *nsize, low);
377 if (cy != 0)
379 assert (*nsize < MPNSIZE);
380 n[*nsize] = cy;
381 ++(*nsize);
384 cnt = 0;
385 low = 0;
388 /* There might be thousands separators or radix characters in
389 the string. But these all can be ignored because we know the
390 format of the number is correct and we have an exact number
391 of characters to read. */
392 #ifdef USE_WIDE_CHAR
393 if (*str < L'0' || *str > L'9')
394 ++str;
395 #else
396 if (*str < '0' || *str > '9')
398 int inner = 0;
399 if (thousands != NULL && *str == *thousands
400 && ({ for (inner = 1; thousands[inner] != '\0'; ++inner)
401 if (thousands[inner] != str[inner])
402 break;
403 thousands[inner] == '\0'; }))
404 str += inner;
405 else
406 str += decimal_len;
408 #endif
409 low = low * 10 + *str++ - L_('0');
410 ++cnt;
412 while (--digcnt > 0);
414 if (*exponent > 0 && *exponent <= MAX_DIG_PER_LIMB - cnt)
416 low *= _tens_in_limb[*exponent];
417 start = _tens_in_limb[cnt + *exponent];
418 *exponent = 0;
420 else
421 start = _tens_in_limb[cnt];
423 if (*nsize == 0)
425 n[0] = low;
426 *nsize = 1;
428 else
430 mp_limb_t cy;
431 cy = __mpn_mul_1 (n, n, *nsize, start);
432 cy += __mpn_add_1 (n, n, *nsize, low);
433 if (cy != 0)
435 assert (*nsize < MPNSIZE);
436 n[(*nsize)++] = cy;
440 return str;
444 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
445 with the COUNT most significant bits of LIMB.
447 Implemented as a macro, so that __builtin_constant_p works even at -O0.
449 Tege doesn't like this macro so I have to write it here myself. :)
450 --drepper */
451 #define __mpn_lshift_1(ptr, size, count, limb) \
452 do \
454 mp_limb_t *__ptr = (ptr); \
455 if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \
457 mp_size_t i; \
458 for (i = (size) - 1; i > 0; --i) \
459 __ptr[i] = __ptr[i - 1]; \
460 __ptr[0] = (limb); \
462 else \
464 /* We assume count > 0 && count < BITS_PER_MP_LIMB here. */ \
465 unsigned int __count = (count); \
466 (void) __mpn_lshift (__ptr, __ptr, size, __count); \
467 __ptr[0] |= (limb) >> (BITS_PER_MP_LIMB - __count); \
470 while (0)
473 #define INTERNAL(x) INTERNAL1(x)
474 #define INTERNAL1(x) __##x##_internal
475 #ifndef ____STRTOF_INTERNAL
476 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
477 #endif
479 /* This file defines a function to check for correct grouping. */
480 #include "grouping.h"
483 /* Return a floating point number with the value of the given string NPTR.
484 Set *ENDPTR to the character after the last used one. If the number is
485 smaller than the smallest representable number, set `errno' to ERANGE and
486 return 0.0. If the number is too big to be represented, set `errno' to
487 ERANGE and return HUGE_VAL with the appropriate sign. */
488 FLOAT
489 ____STRTOF_INTERNAL (nptr, endptr, group, loc)
490 const STRING_TYPE *nptr;
491 STRING_TYPE **endptr;
492 int group;
493 __locale_t loc;
495 int negative; /* The sign of the number. */
496 MPN_VAR (num); /* MP representation of the number. */
497 intmax_t exponent; /* Exponent of the number. */
499 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
500 int base = 10;
502 /* When we have to compute fractional digits we form a fraction with a
503 second multi-precision number (and we sometimes need a second for
504 temporary results). */
505 MPN_VAR (den);
507 /* Representation for the return value. */
508 mp_limb_t retval[RETURN_LIMB_SIZE];
509 /* Number of bits currently in result value. */
510 int bits;
512 /* Running pointer after the last character processed in the string. */
513 const STRING_TYPE *cp, *tp;
514 /* Start of significant part of the number. */
515 const STRING_TYPE *startp, *start_of_digits;
516 /* Points at the character following the integer and fractional digits. */
517 const STRING_TYPE *expp;
518 /* Total number of digit and number of digits in integer part. */
519 size_t dig_no, int_no, lead_zero;
520 /* Contains the last character read. */
521 CHAR_TYPE c;
523 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
524 there. So define it ourselves if it remains undefined. */
525 #ifndef _WINT_T
526 typedef unsigned int wint_t;
527 #endif
528 /* The radix character of the current locale. */
529 #ifdef USE_WIDE_CHAR
530 wchar_t decimal;
531 #else
532 const char *decimal;
533 size_t decimal_len;
534 #endif
535 /* The thousands character of the current locale. */
536 #ifdef USE_WIDE_CHAR
537 wchar_t thousands = L'\0';
538 #else
539 const char *thousands = NULL;
540 #endif
541 /* The numeric grouping specification of the current locale,
542 in the format described in <locale.h>. */
543 const char *grouping;
544 /* Used in several places. */
545 int cnt;
547 struct __locale_data *current = loc->__locales[LC_NUMERIC];
549 if (__builtin_expect (group, 0))
551 grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
552 if (*grouping <= 0 || *grouping == CHAR_MAX)
553 grouping = NULL;
554 else
556 /* Figure out the thousands separator character. */
557 #ifdef USE_WIDE_CHAR
558 thousands = _NL_CURRENT_WORD (LC_NUMERIC,
559 _NL_NUMERIC_THOUSANDS_SEP_WC);
560 if (thousands == L'\0')
561 grouping = NULL;
562 #else
563 thousands = _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
564 if (*thousands == '\0')
566 thousands = NULL;
567 grouping = NULL;
569 #endif
572 else
573 grouping = NULL;
575 /* Find the locale's decimal point character. */
576 #ifdef USE_WIDE_CHAR
577 decimal = _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_DECIMAL_POINT_WC);
578 assert (decimal != L'\0');
579 # define decimal_len 1
580 #else
581 decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
582 decimal_len = strlen (decimal);
583 assert (decimal_len > 0);
584 #endif
586 /* Prepare number representation. */
587 exponent = 0;
588 negative = 0;
589 bits = 0;
591 /* Parse string to get maximal legal prefix. We need the number of
592 characters of the integer part, the fractional part and the exponent. */
593 cp = nptr - 1;
594 /* Ignore leading white space. */
596 c = *++cp;
597 while (ISSPACE (c));
599 /* Get sign of the result. */
600 if (c == L_('-'))
602 negative = 1;
603 c = *++cp;
605 else if (c == L_('+'))
606 c = *++cp;
608 /* Return 0.0 if no legal string is found.
609 No character is used even if a sign was found. */
610 #ifdef USE_WIDE_CHAR
611 if (c == (wint_t) decimal
612 && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9')
614 /* We accept it. This funny construct is here only to indent
615 the code correctly. */
617 #else
618 for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
619 if (cp[cnt] != decimal[cnt])
620 break;
621 if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9')
623 /* We accept it. This funny construct is here only to indent
624 the code correctly. */
626 #endif
627 else if (c < L_('0') || c > L_('9'))
629 /* Check for `INF' or `INFINITY'. */
630 CHAR_TYPE lowc = TOLOWER_C (c);
632 if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0)
634 /* Return +/- infinity. */
635 if (endptr != NULL)
636 *endptr = (STRING_TYPE *)
637 (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0
638 ? 8 : 3));
640 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
643 if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
645 /* Return NaN. */
646 FLOAT retval = NAN;
648 cp += 3;
650 /* Match `(n-char-sequence-digit)'. */
651 if (*cp == L_('('))
653 const STRING_TYPE *startp = cp;
655 ++cp;
656 while ((*cp >= L_('0') && *cp <= L_('9'))
657 || ({ CHAR_TYPE lo = TOLOWER (*cp);
658 lo >= L_('a') && lo <= L_('z'); })
659 || *cp == L_('_'));
661 if (*cp != L_(')'))
662 /* The closing brace is missing. Only match the NAN
663 part. */
664 cp = startp;
665 else
667 /* This is a system-dependent way to specify the
668 bitmask used for the NaN. We expect it to be
669 a number which is put in the mantissa of the
670 number. */
671 STRING_TYPE *endp;
672 unsigned long long int mant;
674 mant = STRTOULL (startp + 1, &endp, 0);
675 if (endp == cp)
676 SET_MANTISSA (retval, mant);
678 /* Consume the closing brace. */
679 ++cp;
683 if (endptr != NULL)
684 *endptr = (STRING_TYPE *) cp;
686 return retval;
689 /* It is really a text we do not recognize. */
690 RETURN (0.0, nptr);
693 /* First look whether we are faced with a hexadecimal number. */
694 if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
696 /* Okay, it is a hexa-decimal number. Remember this and skip
697 the characters. BTW: hexadecimal numbers must not be
698 grouped. */
699 base = 16;
700 cp += 2;
701 c = *cp;
702 grouping = NULL;
705 /* Record the start of the digits, in case we will check their grouping. */
706 start_of_digits = startp = cp;
708 /* Ignore leading zeroes. This helps us to avoid useless computations. */
709 #ifdef USE_WIDE_CHAR
710 while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands))
711 c = *++cp;
712 #else
713 if (__builtin_expect (thousands == NULL, 1))
714 while (c == '0')
715 c = *++cp;
716 else
718 /* We also have the multibyte thousands string. */
719 while (1)
721 if (c != '0')
723 for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
724 if (thousands[cnt] != cp[cnt])
725 break;
726 if (thousands[cnt] != '\0')
727 break;
728 cp += cnt - 1;
730 c = *++cp;
733 #endif
735 /* If no other digit but a '0' is found the result is 0.0.
736 Return current read pointer. */
737 CHAR_TYPE lowc = TOLOWER (c);
738 if (!((c >= L_('0') && c <= L_('9'))
739 || (base == 16 && lowc >= L_('a') && lowc <= L_('f'))
740 || (
741 #ifdef USE_WIDE_CHAR
742 c == (wint_t) decimal
743 #else
744 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
745 if (decimal[cnt] != cp[cnt])
746 break;
747 decimal[cnt] == '\0'; })
748 #endif
749 /* '0x.' alone is not a valid hexadecimal number.
750 '.' alone is not valid either, but that has been checked
751 already earlier. */
752 && (base != 16
753 || cp != start_of_digits
754 || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9'))
755 || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]);
756 lo >= L_('a') && lo <= L_('f'); })))
757 || (base == 16 && (cp != start_of_digits
758 && lowc == L_('p')))
759 || (base != 16 && lowc == L_('e'))))
761 #ifdef USE_WIDE_CHAR
762 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
763 grouping);
764 #else
765 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
766 grouping);
767 #endif
768 /* If TP is at the start of the digits, there was no correctly
769 grouped prefix of the string; so no number found. */
770 RETURN (negative ? -0.0 : 0.0,
771 tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
774 /* Remember first significant digit and read following characters until the
775 decimal point, exponent character or any non-FP number character. */
776 startp = cp;
777 dig_no = 0;
778 while (1)
780 if ((c >= L_('0') && c <= L_('9'))
781 || (base == 16
782 && ({ CHAR_TYPE lo = TOLOWER (c);
783 lo >= L_('a') && lo <= L_('f'); })))
784 ++dig_no;
785 else
787 #ifdef USE_WIDE_CHAR
788 if (__builtin_expect ((wint_t) thousands == L'\0', 1)
789 || c != (wint_t) thousands)
790 /* Not a digit or separator: end of the integer part. */
791 break;
792 #else
793 if (__builtin_expect (thousands == NULL, 1))
794 break;
795 else
797 for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
798 if (thousands[cnt] != cp[cnt])
799 break;
800 if (thousands[cnt] != '\0')
801 break;
802 cp += cnt - 1;
804 #endif
806 c = *++cp;
809 if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits)
811 /* Check the grouping of the digits. */
812 #ifdef USE_WIDE_CHAR
813 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
814 grouping);
815 #else
816 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
817 grouping);
818 #endif
819 if (cp != tp)
821 /* Less than the entire string was correctly grouped. */
823 if (tp == start_of_digits)
824 /* No valid group of numbers at all: no valid number. */
825 RETURN (0.0, nptr);
827 if (tp < startp)
828 /* The number is validly grouped, but consists
829 only of zeroes. The whole value is zero. */
830 RETURN (negative ? -0.0 : 0.0, tp);
832 /* Recompute DIG_NO so we won't read more digits than
833 are properly grouped. */
834 cp = tp;
835 dig_no = 0;
836 for (tp = startp; tp < cp; ++tp)
837 if (*tp >= L_('0') && *tp <= L_('9'))
838 ++dig_no;
840 int_no = dig_no;
841 lead_zero = 0;
843 goto number_parsed;
847 /* We have the number of digits in the integer part. Whether these
848 are all or any is really a fractional digit will be decided
849 later. */
850 int_no = dig_no;
851 lead_zero = int_no == 0 ? (size_t) -1 : 0;
853 /* Read the fractional digits. A special case are the 'american
854 style' numbers like `16.' i.e. with decimal point but without
855 trailing digits. */
856 if (
857 #ifdef USE_WIDE_CHAR
858 c == (wint_t) decimal
859 #else
860 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
861 if (decimal[cnt] != cp[cnt])
862 break;
863 decimal[cnt] == '\0'; })
864 #endif
867 cp += decimal_len;
868 c = *cp;
869 while ((c >= L_('0') && c <= L_('9')) ||
870 (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c);
871 lo >= L_('a') && lo <= L_('f'); })))
873 if (c != L_('0') && lead_zero == (size_t) -1)
874 lead_zero = dig_no - int_no;
875 ++dig_no;
876 c = *++cp;
879 assert (dig_no <= (uintmax_t) INTMAX_MAX);
881 /* Remember start of exponent (if any). */
882 expp = cp;
884 /* Read exponent. */
885 lowc = TOLOWER (c);
886 if ((base == 16 && lowc == L_('p'))
887 || (base != 16 && lowc == L_('e')))
889 int exp_negative = 0;
891 c = *++cp;
892 if (c == L_('-'))
894 exp_negative = 1;
895 c = *++cp;
897 else if (c == L_('+'))
898 c = *++cp;
900 if (c >= L_('0') && c <= L_('9'))
902 intmax_t exp_limit;
904 /* Get the exponent limit. */
905 if (base == 16)
907 if (exp_negative)
909 assert (int_no <= (uintmax_t) (INTMAX_MAX
910 + MIN_EXP - MANT_DIG) / 4);
911 exp_limit = -MIN_EXP + MANT_DIG + 4 * (intmax_t) int_no;
913 else
915 if (int_no)
917 assert (lead_zero == 0
918 && int_no <= (uintmax_t) INTMAX_MAX / 4);
919 exp_limit = MAX_EXP - 4 * (intmax_t) int_no + 3;
921 else if (lead_zero == (size_t) -1)
923 /* The number is zero and this limit is
924 arbitrary. */
925 exp_limit = MAX_EXP + 3;
927 else
929 assert (lead_zero
930 <= (uintmax_t) (INTMAX_MAX - MAX_EXP - 3) / 4);
931 exp_limit = (MAX_EXP
932 + 4 * (intmax_t) lead_zero
933 + 3);
937 else
939 if (exp_negative)
941 assert (int_no
942 <= (uintmax_t) (INTMAX_MAX + MIN_10_EXP - MANT_DIG));
943 exp_limit = -MIN_10_EXP + MANT_DIG + (intmax_t) int_no;
945 else
947 if (int_no)
949 assert (lead_zero == 0
950 && int_no <= (uintmax_t) INTMAX_MAX);
951 exp_limit = MAX_10_EXP - (intmax_t) int_no + 1;
953 else if (lead_zero == (size_t) -1)
955 /* The number is zero and this limit is
956 arbitrary. */
957 exp_limit = MAX_10_EXP + 1;
959 else
961 assert (lead_zero
962 <= (uintmax_t) (INTMAX_MAX - MAX_10_EXP - 1));
963 exp_limit = MAX_10_EXP + (intmax_t) lead_zero + 1;
968 if (exp_limit < 0)
969 exp_limit = 0;
973 if (__builtin_expect ((exponent > exp_limit / 10
974 || (exponent == exp_limit / 10
975 && c - L_('0') > exp_limit % 10)), 0))
976 /* The exponent is too large/small to represent a valid
977 number. */
979 FLOAT result;
981 /* We have to take care for special situation: a joker
982 might have written "0.0e100000" which is in fact
983 zero. */
984 if (lead_zero == (size_t) -1)
985 result = negative ? -0.0 : 0.0;
986 else
988 /* Overflow or underflow. */
989 result = (exp_negative
990 ? underflow_value (negative)
991 : overflow_value (negative));
994 /* Accept all following digits as part of the exponent. */
996 ++cp;
997 while (*cp >= L_('0') && *cp <= L_('9'));
999 RETURN (result, cp);
1000 /* NOTREACHED */
1003 exponent *= 10;
1004 exponent += c - L_('0');
1006 c = *++cp;
1008 while (c >= L_('0') && c <= L_('9'));
1010 if (exp_negative)
1011 exponent = -exponent;
1013 else
1014 cp = expp;
1017 /* We don't want to have to work with trailing zeroes after the radix. */
1018 if (dig_no > int_no)
1020 while (expp[-1] == L_('0'))
1022 --expp;
1023 --dig_no;
1025 assert (dig_no >= int_no);
1028 if (dig_no == int_no && dig_no > 0 && exponent < 0)
1031 while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1])))
1032 --expp;
1034 if (expp[-1] != L_('0'))
1035 break;
1037 --expp;
1038 --dig_no;
1039 --int_no;
1040 exponent += base == 16 ? 4 : 1;
1042 while (dig_no > 0 && exponent < 0);
1044 number_parsed:
1046 /* The whole string is parsed. Store the address of the next character. */
1047 if (endptr)
1048 *endptr = (STRING_TYPE *) cp;
1050 if (dig_no == 0)
1051 return negative ? -0.0 : 0.0;
1053 if (lead_zero)
1055 /* Find the decimal point */
1056 #ifdef USE_WIDE_CHAR
1057 while (*startp != decimal)
1058 ++startp;
1059 #else
1060 while (1)
1062 if (*startp == decimal[0])
1064 for (cnt = 1; decimal[cnt] != '\0'; ++cnt)
1065 if (decimal[cnt] != startp[cnt])
1066 break;
1067 if (decimal[cnt] == '\0')
1068 break;
1070 ++startp;
1072 #endif
1073 startp += lead_zero + decimal_len;
1074 assert (lead_zero <= (base == 16
1075 ? (uintmax_t) INTMAX_MAX / 4
1076 : (uintmax_t) INTMAX_MAX));
1077 assert (lead_zero <= (base == 16
1078 ? ((uintmax_t) exponent
1079 - (uintmax_t) INTMAX_MIN) / 4
1080 : ((uintmax_t) exponent - (uintmax_t) INTMAX_MIN)));
1081 exponent -= base == 16 ? 4 * (intmax_t) lead_zero : (intmax_t) lead_zero;
1082 dig_no -= lead_zero;
1085 /* If the BASE is 16 we can use a simpler algorithm. */
1086 if (base == 16)
1088 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
1089 4, 4, 4, 4, 4, 4, 4, 4 };
1090 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
1091 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
1092 mp_limb_t val;
1094 while (!ISXDIGIT (*startp))
1095 ++startp;
1096 while (*startp == L_('0'))
1097 ++startp;
1098 if (ISDIGIT (*startp))
1099 val = *startp++ - L_('0');
1100 else
1101 val = 10 + TOLOWER (*startp++) - L_('a');
1102 bits = nbits[val];
1103 /* We cannot have a leading zero. */
1104 assert (bits != 0);
1106 if (pos + 1 >= 4 || pos + 1 >= bits)
1108 /* We don't have to care for wrapping. This is the normal
1109 case so we add the first clause in the `if' expression as
1110 an optimization. It is a compile-time constant and so does
1111 not cost anything. */
1112 retval[idx] = val << (pos - bits + 1);
1113 pos -= bits;
1115 else
1117 retval[idx--] = val >> (bits - pos - 1);
1118 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
1119 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
1122 /* Adjust the exponent for the bits we are shifting in. */
1123 assert (int_no <= (uintmax_t) (exponent < 0
1124 ? (INTMAX_MAX - bits + 1) / 4
1125 : (INTMAX_MAX - exponent - bits + 1) / 4));
1126 exponent += bits - 1 + ((intmax_t) int_no - 1) * 4;
1128 while (--dig_no > 0 && idx >= 0)
1130 if (!ISXDIGIT (*startp))
1131 startp += decimal_len;
1132 if (ISDIGIT (*startp))
1133 val = *startp++ - L_('0');
1134 else
1135 val = 10 + TOLOWER (*startp++) - L_('a');
1137 if (pos + 1 >= 4)
1139 retval[idx] |= val << (pos - 4 + 1);
1140 pos -= 4;
1142 else
1144 retval[idx--] |= val >> (4 - pos - 1);
1145 val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
1146 if (idx < 0)
1148 int rest_nonzero = 0;
1149 while (--dig_no > 0)
1151 if (*startp != L_('0'))
1153 rest_nonzero = 1;
1154 break;
1156 startp++;
1158 return round_and_return (retval, exponent, negative, val,
1159 BITS_PER_MP_LIMB - 1, rest_nonzero);
1162 retval[idx] = val;
1163 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
1167 /* We ran out of digits. */
1168 MPN_ZERO (retval, idx);
1170 return round_and_return (retval, exponent, negative, 0, 0, 0);
1173 /* Now we have the number of digits in total and the integer digits as well
1174 as the exponent and its sign. We can decide whether the read digits are
1175 really integer digits or belong to the fractional part; i.e. we normalize
1176 123e-2 to 1.23. */
1178 register intmax_t incr = (exponent < 0
1179 ? MAX (-(intmax_t) int_no, exponent)
1180 : MIN ((intmax_t) dig_no - (intmax_t) int_no,
1181 exponent));
1182 int_no += incr;
1183 exponent -= incr;
1186 if (__builtin_expect (exponent > MAX_10_EXP + 1 - (intmax_t) int_no, 0))
1187 return overflow_value (negative);
1189 if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0))
1190 return underflow_value (negative);
1192 if (int_no > 0)
1194 /* Read the integer part as a multi-precision number to NUM. */
1195 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent
1196 #ifndef USE_WIDE_CHAR
1197 , decimal, decimal_len, thousands
1198 #endif
1201 if (exponent > 0)
1203 /* We now multiply the gained number by the given power of ten. */
1204 mp_limb_t *psrc = num;
1205 mp_limb_t *pdest = den;
1206 int expbit = 1;
1207 const struct mp_power *ttab = &_fpioconst_pow10[0];
1211 if ((exponent & expbit) != 0)
1213 size_t size = ttab->arraysize - _FPIO_CONST_OFFSET;
1214 mp_limb_t cy;
1215 exponent ^= expbit;
1217 /* FIXME: not the whole multiplication has to be
1218 done. If we have the needed number of bits we
1219 only need the information whether more non-zero
1220 bits follow. */
1221 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
1222 cy = __mpn_mul (pdest, psrc, numsize,
1223 &__tens[ttab->arrayoff
1224 + _FPIO_CONST_OFFSET],
1225 size);
1226 else
1227 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
1228 + _FPIO_CONST_OFFSET],
1229 size, psrc, numsize);
1230 numsize += size;
1231 if (cy == 0)
1232 --numsize;
1233 (void) SWAP (psrc, pdest);
1235 expbit <<= 1;
1236 ++ttab;
1238 while (exponent != 0);
1240 if (psrc == den)
1241 memcpy (num, den, numsize * sizeof (mp_limb_t));
1244 /* Determine how many bits of the result we already have. */
1245 count_leading_zeros (bits, num[numsize - 1]);
1246 bits = numsize * BITS_PER_MP_LIMB - bits;
1248 /* Now we know the exponent of the number in base two.
1249 Check it against the maximum possible exponent. */
1250 if (__builtin_expect (bits > MAX_EXP, 0))
1251 return overflow_value (negative);
1253 /* We have already the first BITS bits of the result. Together with
1254 the information whether more non-zero bits follow this is enough
1255 to determine the result. */
1256 if (bits > MANT_DIG)
1258 int i;
1259 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
1260 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
1261 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
1262 : least_idx;
1263 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
1264 : least_bit - 1;
1266 if (least_bit == 0)
1267 memcpy (retval, &num[least_idx],
1268 RETURN_LIMB_SIZE * sizeof (mp_limb_t));
1269 else
1271 for (i = least_idx; i < numsize - 1; ++i)
1272 retval[i - least_idx] = (num[i] >> least_bit)
1273 | (num[i + 1]
1274 << (BITS_PER_MP_LIMB - least_bit));
1275 if (i - least_idx < RETURN_LIMB_SIZE)
1276 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
1279 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1280 for (i = 0; num[i] == 0; ++i)
1283 return round_and_return (retval, bits - 1, negative,
1284 num[round_idx], round_bit,
1285 int_no < dig_no || i < round_idx);
1286 /* NOTREACHED */
1288 else if (dig_no == int_no)
1290 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
1291 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
1293 if (target_bit == is_bit)
1295 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
1296 numsize * sizeof (mp_limb_t));
1297 /* FIXME: the following loop can be avoided if we assume a
1298 maximal MANT_DIG value. */
1299 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1301 else if (target_bit > is_bit)
1303 (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
1304 num, numsize, target_bit - is_bit);
1305 /* FIXME: the following loop can be avoided if we assume a
1306 maximal MANT_DIG value. */
1307 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1309 else
1311 mp_limb_t cy;
1312 assert (numsize < RETURN_LIMB_SIZE);
1314 cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
1315 num, numsize, is_bit - target_bit);
1316 retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
1317 /* FIXME: the following loop can be avoided if we assume a
1318 maximal MANT_DIG value. */
1319 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
1322 return round_and_return (retval, bits - 1, negative, 0, 0, 0);
1323 /* NOTREACHED */
1326 /* Store the bits we already have. */
1327 memcpy (retval, num, numsize * sizeof (mp_limb_t));
1328 #if RETURN_LIMB_SIZE > 1
1329 if (numsize < RETURN_LIMB_SIZE)
1330 # if RETURN_LIMB_SIZE == 2
1331 retval[numsize] = 0;
1332 # else
1333 MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize);
1334 # endif
1335 #endif
1338 /* We have to compute at least some of the fractional digits. */
1340 /* We construct a fraction and the result of the division gives us
1341 the needed digits. The denominator is 1.0 multiplied by the
1342 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1343 123e-6 gives 123 / 1000000. */
1345 int expbit;
1346 int neg_exp;
1347 int more_bits;
1348 int need_frac_digits;
1349 mp_limb_t cy;
1350 mp_limb_t *psrc = den;
1351 mp_limb_t *pdest = num;
1352 const struct mp_power *ttab = &_fpioconst_pow10[0];
1354 assert (dig_no > int_no
1355 && exponent <= 0
1356 && exponent >= MIN_10_EXP - (DIG + 1));
1358 /* We need to compute MANT_DIG - BITS fractional bits that lie
1359 within the mantissa of the result, the following bit for
1360 rounding, and to know whether any subsequent bit is 0.
1361 Computing a bit with value 2^-n means looking at n digits after
1362 the decimal point. */
1363 if (bits > 0)
1365 /* The bits required are those immediately after the point. */
1366 assert (int_no > 0 && exponent == 0);
1367 need_frac_digits = 1 + MANT_DIG - bits;
1369 else
1371 /* The number is in the form .123eEXPONENT. */
1372 assert (int_no == 0 && *startp != L_('0'));
1373 /* The number is at least 10^(EXPONENT-1), and 10^3 <
1374 2^10. */
1375 int neg_exp_2 = ((1 - exponent) * 10) / 3 + 1;
1376 /* The number is at least 2^-NEG_EXP_2. We need up to
1377 MANT_DIG bits following that bit. */
1378 need_frac_digits = neg_exp_2 + MANT_DIG;
1379 /* However, we never need bits beyond 1/4 ulp of the smallest
1380 representable value. (That 1/4 ulp bit is only needed to
1381 determine tinyness on machines where tinyness is determined
1382 after rounding.) */
1383 if (need_frac_digits > MANT_DIG - MIN_EXP + 2)
1384 need_frac_digits = MANT_DIG - MIN_EXP + 2;
1385 /* At this point, NEED_FRAC_DIGITS is the total number of
1386 digits needed after the point, but some of those may be
1387 leading 0s. */
1388 need_frac_digits += exponent;
1389 /* Any cases underflowing enough that none of the fractional
1390 digits are needed should have been caught earlier (such
1391 cases are on the order of 10^-n or smaller where 2^-n is
1392 the least subnormal). */
1393 assert (need_frac_digits > 0);
1396 if (need_frac_digits > (intmax_t) dig_no - (intmax_t) int_no)
1397 need_frac_digits = (intmax_t) dig_no - (intmax_t) int_no;
1399 if ((intmax_t) dig_no > (intmax_t) int_no + need_frac_digits)
1401 dig_no = int_no + need_frac_digits;
1402 more_bits = 1;
1404 else
1405 more_bits = 0;
1407 neg_exp = (intmax_t) dig_no - (intmax_t) int_no - exponent;
1409 /* Construct the denominator. */
1410 densize = 0;
1411 expbit = 1;
1414 if ((neg_exp & expbit) != 0)
1416 mp_limb_t cy;
1417 neg_exp ^= expbit;
1419 if (densize == 0)
1421 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
1422 memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET],
1423 densize * sizeof (mp_limb_t));
1425 else
1427 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
1428 + _FPIO_CONST_OFFSET],
1429 ttab->arraysize - _FPIO_CONST_OFFSET,
1430 psrc, densize);
1431 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
1432 if (cy == 0)
1433 --densize;
1434 (void) SWAP (psrc, pdest);
1437 expbit <<= 1;
1438 ++ttab;
1440 while (neg_exp != 0);
1442 if (psrc == num)
1443 memcpy (den, num, densize * sizeof (mp_limb_t));
1445 /* Read the fractional digits from the string. */
1446 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent
1447 #ifndef USE_WIDE_CHAR
1448 , decimal, decimal_len, thousands
1449 #endif
1452 /* We now have to shift both numbers so that the highest bit in the
1453 denominator is set. In the same process we copy the numerator to
1454 a high place in the array so that the division constructs the wanted
1455 digits. This is done by a "quasi fix point" number representation.
1457 num: ddddddddddd . 0000000000000000000000
1458 |--- m ---|
1459 den: ddddddddddd n >= m
1460 |--- n ---|
1463 count_leading_zeros (cnt, den[densize - 1]);
1465 if (cnt > 0)
1467 /* Don't call `mpn_shift' with a count of zero since the specification
1468 does not allow this. */
1469 (void) __mpn_lshift (den, den, densize, cnt);
1470 cy = __mpn_lshift (num, num, numsize, cnt);
1471 if (cy != 0)
1472 num[numsize++] = cy;
1475 /* Now we are ready for the division. But it is not necessary to
1476 do a full multi-precision division because we only need a small
1477 number of bits for the result. So we do not use __mpn_divmod
1478 here but instead do the division here by hand and stop whenever
1479 the needed number of bits is reached. The code itself comes
1480 from the GNU MP Library by Torbj\"orn Granlund. */
1482 exponent = bits;
1484 switch (densize)
1486 case 1:
1488 mp_limb_t d, n, quot;
1489 int used = 0;
1491 n = num[0];
1492 d = den[0];
1493 assert (numsize == 1 && n < d);
1497 udiv_qrnnd (quot, n, n, 0, d);
1499 #define got_limb \
1500 if (bits == 0) \
1502 register int cnt; \
1503 if (quot == 0) \
1504 cnt = BITS_PER_MP_LIMB; \
1505 else \
1506 count_leading_zeros (cnt, quot); \
1507 exponent -= cnt; \
1508 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1510 used = MANT_DIG + cnt; \
1511 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1512 bits = MANT_DIG + 1; \
1514 else \
1516 /* Note that we only clear the second element. */ \
1517 /* The conditional is determined at compile time. */ \
1518 if (RETURN_LIMB_SIZE > 1) \
1519 retval[1] = 0; \
1520 retval[0] = quot; \
1521 bits = -cnt; \
1524 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1525 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1526 quot); \
1527 else \
1529 used = MANT_DIG - bits; \
1530 if (used > 0) \
1531 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1533 bits += BITS_PER_MP_LIMB
1535 got_limb;
1537 while (bits <= MANT_DIG);
1539 return round_and_return (retval, exponent - 1, negative,
1540 quot, BITS_PER_MP_LIMB - 1 - used,
1541 more_bits || n != 0);
1543 case 2:
1545 mp_limb_t d0, d1, n0, n1;
1546 mp_limb_t quot = 0;
1547 int used = 0;
1549 d0 = den[0];
1550 d1 = den[1];
1552 if (numsize < densize)
1554 if (num[0] >= d1)
1556 /* The numerator of the number occupies fewer bits than
1557 the denominator but the one limb is bigger than the
1558 high limb of the numerator. */
1559 n1 = 0;
1560 n0 = num[0];
1562 else
1564 if (bits <= 0)
1565 exponent -= BITS_PER_MP_LIMB;
1566 else
1568 if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
1569 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1570 BITS_PER_MP_LIMB, 0);
1571 else
1573 used = MANT_DIG - bits;
1574 if (used > 0)
1575 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1577 bits += BITS_PER_MP_LIMB;
1579 n1 = num[0];
1580 n0 = 0;
1583 else
1585 n1 = num[1];
1586 n0 = num[0];
1589 while (bits <= MANT_DIG)
1591 mp_limb_t r;
1593 if (n1 == d1)
1595 /* QUOT should be either 111..111 or 111..110. We need
1596 special treatment of this rare case as normal division
1597 would give overflow. */
1598 quot = ~(mp_limb_t) 0;
1600 r = n0 + d1;
1601 if (r < d1) /* Carry in the addition? */
1603 add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1604 goto have_quot;
1606 n1 = d0 - (d0 != 0);
1607 n0 = -d0;
1609 else
1611 udiv_qrnnd (quot, r, n1, n0, d1);
1612 umul_ppmm (n1, n0, d0, quot);
1615 q_test:
1616 if (n1 > r || (n1 == r && n0 > 0))
1618 /* The estimated QUOT was too large. */
1619 --quot;
1621 sub_ddmmss (n1, n0, n1, n0, 0, d0);
1622 r += d1;
1623 if (r >= d1) /* If not carry, test QUOT again. */
1624 goto q_test;
1626 sub_ddmmss (n1, n0, r, 0, n1, n0);
1628 have_quot:
1629 got_limb;
1632 return round_and_return (retval, exponent - 1, negative,
1633 quot, BITS_PER_MP_LIMB - 1 - used,
1634 more_bits || n1 != 0 || n0 != 0);
1636 default:
1638 int i;
1639 mp_limb_t cy, dX, d1, n0, n1;
1640 mp_limb_t quot = 0;
1641 int used = 0;
1643 dX = den[densize - 1];
1644 d1 = den[densize - 2];
1646 /* The division does not work if the upper limb of the two-limb
1647 numerator is greater than the denominator. */
1648 if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1649 num[numsize++] = 0;
1651 if (numsize < densize)
1653 mp_size_t empty = densize - numsize;
1654 register int i;
1656 if (bits <= 0)
1657 exponent -= empty * BITS_PER_MP_LIMB;
1658 else
1660 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1662 /* We make a difference here because the compiler
1663 cannot optimize the `else' case that good and
1664 this reflects all currently used FLOAT types
1665 and GMP implementations. */
1666 #if RETURN_LIMB_SIZE <= 2
1667 assert (empty == 1);
1668 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1669 BITS_PER_MP_LIMB, 0);
1670 #else
1671 for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i)
1672 retval[i] = retval[i - empty];
1673 while (i >= 0)
1674 retval[i--] = 0;
1675 #endif
1677 else
1679 used = MANT_DIG - bits;
1680 if (used >= BITS_PER_MP_LIMB)
1682 register int i;
1683 (void) __mpn_lshift (&retval[used
1684 / BITS_PER_MP_LIMB],
1685 retval,
1686 (RETURN_LIMB_SIZE
1687 - used / BITS_PER_MP_LIMB),
1688 used % BITS_PER_MP_LIMB);
1689 for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i)
1690 retval[i] = 0;
1692 else if (used > 0)
1693 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1695 bits += empty * BITS_PER_MP_LIMB;
1697 for (i = numsize; i > 0; --i)
1698 num[i + empty] = num[i - 1];
1699 MPN_ZERO (num, empty + 1);
1701 else
1703 int i;
1704 assert (numsize == densize);
1705 for (i = numsize; i > 0; --i)
1706 num[i] = num[i - 1];
1707 num[0] = 0;
1710 den[densize] = 0;
1711 n0 = num[densize];
1713 while (bits <= MANT_DIG)
1715 if (n0 == dX)
1716 /* This might over-estimate QUOT, but it's probably not
1717 worth the extra code here to find out. */
1718 quot = ~(mp_limb_t) 0;
1719 else
1721 mp_limb_t r;
1723 udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1724 umul_ppmm (n1, n0, d1, quot);
1726 while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1728 --quot;
1729 r += dX;
1730 if (r < dX) /* I.e. "carry in previous addition?" */
1731 break;
1732 n1 -= n0 < d1;
1733 n0 -= d1;
1737 /* Possible optimization: We already have (q * n0) and (1 * n1)
1738 after the calculation of QUOT. Taking advantage of this, we
1739 could make this loop make two iterations less. */
1741 cy = __mpn_submul_1 (num, den, densize + 1, quot);
1743 if (num[densize] != cy)
1745 cy = __mpn_add_n (num, num, den, densize);
1746 assert (cy != 0);
1747 --quot;
1749 n0 = num[densize] = num[densize - 1];
1750 for (i = densize - 1; i > 0; --i)
1751 num[i] = num[i - 1];
1752 num[0] = 0;
1754 got_limb;
1757 for (i = densize; num[i] == 0 && i >= 0; --i)
1759 return round_and_return (retval, exponent - 1, negative,
1760 quot, BITS_PER_MP_LIMB - 1 - used,
1761 more_bits || i >= 0);
1766 /* NOTREACHED */
1768 #if defined _LIBC && !defined USE_WIDE_CHAR
1769 libc_hidden_def (____STRTOF_INTERNAL)
1770 #endif
1772 /* External user entry point. */
1774 FLOAT
1775 #ifdef weak_function
1776 weak_function
1777 #endif
1778 __STRTOF (nptr, endptr, loc)
1779 const STRING_TYPE *nptr;
1780 STRING_TYPE **endptr;
1781 __locale_t loc;
1783 return ____STRTOF_INTERNAL (nptr, endptr, 0, loc);
1785 #if defined _LIBC
1786 libc_hidden_def (__STRTOF)
1787 libc_hidden_ver (__STRTOF, STRTOF)
1788 #endif
1789 weak_alias (__STRTOF, STRTOF)
1791 #ifdef LONG_DOUBLE_COMPAT
1792 # if LONG_DOUBLE_COMPAT(libc, GLIBC_2_1)
1793 # ifdef USE_WIDE_CHAR
1794 compat_symbol (libc, __wcstod_l, __wcstold_l, GLIBC_2_1);
1795 # else
1796 compat_symbol (libc, __strtod_l, __strtold_l, GLIBC_2_1);
1797 # endif
1798 # endif
1799 # if LONG_DOUBLE_COMPAT(libc, GLIBC_2_3)
1800 # ifdef USE_WIDE_CHAR
1801 compat_symbol (libc, wcstod_l, wcstold_l, GLIBC_2_3);
1802 # else
1803 compat_symbol (libc, strtod_l, strtold_l, GLIBC_2_3);
1804 # endif
1805 # endif
1806 #endif