1/* Convert string representing a number to float value, using given locale.
2 Copyright (C) 1997-2017 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
5
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.
10
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.
15
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/>. */
19
20#include <locale.h>
21
22extern double ____strtod_l_internal (const char *, char **, int, locale_t);
23
24/* Configuration part. These macros are defined by `strtold.c',
25 `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the
26 `long double' and `float' versions of the reader. */
27#ifndef FLOAT
28# include <math_ldbl_opt.h>
29# define FLOAT double
30# define FLT DBL
31# ifdef USE_WIDE_CHAR
32# define STRTOF wcstod_l
33# define __STRTOF __wcstod_l
34# define STRTOF_NAN __wcstod_nan
35# else
36# define STRTOF strtod_l
37# define __STRTOF __strtod_l
38# define STRTOF_NAN __strtod_nan
39# endif
40# define MPN2FLOAT __mpn_construct_double
41# define FLOAT_HUGE_VAL HUGE_VAL
42#endif
43/* End of configuration part. */
44
45#include <ctype.h>
46#include <errno.h>
47#include <float.h>
48#include "../locale/localeinfo.h"
49#include <math.h>
50#include <math_private.h>
51#include <stdlib.h>
52#include <string.h>
53#include <stdint.h>
54#include <rounding-mode.h>
55#include <tininess.h>
56
57/* The gmp headers need some configuration frobs. */
58#define HAVE_ALLOCA 1
59
60/* Include gmp-mparam.h first, such that definitions of _SHORT_LIMB
61 and _LONG_LONG_LIMB in it can take effect into gmp.h. */
62#include <gmp-mparam.h>
63#include <gmp.h>
64#include "gmp-impl.h"
65#include "longlong.h"
66#include "fpioconst.h"
67
68#include <assert.h>
69
70
71/* We use this code for the extended locale handling where the
72 function gets as an additional argument the locale which has to be
73 used. To access the values we have to redefine the _NL_CURRENT and
74 _NL_CURRENT_WORD macros. */
75#undef _NL_CURRENT
76#define _NL_CURRENT(category, item) \
77 (current->values[_NL_ITEM_INDEX (item)].string)
78#undef _NL_CURRENT_WORD
79#define _NL_CURRENT_WORD(category, item) \
80 ((uint32_t) current->values[_NL_ITEM_INDEX (item)].word)
81
82#if defined _LIBC || defined HAVE_WCHAR_H
83# include <wchar.h>
84#endif
85
86#ifdef USE_WIDE_CHAR
87# include <wctype.h>
88# define STRING_TYPE wchar_t
89# define CHAR_TYPE wint_t
90# define L_(Ch) L##Ch
91# define ISSPACE(Ch) __iswspace_l ((Ch), loc)
92# define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
93# define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
94# define TOLOWER(Ch) __towlower_l ((Ch), loc)
95# define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
96# define STRNCASECMP(S1, S2, N) \
97 __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
98#else
99# define STRING_TYPE char
100# define CHAR_TYPE char
101# define L_(Ch) Ch
102# define ISSPACE(Ch) __isspace_l ((Ch), loc)
103# define ISDIGIT(Ch) __isdigit_l ((Ch), loc)
104# define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc)
105# define TOLOWER(Ch) __tolower_l ((Ch), loc)
106# define TOLOWER_C(Ch) __tolower_l ((Ch), _nl_C_locobj_ptr)
107# define STRNCASECMP(S1, S2, N) \
108 __strncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
109#endif
110
111
112/* Constants we need from float.h; select the set for the FLOAT precision. */
113#define MANT_DIG PASTE(FLT,_MANT_DIG)
114#define DIG PASTE(FLT,_DIG)
115#define MAX_EXP PASTE(FLT,_MAX_EXP)
116#define MIN_EXP PASTE(FLT,_MIN_EXP)
117#define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
118#define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
119#define MAX_VALUE PASTE(FLT,_MAX)
120#define MIN_VALUE PASTE(FLT,_MIN)
121
122/* Extra macros required to get FLT expanded before the pasting. */
123#define PASTE(a,b) PASTE1(a,b)
124#define PASTE1(a,b) a##b
125
126/* Function to construct a floating point number from an MP integer
127 containing the fraction bits, a base 2 exponent, and a sign flag. */
128extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
129
130/* Definitions according to limb size used. */
131#if BITS_PER_MP_LIMB == 32
132# define MAX_DIG_PER_LIMB 9
133# define MAX_FAC_PER_LIMB 1000000000UL
134#elif BITS_PER_MP_LIMB == 64
135# define MAX_DIG_PER_LIMB 19
136# define MAX_FAC_PER_LIMB 10000000000000000000ULL
137#else
138# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
139#endif
140
141extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1];
142
143#ifndef howmany
144#define howmany(x,y) (((x)+((y)-1))/(y))
145#endif
146#define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
147
148#define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
149
150#define RETURN(val,end) \
151 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
152 return val; } while (0)
153
154/* Maximum size necessary for mpn integers to hold floating point
155 numbers. The largest number we need to hold is 10^n where 2^-n is
156 1/4 ulp of the smallest representable value (that is, n = MANT_DIG
157 - MIN_EXP + 2). Approximate using 10^3 < 2^10. */
158#define MPNSIZE (howmany (1 + ((MANT_DIG - MIN_EXP + 2) * 10) / 3, \
159 BITS_PER_MP_LIMB) + 2)
160/* Declare an mpn integer variable that big. */
161#define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
162/* Copy an mpn integer value. */
163#define MPN_ASSIGN(dst, src) \
164 memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
165
166
167/* Set errno and return an overflowing value with sign specified by
168 NEGATIVE. */
169static FLOAT
170overflow_value (int negative)
171{
172 __set_errno (ERANGE);
173 FLOAT result = math_narrow_eval ((negative ? -MAX_VALUE : MAX_VALUE)
174 * MAX_VALUE);
175 return result;
176}
177
178
179/* Set errno and return an underflowing value with sign specified by
180 NEGATIVE. */
181static FLOAT
182underflow_value (int negative)
183{
184 __set_errno (ERANGE);
185 FLOAT result = math_narrow_eval ((negative ? -MIN_VALUE : MIN_VALUE)
186 * MIN_VALUE);
187 return result;
188}
189
190
191/* Return a floating point number of the needed type according to the given
192 multi-precision number after possible rounding. */
193static FLOAT
194round_and_return (mp_limb_t *retval, intmax_t exponent, int negative,
195 mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
196{
197 int mode = get_rounding_mode ();
198
199 if (exponent < MIN_EXP - 1)
200 {
201 if (exponent < MIN_EXP - 1 - MANT_DIG)
202 return underflow_value (negative);
203
204 mp_size_t shift = MIN_EXP - 1 - exponent;
205 bool is_tiny = true;
206
207 more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
208 if (shift == MANT_DIG)
209 /* This is a special case to handle the very seldom case where
210 the mantissa will be empty after the shift. */
211 {
212 int i;
213
214 round_limb = retval[RETURN_LIMB_SIZE - 1];
215 round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
216 for (i = 0; i < RETURN_LIMB_SIZE - 1; ++i)
217 more_bits |= retval[i] != 0;
218 MPN_ZERO (retval, RETURN_LIMB_SIZE);
219 }
220 else if (shift >= BITS_PER_MP_LIMB)
221 {
222 int i;
223
224 round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
225 round_bit = (shift - 1) % BITS_PER_MP_LIMB;
226 for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
227 more_bits |= retval[i] != 0;
228 more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
229 != 0);
230
231 /* __mpn_rshift requires 0 < shift < BITS_PER_MP_LIMB. */
232 if ((shift % BITS_PER_MP_LIMB) != 0)
233 (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
234 RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
235 shift % BITS_PER_MP_LIMB);
236 else
237 for (i = 0; i < RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB); i++)
238 retval[i] = retval[i + (shift / BITS_PER_MP_LIMB)];
239 MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
240 shift / BITS_PER_MP_LIMB);
241 }
242 else if (shift > 0)
243 {
244 if (TININESS_AFTER_ROUNDING && shift == 1)
245 {
246 /* Whether the result counts as tiny depends on whether,
247 after rounding to the normal precision, it still has
248 a subnormal exponent. */
249 mp_limb_t retval_normal[RETURN_LIMB_SIZE];
250 if (round_away (negative,
251 (retval[0] & 1) != 0,
252 (round_limb
253 & (((mp_limb_t) 1) << round_bit)) != 0,
254 (more_bits
255 || ((round_limb
256 & ((((mp_limb_t) 1) << round_bit) - 1))
257 != 0)),
258 mode))
259 {
260 mp_limb_t cy = __mpn_add_1 (retval_normal, retval,
261 RETURN_LIMB_SIZE, 1);
262
263 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
264 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
265 ((retval_normal[RETURN_LIMB_SIZE - 1]
266 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB)))
267 != 0)))
268 is_tiny = false;
269 }
270 }
271 round_limb = retval[0];
272 round_bit = shift - 1;
273 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
274 }
275 /* This is a hook for the m68k long double format, where the
276 exponent bias is the same for normalized and denormalized
277 numbers. */
278#ifndef DENORM_EXP
279# define DENORM_EXP (MIN_EXP - 2)
280#endif
281 exponent = DENORM_EXP;
282 if (is_tiny
283 && ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
284 || more_bits
285 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
286 {
287 __set_errno (ERANGE);
288 FLOAT force_underflow = MIN_VALUE * MIN_VALUE;
289 math_force_eval (force_underflow);
290 }
291 }
292
293 if (exponent > MAX_EXP)
294 goto overflow;
295
296 bool half_bit = (round_limb & (((mp_limb_t) 1) << round_bit)) != 0;
297 bool more_bits_nonzero
298 = (more_bits
299 || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0);
300 if (round_away (negative,
301 (retval[0] & 1) != 0,
302 half_bit,
303 more_bits_nonzero,
304 mode))
305 {
306 mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
307
308 if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
309 ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
310 (retval[RETURN_LIMB_SIZE - 1]
311 & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
312 {
313 ++exponent;
314 (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
315 retval[RETURN_LIMB_SIZE - 1]
316 |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
317 }
318 else if (exponent == DENORM_EXP
319 && (retval[RETURN_LIMB_SIZE - 1]
320 & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
321 != 0)
322 /* The number was denormalized but now normalized. */
323 exponent = MIN_EXP - 1;
324 }
325
326 if (exponent > MAX_EXP)
327 overflow:
328 return overflow_value (negative);
329
330 if (half_bit || more_bits_nonzero)
331 {
332 FLOAT force_inexact = (FLOAT) 1 + MIN_VALUE;
333 math_force_eval (force_inexact);
334 }
335 return MPN2FLOAT (retval, exponent, negative);
336}
337
338
339/* Read a multi-precision integer starting at STR with exactly DIGCNT digits
340 into N. Return the size of the number limbs in NSIZE at the first
341 character od the string that is not part of the integer as the function
342 value. If the EXPONENT is small enough to be taken as an additional
343 factor for the resulting number (see code) multiply by it. */
344static const STRING_TYPE *
345str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
346 intmax_t *exponent
347#ifndef USE_WIDE_CHAR
348 , const char *decimal, size_t decimal_len, const char *thousands
349#endif
350
351 )
352{
353 /* Number of digits for actual limb. */
354 int cnt = 0;
355 mp_limb_t low = 0;
356 mp_limb_t start;
357
358 *nsize = 0;
359 assert (digcnt > 0);
360 do
361 {
362 if (cnt == MAX_DIG_PER_LIMB)
363 {
364 if (*nsize == 0)
365 {
366 n[0] = low;
367 *nsize = 1;
368 }
369 else
370 {
371 mp_limb_t cy;
372 cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
373 cy += __mpn_add_1 (n, n, *nsize, low);
374 if (cy != 0)
375 {
376 assert (*nsize < MPNSIZE);
377 n[*nsize] = cy;
378 ++(*nsize);
379 }
380 }
381 cnt = 0;
382 low = 0;
383 }
384
385 /* There might be thousands separators or radix characters in
386 the string. But these all can be ignored because we know the
387 format of the number is correct and we have an exact number
388 of characters to read. */
389#ifdef USE_WIDE_CHAR
390 if (*str < L'0' || *str > L'9')
391 ++str;
392#else
393 if (*str < '0' || *str > '9')
394 {
395 int inner = 0;
396 if (thousands != NULL && *str == *thousands
397 && ({ for (inner = 1; thousands[inner] != '\0'; ++inner)
398 if (thousands[inner] != str[inner])
399 break;
400 thousands[inner] == '\0'; }))
401 str += inner;
402 else
403 str += decimal_len;
404 }
405#endif
406 low = low * 10 + *str++ - L_('0');
407 ++cnt;
408 }
409 while (--digcnt > 0);
410
411 if (*exponent > 0 && *exponent <= MAX_DIG_PER_LIMB - cnt)
412 {
413 low *= _tens_in_limb[*exponent];
414 start = _tens_in_limb[cnt + *exponent];
415 *exponent = 0;
416 }
417 else
418 start = _tens_in_limb[cnt];
419
420 if (*nsize == 0)
421 {
422 n[0] = low;
423 *nsize = 1;
424 }
425 else
426 {
427 mp_limb_t cy;
428 cy = __mpn_mul_1 (n, n, *nsize, start);
429 cy += __mpn_add_1 (n, n, *nsize, low);
430 if (cy != 0)
431 {
432 assert (*nsize < MPNSIZE);
433 n[(*nsize)++] = cy;
434 }
435 }
436
437 return str;
438}
439
440
441/* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
442 with the COUNT most significant bits of LIMB.
443
444 Implemented as a macro, so that __builtin_constant_p works even at -O0.
445
446 Tege doesn't like this macro so I have to write it here myself. :)
447 --drepper */
448#define __mpn_lshift_1(ptr, size, count, limb) \
449 do \
450 { \
451 mp_limb_t *__ptr = (ptr); \
452 if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \
453 { \
454 mp_size_t i; \
455 for (i = (size) - 1; i > 0; --i) \
456 __ptr[i] = __ptr[i - 1]; \
457 __ptr[0] = (limb); \
458 } \
459 else \
460 { \
461 /* We assume count > 0 && count < BITS_PER_MP_LIMB here. */ \
462 unsigned int __count = (count); \
463 (void) __mpn_lshift (__ptr, __ptr, size, __count); \
464 __ptr[0] |= (limb) >> (BITS_PER_MP_LIMB - __count); \
465 } \
466 } \
467 while (0)
468
469
470#define INTERNAL(x) INTERNAL1(x)
471#define INTERNAL1(x) __##x##_internal
472#ifndef ____STRTOF_INTERNAL
473# define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
474#endif
475
476/* This file defines a function to check for correct grouping. */
477#include "grouping.h"
478
479
480/* Return a floating point number with the value of the given string NPTR.
481 Set *ENDPTR to the character after the last used one. If the number is
482 smaller than the smallest representable number, set `errno' to ERANGE and
483 return 0.0. If the number is too big to be represented, set `errno' to
484 ERANGE and return HUGE_VAL with the appropriate sign. */
485FLOAT
486____STRTOF_INTERNAL (const STRING_TYPE *nptr, STRING_TYPE **endptr, int group,
487 locale_t loc)
488{
489 int negative; /* The sign of the number. */
490 MPN_VAR (num); /* MP representation of the number. */
491 intmax_t exponent; /* Exponent of the number. */
492
493 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
494 int base = 10;
495
496 /* When we have to compute fractional digits we form a fraction with a
497 second multi-precision number (and we sometimes need a second for
498 temporary results). */
499 MPN_VAR (den);
500
501 /* Representation for the return value. */
502 mp_limb_t retval[RETURN_LIMB_SIZE];
503 /* Number of bits currently in result value. */
504 int bits;
505
506 /* Running pointer after the last character processed in the string. */
507 const STRING_TYPE *cp, *tp;
508 /* Start of significant part of the number. */
509 const STRING_TYPE *startp, *start_of_digits;
510 /* Points at the character following the integer and fractional digits. */
511 const STRING_TYPE *expp;
512 /* Total number of digit and number of digits in integer part. */
513 size_t dig_no, int_no, lead_zero;
514 /* Contains the last character read. */
515 CHAR_TYPE c;
516
517/* We should get wint_t from <stddef.h>, but not all GCC versions define it
518 there. So define it ourselves if it remains undefined. */
519#ifndef _WINT_T
520 typedef unsigned int wint_t;
521#endif
522 /* The radix character of the current locale. */
523#ifdef USE_WIDE_CHAR
524 wchar_t decimal;
525#else
526 const char *decimal;
527 size_t decimal_len;
528#endif
529 /* The thousands character of the current locale. */
530#ifdef USE_WIDE_CHAR
531 wchar_t thousands = L'\0';
532#else
533 const char *thousands = NULL;
534#endif
535 /* The numeric grouping specification of the current locale,
536 in the format described in <locale.h>. */
537 const char *grouping;
538 /* Used in several places. */
539 int cnt;
540
541 struct __locale_data *current = loc->__locales[LC_NUMERIC];
542
543 if (__glibc_unlikely (group))
544 {
545 grouping = _NL_CURRENT (LC_NUMERIC, GROUPING);
546 if (*grouping <= 0 || *grouping == CHAR_MAX)
547 grouping = NULL;
548 else
549 {
550 /* Figure out the thousands separator character. */
551#ifdef USE_WIDE_CHAR
552 thousands = _NL_CURRENT_WORD (LC_NUMERIC,
553 _NL_NUMERIC_THOUSANDS_SEP_WC);
554 if (thousands == L'\0')
555 grouping = NULL;
556#else
557 thousands = _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP);
558 if (*thousands == '\0')
559 {
560 thousands = NULL;
561 grouping = NULL;
562 }
563#endif
564 }
565 }
566 else
567 grouping = NULL;
568
569 /* Find the locale's decimal point character. */
570#ifdef USE_WIDE_CHAR
571 decimal = _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_DECIMAL_POINT_WC);
572 assert (decimal != L'\0');
573# define decimal_len 1
574#else
575 decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT);
576 decimal_len = strlen (decimal);
577 assert (decimal_len > 0);
578#endif
579
580 /* Prepare number representation. */
581 exponent = 0;
582 negative = 0;
583 bits = 0;
584
585 /* Parse string to get maximal legal prefix. We need the number of
586 characters of the integer part, the fractional part and the exponent. */
587 cp = nptr - 1;
588 /* Ignore leading white space. */
589 do
590 c = *++cp;
591 while (ISSPACE (c));
592
593 /* Get sign of the result. */
594 if (c == L_('-'))
595 {
596 negative = 1;
597 c = *++cp;
598 }
599 else if (c == L_('+'))
600 c = *++cp;
601
602 /* Return 0.0 if no legal string is found.
603 No character is used even if a sign was found. */
604#ifdef USE_WIDE_CHAR
605 if (c == (wint_t) decimal
606 && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9')
607 {
608 /* We accept it. This funny construct is here only to indent
609 the code correctly. */
610 }
611#else
612 for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
613 if (cp[cnt] != decimal[cnt])
614 break;
615 if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9')
616 {
617 /* We accept it. This funny construct is here only to indent
618 the code correctly. */
619 }
620#endif
621 else if (c < L_('0') || c > L_('9'))
622 {
623 /* Check for `INF' or `INFINITY'. */
624 CHAR_TYPE lowc = TOLOWER_C (c);
625
626 if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0)
627 {
628 /* Return +/- infinity. */
629 if (endptr != NULL)
630 *endptr = (STRING_TYPE *)
631 (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0
632 ? 8 : 3));
633
634 return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
635 }
636
637 if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
638 {
639 /* Return NaN. */
640 FLOAT retval = NAN;
641
642 cp += 3;
643
644 /* Match `(n-char-sequence-digit)'. */
645 if (*cp == L_('('))
646 {
647 const STRING_TYPE *startp = cp;
648 STRING_TYPE *endp;
649 retval = STRTOF_NAN (cp + 1, &endp, L_(')'));
650 if (*endp == L_(')'))
651 /* Consume the closing parenthesis. */
652 cp = endp + 1;
653 else
654 /* Only match the NAN part. */
655 cp = startp;
656 }
657
658 if (endptr != NULL)
659 *endptr = (STRING_TYPE *) cp;
660
661 return retval;
662 }
663
664 /* It is really a text we do not recognize. */
665 RETURN (0.0, nptr);
666 }
667
668 /* First look whether we are faced with a hexadecimal number. */
669 if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
670 {
671 /* Okay, it is a hexa-decimal number. Remember this and skip
672 the characters. BTW: hexadecimal numbers must not be
673 grouped. */
674 base = 16;
675 cp += 2;
676 c = *cp;
677 grouping = NULL;
678 }
679
680 /* Record the start of the digits, in case we will check their grouping. */
681 start_of_digits = startp = cp;
682
683 /* Ignore leading zeroes. This helps us to avoid useless computations. */
684#ifdef USE_WIDE_CHAR
685 while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands))
686 c = *++cp;
687#else
688 if (__glibc_likely (thousands == NULL))
689 while (c == '0')
690 c = *++cp;
691 else
692 {
693 /* We also have the multibyte thousands string. */
694 while (1)
695 {
696 if (c != '0')
697 {
698 for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
699 if (thousands[cnt] != cp[cnt])
700 break;
701 if (thousands[cnt] != '\0')
702 break;
703 cp += cnt - 1;
704 }
705 c = *++cp;
706 }
707 }
708#endif
709
710 /* If no other digit but a '0' is found the result is 0.0.
711 Return current read pointer. */
712 CHAR_TYPE lowc = TOLOWER (c);
713 if (!((c >= L_('0') && c <= L_('9'))
714 || (base == 16 && lowc >= L_('a') && lowc <= L_('f'))
715 || (
716#ifdef USE_WIDE_CHAR
717 c == (wint_t) decimal
718#else
719 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
720 if (decimal[cnt] != cp[cnt])
721 break;
722 decimal[cnt] == '\0'; })
723#endif
724 /* '0x.' alone is not a valid hexadecimal number.
725 '.' alone is not valid either, but that has been checked
726 already earlier. */
727 && (base != 16
728 || cp != start_of_digits
729 || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9'))
730 || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]);
731 lo >= L_('a') && lo <= L_('f'); })))
732 || (base == 16 && (cp != start_of_digits
733 && lowc == L_('p')))
734 || (base != 16 && lowc == L_('e'))))
735 {
736#ifdef USE_WIDE_CHAR
737 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
738 grouping);
739#else
740 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
741 grouping);
742#endif
743 /* If TP is at the start of the digits, there was no correctly
744 grouped prefix of the string; so no number found. */
745 RETURN (negative ? -0.0 : 0.0,
746 tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
747 }
748
749 /* Remember first significant digit and read following characters until the
750 decimal point, exponent character or any non-FP number character. */
751 startp = cp;
752 dig_no = 0;
753 while (1)
754 {
755 if ((c >= L_('0') && c <= L_('9'))
756 || (base == 16
757 && ({ CHAR_TYPE lo = TOLOWER (c);
758 lo >= L_('a') && lo <= L_('f'); })))
759 ++dig_no;
760 else
761 {
762#ifdef USE_WIDE_CHAR
763 if (__builtin_expect ((wint_t) thousands == L'\0', 1)
764 || c != (wint_t) thousands)
765 /* Not a digit or separator: end of the integer part. */
766 break;
767#else
768 if (__glibc_likely (thousands == NULL))
769 break;
770 else
771 {
772 for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
773 if (thousands[cnt] != cp[cnt])
774 break;
775 if (thousands[cnt] != '\0')
776 break;
777 cp += cnt - 1;
778 }
779#endif
780 }
781 c = *++cp;
782 }
783
784 if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits)
785 {
786 /* Check the grouping of the digits. */
787#ifdef USE_WIDE_CHAR
788 tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
789 grouping);
790#else
791 tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
792 grouping);
793#endif
794 if (cp != tp)
795 {
796 /* Less than the entire string was correctly grouped. */
797
798 if (tp == start_of_digits)
799 /* No valid group of numbers at all: no valid number. */
800 RETURN (0.0, nptr);
801
802 if (tp < startp)
803 /* The number is validly grouped, but consists
804 only of zeroes. The whole value is zero. */
805 RETURN (negative ? -0.0 : 0.0, tp);
806
807 /* Recompute DIG_NO so we won't read more digits than
808 are properly grouped. */
809 cp = tp;
810 dig_no = 0;
811 for (tp = startp; tp < cp; ++tp)
812 if (*tp >= L_('0') && *tp <= L_('9'))
813 ++dig_no;
814
815 int_no = dig_no;
816 lead_zero = 0;
817
818 goto number_parsed;
819 }
820 }
821
822 /* We have the number of digits in the integer part. Whether these
823 are all or any is really a fractional digit will be decided
824 later. */
825 int_no = dig_no;
826 lead_zero = int_no == 0 ? (size_t) -1 : 0;
827
828 /* Read the fractional digits. A special case are the 'american
829 style' numbers like `16.' i.e. with decimal point but without
830 trailing digits. */
831 if (
832#ifdef USE_WIDE_CHAR
833 c == (wint_t) decimal
834#else
835 ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
836 if (decimal[cnt] != cp[cnt])
837 break;
838 decimal[cnt] == '\0'; })
839#endif
840 )
841 {
842 cp += decimal_len;
843 c = *cp;
844 while ((c >= L_('0') && c <= L_('9')) ||
845 (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c);
846 lo >= L_('a') && lo <= L_('f'); })))
847 {
848 if (c != L_('0') && lead_zero == (size_t) -1)
849 lead_zero = dig_no - int_no;
850 ++dig_no;
851 c = *++cp;
852 }
853 }
854 assert (dig_no <= (uintmax_t) INTMAX_MAX);
855
856 /* Remember start of exponent (if any). */
857 expp = cp;
858
859 /* Read exponent. */
860 lowc = TOLOWER (c);
861 if ((base == 16 && lowc == L_('p'))
862 || (base != 16 && lowc == L_('e')))
863 {
864 int exp_negative = 0;
865
866 c = *++cp;
867 if (c == L_('-'))
868 {
869 exp_negative = 1;
870 c = *++cp;
871 }
872 else if (c == L_('+'))
873 c = *++cp;
874
875 if (c >= L_('0') && c <= L_('9'))
876 {
877 intmax_t exp_limit;
878
879 /* Get the exponent limit. */
880 if (base == 16)
881 {
882 if (exp_negative)
883 {
884 assert (int_no <= (uintmax_t) (INTMAX_MAX
885 + MIN_EXP - MANT_DIG) / 4);
886 exp_limit = -MIN_EXP + MANT_DIG + 4 * (intmax_t) int_no;
887 }
888 else
889 {
890 if (int_no)
891 {
892 assert (lead_zero == 0
893 && int_no <= (uintmax_t) INTMAX_MAX / 4);
894 exp_limit = MAX_EXP - 4 * (intmax_t) int_no + 3;
895 }
896 else if (lead_zero == (size_t) -1)
897 {
898 /* The number is zero and this limit is
899 arbitrary. */
900 exp_limit = MAX_EXP + 3;
901 }
902 else
903 {
904 assert (lead_zero
905 <= (uintmax_t) (INTMAX_MAX - MAX_EXP - 3) / 4);
906 exp_limit = (MAX_EXP
907 + 4 * (intmax_t) lead_zero
908 + 3);
909 }
910 }
911 }
912 else
913 {
914 if (exp_negative)
915 {
916 assert (int_no
917 <= (uintmax_t) (INTMAX_MAX + MIN_10_EXP - MANT_DIG));
918 exp_limit = -MIN_10_EXP + MANT_DIG + (intmax_t) int_no;
919 }
920 else
921 {
922 if (int_no)
923 {
924 assert (lead_zero == 0
925 && int_no <= (uintmax_t) INTMAX_MAX);
926 exp_limit = MAX_10_EXP - (intmax_t) int_no + 1;
927 }
928 else if (lead_zero == (size_t) -1)
929 {
930 /* The number is zero and this limit is
931 arbitrary. */
932 exp_limit = MAX_10_EXP + 1;
933 }
934 else
935 {
936 assert (lead_zero
937 <= (uintmax_t) (INTMAX_MAX - MAX_10_EXP - 1));
938 exp_limit = MAX_10_EXP + (intmax_t) lead_zero + 1;
939 }
940 }
941 }
942
943 if (exp_limit < 0)
944 exp_limit = 0;
945
946 do
947 {
948 if (__builtin_expect ((exponent > exp_limit / 10
949 || (exponent == exp_limit / 10
950 && c - L_('0') > exp_limit % 10)), 0))
951 /* The exponent is too large/small to represent a valid
952 number. */
953 {
954 FLOAT result;
955
956 /* We have to take care for special situation: a joker
957 might have written "0.0e100000" which is in fact
958 zero. */
959 if (lead_zero == (size_t) -1)
960 result = negative ? -0.0 : 0.0;
961 else
962 {
963 /* Overflow or underflow. */
964 result = (exp_negative
965 ? underflow_value (negative)
966 : overflow_value (negative));
967 }
968
969 /* Accept all following digits as part of the exponent. */
970 do
971 ++cp;
972 while (*cp >= L_('0') && *cp <= L_('9'));
973
974 RETURN (result, cp);
975 /* NOTREACHED */
976 }
977
978 exponent *= 10;
979 exponent += c - L_('0');
980
981 c = *++cp;
982 }
983 while (c >= L_('0') && c <= L_('9'));
984
985 if (exp_negative)
986 exponent = -exponent;
987 }
988 else
989 cp = expp;
990 }
991
992 /* We don't want to have to work with trailing zeroes after the radix. */
993 if (dig_no > int_no)
994 {
995 while (expp[-1] == L_('0'))
996 {
997 --expp;
998 --dig_no;
999 }
1000 assert (dig_no >= int_no);
1001 }
1002
1003 if (dig_no == int_no && dig_no > 0 && exponent < 0)
1004 do
1005 {
1006 while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1])))
1007 --expp;
1008
1009 if (expp[-1] != L_('0'))
1010 break;
1011
1012 --expp;
1013 --dig_no;
1014 --int_no;
1015 exponent += base == 16 ? 4 : 1;
1016 }
1017 while (dig_no > 0 && exponent < 0);
1018
1019 number_parsed:
1020
1021 /* The whole string is parsed. Store the address of the next character. */
1022 if (endptr)
1023 *endptr = (STRING_TYPE *) cp;
1024
1025 if (dig_no == 0)
1026 return negative ? -0.0 : 0.0;
1027
1028 if (lead_zero)
1029 {
1030 /* Find the decimal point */
1031#ifdef USE_WIDE_CHAR
1032 while (*startp != decimal)
1033 ++startp;
1034#else
1035 while (1)
1036 {
1037 if (*startp == decimal[0])
1038 {
1039 for (cnt = 1; decimal[cnt] != '\0'; ++cnt)
1040 if (decimal[cnt] != startp[cnt])
1041 break;
1042 if (decimal[cnt] == '\0')
1043 break;
1044 }
1045 ++startp;
1046 }
1047#endif
1048 startp += lead_zero + decimal_len;
1049 assert (lead_zero <= (base == 16
1050 ? (uintmax_t) INTMAX_MAX / 4
1051 : (uintmax_t) INTMAX_MAX));
1052 assert (lead_zero <= (base == 16
1053 ? ((uintmax_t) exponent
1054 - (uintmax_t) INTMAX_MIN) / 4
1055 : ((uintmax_t) exponent - (uintmax_t) INTMAX_MIN)));
1056 exponent -= base == 16 ? 4 * (intmax_t) lead_zero : (intmax_t) lead_zero;
1057 dig_no -= lead_zero;
1058 }
1059
1060 /* If the BASE is 16 we can use a simpler algorithm. */
1061 if (base == 16)
1062 {
1063 static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
1064 4, 4, 4, 4, 4, 4, 4, 4 };
1065 int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
1066 int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
1067 mp_limb_t val;
1068
1069 while (!ISXDIGIT (*startp))
1070 ++startp;
1071 while (*startp == L_('0'))
1072 ++startp;
1073 if (ISDIGIT (*startp))
1074 val = *startp++ - L_('0');
1075 else
1076 val = 10 + TOLOWER (*startp++) - L_('a');
1077 bits = nbits[val];
1078 /* We cannot have a leading zero. */
1079 assert (bits != 0);
1080
1081 if (pos + 1 >= 4 || pos + 1 >= bits)
1082 {
1083 /* We don't have to care for wrapping. This is the normal
1084 case so we add the first clause in the `if' expression as
1085 an optimization. It is a compile-time constant and so does
1086 not cost anything. */
1087 retval[idx] = val << (pos - bits + 1);
1088 pos -= bits;
1089 }
1090 else
1091 {
1092 retval[idx--] = val >> (bits - pos - 1);
1093 retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
1094 pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
1095 }
1096
1097 /* Adjust the exponent for the bits we are shifting in. */
1098 assert (int_no <= (uintmax_t) (exponent < 0
1099 ? (INTMAX_MAX - bits + 1) / 4
1100 : (INTMAX_MAX - exponent - bits + 1) / 4));
1101 exponent += bits - 1 + ((intmax_t) int_no - 1) * 4;
1102
1103 while (--dig_no > 0 && idx >= 0)
1104 {
1105 if (!ISXDIGIT (*startp))
1106 startp += decimal_len;
1107 if (ISDIGIT (*startp))
1108 val = *startp++ - L_('0');
1109 else
1110 val = 10 + TOLOWER (*startp++) - L_('a');
1111
1112 if (pos + 1 >= 4)
1113 {
1114 retval[idx] |= val << (pos - 4 + 1);
1115 pos -= 4;
1116 }
1117 else
1118 {
1119 retval[idx--] |= val >> (4 - pos - 1);
1120 val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
1121 if (idx < 0)
1122 {
1123 int rest_nonzero = 0;
1124 while (--dig_no > 0)
1125 {
1126 if (*startp != L_('0'))
1127 {
1128 rest_nonzero = 1;
1129 break;
1130 }
1131 startp++;
1132 }
1133 return round_and_return (retval, exponent, negative, val,
1134 BITS_PER_MP_LIMB - 1, rest_nonzero);
1135 }
1136
1137 retval[idx] = val;
1138 pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
1139 }
1140 }
1141
1142 /* We ran out of digits. */
1143 MPN_ZERO (retval, idx);
1144
1145 return round_and_return (retval, exponent, negative, 0, 0, 0);
1146 }
1147
1148 /* Now we have the number of digits in total and the integer digits as well
1149 as the exponent and its sign. We can decide whether the read digits are
1150 really integer digits or belong to the fractional part; i.e. we normalize
1151 123e-2 to 1.23. */
1152 {
1153 intmax_t incr = (exponent < 0
1154 ? MAX (-(intmax_t) int_no, exponent)
1155 : MIN ((intmax_t) dig_no - (intmax_t) int_no, exponent));
1156 int_no += incr;
1157 exponent -= incr;
1158 }
1159
1160 if (__glibc_unlikely (exponent > MAX_10_EXP + 1 - (intmax_t) int_no))
1161 return overflow_value (negative);
1162
1163 /* 10^(MIN_10_EXP-1) is not normal. Thus, 10^(MIN_10_EXP-1) /
1164 2^MANT_DIG is below half the least subnormal, so anything with a
1165 base-10 exponent less than the base-10 exponent (which is
1166 MIN_10_EXP - 1 - ceil(MANT_DIG*log10(2))) of that value
1167 underflows. DIG is floor((MANT_DIG-1)log10(2)), so an exponent
1168 below MIN_10_EXP - (DIG + 3) underflows. But EXPONENT is
1169 actually an exponent multiplied only by a fractional part, not an
1170 integer part, so an exponent below MIN_10_EXP - (DIG + 2)
1171 underflows. */
1172 if (__glibc_unlikely (exponent < MIN_10_EXP - (DIG + 2)))
1173 return underflow_value (negative);
1174
1175 if (int_no > 0)
1176 {
1177 /* Read the integer part as a multi-precision number to NUM. */
1178 startp = str_to_mpn (startp, int_no, num, &numsize, &exponent
1179#ifndef USE_WIDE_CHAR
1180 , decimal, decimal_len, thousands
1181#endif
1182 );
1183
1184 if (exponent > 0)
1185 {
1186 /* We now multiply the gained number by the given power of ten. */
1187 mp_limb_t *psrc = num;
1188 mp_limb_t *pdest = den;
1189 int expbit = 1;
1190 const struct mp_power *ttab = &_fpioconst_pow10[0];
1191
1192 do
1193 {
1194 if ((exponent & expbit) != 0)
1195 {
1196 size_t size = ttab->arraysize - _FPIO_CONST_OFFSET;
1197 mp_limb_t cy;
1198 exponent ^= expbit;
1199
1200 /* FIXME: not the whole multiplication has to be
1201 done. If we have the needed number of bits we
1202 only need the information whether more non-zero
1203 bits follow. */
1204 if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
1205 cy = __mpn_mul (pdest, psrc, numsize,
1206 &__tens[ttab->arrayoff
1207 + _FPIO_CONST_OFFSET],
1208 size);
1209 else
1210 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
1211 + _FPIO_CONST_OFFSET],
1212 size, psrc, numsize);
1213 numsize += size;
1214 if (cy == 0)
1215 --numsize;
1216 (void) SWAP (psrc, pdest);
1217 }
1218 expbit <<= 1;
1219 ++ttab;
1220 }
1221 while (exponent != 0);
1222
1223 if (psrc == den)
1224 memcpy (num, den, numsize * sizeof (mp_limb_t));
1225 }
1226
1227 /* Determine how many bits of the result we already have. */
1228 count_leading_zeros (bits, num[numsize - 1]);
1229 bits = numsize * BITS_PER_MP_LIMB - bits;
1230
1231 /* Now we know the exponent of the number in base two.
1232 Check it against the maximum possible exponent. */
1233 if (__glibc_unlikely (bits > MAX_EXP))
1234 return overflow_value (negative);
1235
1236 /* We have already the first BITS bits of the result. Together with
1237 the information whether more non-zero bits follow this is enough
1238 to determine the result. */
1239 if (bits > MANT_DIG)
1240 {
1241 int i;
1242 const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
1243 const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
1244 const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
1245 : least_idx;
1246 const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
1247 : least_bit - 1;
1248
1249 if (least_bit == 0)
1250 memcpy (retval, &num[least_idx],
1251 RETURN_LIMB_SIZE * sizeof (mp_limb_t));
1252 else
1253 {
1254 for (i = least_idx; i < numsize - 1; ++i)
1255 retval[i - least_idx] = (num[i] >> least_bit)
1256 | (num[i + 1]
1257 << (BITS_PER_MP_LIMB - least_bit));
1258 if (i - least_idx < RETURN_LIMB_SIZE)
1259 retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
1260 }
1261
1262 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1263 for (i = 0; num[i] == 0; ++i)
1264 ;
1265
1266 return round_and_return (retval, bits - 1, negative,
1267 num[round_idx], round_bit,
1268 int_no < dig_no || i < round_idx);
1269 /* NOTREACHED */
1270 }
1271 else if (dig_no == int_no)
1272 {
1273 const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
1274 const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
1275
1276 if (target_bit == is_bit)
1277 {
1278 memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
1279 numsize * sizeof (mp_limb_t));
1280 /* FIXME: the following loop can be avoided if we assume a
1281 maximal MANT_DIG value. */
1282 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1283 }
1284 else if (target_bit > is_bit)
1285 {
1286 (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
1287 num, numsize, target_bit - is_bit);
1288 /* FIXME: the following loop can be avoided if we assume a
1289 maximal MANT_DIG value. */
1290 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
1291 }
1292 else
1293 {
1294 mp_limb_t cy;
1295 assert (numsize < RETURN_LIMB_SIZE);
1296
1297 cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
1298 num, numsize, is_bit - target_bit);
1299 retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
1300 /* FIXME: the following loop can be avoided if we assume a
1301 maximal MANT_DIG value. */
1302 MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
1303 }
1304
1305 return round_and_return (retval, bits - 1, negative, 0, 0, 0);
1306 /* NOTREACHED */
1307 }
1308
1309 /* Store the bits we already have. */
1310 memcpy (retval, num, numsize * sizeof (mp_limb_t));
1311#if RETURN_LIMB_SIZE > 1
1312 if (numsize < RETURN_LIMB_SIZE)
1313# if RETURN_LIMB_SIZE == 2
1314 retval[numsize] = 0;
1315# else
1316 MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize);
1317# endif
1318#endif
1319 }
1320
1321 /* We have to compute at least some of the fractional digits. */
1322 {
1323 /* We construct a fraction and the result of the division gives us
1324 the needed digits. The denominator is 1.0 multiplied by the
1325 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1326 123e-6 gives 123 / 1000000. */
1327
1328 int expbit;
1329 int neg_exp;
1330 int more_bits;
1331 int need_frac_digits;
1332 mp_limb_t cy;
1333 mp_limb_t *psrc = den;
1334 mp_limb_t *pdest = num;
1335 const struct mp_power *ttab = &_fpioconst_pow10[0];
1336
1337 assert (dig_no > int_no
1338 && exponent <= 0
1339 && exponent >= MIN_10_EXP - (DIG + 2));
1340
1341 /* We need to compute MANT_DIG - BITS fractional bits that lie
1342 within the mantissa of the result, the following bit for
1343 rounding, and to know whether any subsequent bit is 0.
1344 Computing a bit with value 2^-n means looking at n digits after
1345 the decimal point. */
1346 if (bits > 0)
1347 {
1348 /* The bits required are those immediately after the point. */
1349 assert (int_no > 0 && exponent == 0);
1350 need_frac_digits = 1 + MANT_DIG - bits;
1351 }
1352 else
1353 {
1354 /* The number is in the form .123eEXPONENT. */
1355 assert (int_no == 0 && *startp != L_('0'));
1356 /* The number is at least 10^(EXPONENT-1), and 10^3 <
1357 2^10. */
1358 int neg_exp_2 = ((1 - exponent) * 10) / 3 + 1;
1359 /* The number is at least 2^-NEG_EXP_2. We need up to
1360 MANT_DIG bits following that bit. */
1361 need_frac_digits = neg_exp_2 + MANT_DIG;
1362 /* However, we never need bits beyond 1/4 ulp of the smallest
1363 representable value. (That 1/4 ulp bit is only needed to
1364 determine tinyness on machines where tinyness is determined
1365 after rounding.) */
1366 if (need_frac_digits > MANT_DIG - MIN_EXP + 2)
1367 need_frac_digits = MANT_DIG - MIN_EXP + 2;
1368 /* At this point, NEED_FRAC_DIGITS is the total number of
1369 digits needed after the point, but some of those may be
1370 leading 0s. */
1371 need_frac_digits += exponent;
1372 /* Any cases underflowing enough that none of the fractional
1373 digits are needed should have been caught earlier (such
1374 cases are on the order of 10^-n or smaller where 2^-n is
1375 the least subnormal). */
1376 assert (need_frac_digits > 0);
1377 }
1378
1379 if (need_frac_digits > (intmax_t) dig_no - (intmax_t) int_no)
1380 need_frac_digits = (intmax_t) dig_no - (intmax_t) int_no;
1381
1382 if ((intmax_t) dig_no > (intmax_t) int_no + need_frac_digits)
1383 {
1384 dig_no = int_no + need_frac_digits;
1385 more_bits = 1;
1386 }
1387 else
1388 more_bits = 0;
1389
1390 neg_exp = (intmax_t) dig_no - (intmax_t) int_no - exponent;
1391
1392 /* Construct the denominator. */
1393 densize = 0;
1394 expbit = 1;
1395 do
1396 {
1397 if ((neg_exp & expbit) != 0)
1398 {
1399 mp_limb_t cy;
1400 neg_exp ^= expbit;
1401
1402 if (densize == 0)
1403 {
1404 densize = ttab->arraysize - _FPIO_CONST_OFFSET;
1405 memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET],
1406 densize * sizeof (mp_limb_t));
1407 }
1408 else
1409 {
1410 cy = __mpn_mul (pdest, &__tens[ttab->arrayoff
1411 + _FPIO_CONST_OFFSET],
1412 ttab->arraysize - _FPIO_CONST_OFFSET,
1413 psrc, densize);
1414 densize += ttab->arraysize - _FPIO_CONST_OFFSET;
1415 if (cy == 0)
1416 --densize;
1417 (void) SWAP (psrc, pdest);
1418 }
1419 }
1420 expbit <<= 1;
1421 ++ttab;
1422 }
1423 while (neg_exp != 0);
1424
1425 if (psrc == num)
1426 memcpy (den, num, densize * sizeof (mp_limb_t));
1427
1428 /* Read the fractional digits from the string. */
1429 (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent
1430#ifndef USE_WIDE_CHAR
1431 , decimal, decimal_len, thousands
1432#endif
1433 );
1434
1435 /* We now have to shift both numbers so that the highest bit in the
1436 denominator is set. In the same process we copy the numerator to
1437 a high place in the array so that the division constructs the wanted
1438 digits. This is done by a "quasi fix point" number representation.
1439
1440 num: ddddddddddd . 0000000000000000000000
1441 |--- m ---|
1442 den: ddddddddddd n >= m
1443 |--- n ---|
1444 */
1445
1446 count_leading_zeros (cnt, den[densize - 1]);
1447
1448 if (cnt > 0)
1449 {
1450 /* Don't call `mpn_shift' with a count of zero since the specification
1451 does not allow this. */
1452 (void) __mpn_lshift (den, den, densize, cnt);
1453 cy = __mpn_lshift (num, num, numsize, cnt);
1454 if (cy != 0)
1455 num[numsize++] = cy;
1456 }
1457
1458 /* Now we are ready for the division. But it is not necessary to
1459 do a full multi-precision division because we only need a small
1460 number of bits for the result. So we do not use __mpn_divmod
1461 here but instead do the division here by hand and stop whenever
1462 the needed number of bits is reached. The code itself comes
1463 from the GNU MP Library by Torbj\"orn Granlund. */
1464
1465 exponent = bits;
1466
1467 switch (densize)
1468 {
1469 case 1:
1470 {
1471 mp_limb_t d, n, quot;
1472 int used = 0;
1473
1474 n = num[0];
1475 d = den[0];
1476 assert (numsize == 1 && n < d);
1477
1478 do
1479 {
1480 udiv_qrnnd (quot, n, n, 0, d);
1481
1482#define got_limb \
1483 if (bits == 0) \
1484 { \
1485 int cnt; \
1486 if (quot == 0) \
1487 cnt = BITS_PER_MP_LIMB; \
1488 else \
1489 count_leading_zeros (cnt, quot); \
1490 exponent -= cnt; \
1491 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1492 { \
1493 used = MANT_DIG + cnt; \
1494 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1495 bits = MANT_DIG + 1; \
1496 } \
1497 else \
1498 { \
1499 /* Note that we only clear the second element. */ \
1500 /* The conditional is determined at compile time. */ \
1501 if (RETURN_LIMB_SIZE > 1) \
1502 retval[1] = 0; \
1503 retval[0] = quot; \
1504 bits = -cnt; \
1505 } \
1506 } \
1507 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1508 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1509 quot); \
1510 else \
1511 { \
1512 used = MANT_DIG - bits; \
1513 if (used > 0) \
1514 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1515 } \
1516 bits += BITS_PER_MP_LIMB
1517
1518 got_limb;
1519 }
1520 while (bits <= MANT_DIG);
1521
1522 return round_and_return (retval, exponent - 1, negative,
1523 quot, BITS_PER_MP_LIMB - 1 - used,
1524 more_bits || n != 0);
1525 }
1526 case 2:
1527 {
1528 mp_limb_t d0, d1, n0, n1;
1529 mp_limb_t quot = 0;
1530 int used = 0;
1531
1532 d0 = den[0];
1533 d1 = den[1];
1534
1535 if (numsize < densize)
1536 {
1537 if (num[0] >= d1)
1538 {
1539 /* The numerator of the number occupies fewer bits than
1540 the denominator but the one limb is bigger than the
1541 high limb of the numerator. */
1542 n1 = 0;
1543 n0 = num[0];
1544 }
1545 else
1546 {
1547 if (bits <= 0)
1548 exponent -= BITS_PER_MP_LIMB;
1549 else
1550 {
1551 if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
1552 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1553 BITS_PER_MP_LIMB, 0);
1554 else
1555 {
1556 used = MANT_DIG - bits;
1557 if (used > 0)
1558 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1559 }
1560 bits += BITS_PER_MP_LIMB;
1561 }
1562 n1 = num[0];
1563 n0 = 0;
1564 }
1565 }
1566 else
1567 {
1568 n1 = num[1];
1569 n0 = num[0];
1570 }
1571
1572 while (bits <= MANT_DIG)
1573 {
1574 mp_limb_t r;
1575
1576 if (n1 == d1)
1577 {
1578 /* QUOT should be either 111..111 or 111..110. We need
1579 special treatment of this rare case as normal division
1580 would give overflow. */
1581 quot = ~(mp_limb_t) 0;
1582
1583 r = n0 + d1;
1584 if (r < d1) /* Carry in the addition? */
1585 {
1586 add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
1587 goto have_quot;
1588 }
1589 n1 = d0 - (d0 != 0);
1590 n0 = -d0;
1591 }
1592 else
1593 {
1594 udiv_qrnnd (quot, r, n1, n0, d1);
1595 umul_ppmm (n1, n0, d0, quot);
1596 }
1597
1598 q_test:
1599 if (n1 > r || (n1 == r && n0 > 0))
1600 {
1601 /* The estimated QUOT was too large. */
1602 --quot;
1603
1604 sub_ddmmss (n1, n0, n1, n0, 0, d0);
1605 r += d1;
1606 if (r >= d1) /* If not carry, test QUOT again. */
1607 goto q_test;
1608 }
1609 sub_ddmmss (n1, n0, r, 0, n1, n0);
1610
1611 have_quot:
1612 got_limb;
1613 }
1614
1615 return round_and_return (retval, exponent - 1, negative,
1616 quot, BITS_PER_MP_LIMB - 1 - used,
1617 more_bits || n1 != 0 || n0 != 0);
1618 }
1619 default:
1620 {
1621 int i;
1622 mp_limb_t cy, dX, d1, n0, n1;
1623 mp_limb_t quot = 0;
1624 int used = 0;
1625
1626 dX = den[densize - 1];
1627 d1 = den[densize - 2];
1628
1629 /* The division does not work if the upper limb of the two-limb
1630 numerator is greater than the denominator. */
1631 if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0)
1632 num[numsize++] = 0;
1633
1634 if (numsize < densize)
1635 {
1636 mp_size_t empty = densize - numsize;
1637 int i;
1638
1639 if (bits <= 0)
1640 exponent -= empty * BITS_PER_MP_LIMB;
1641 else
1642 {
1643 if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
1644 {
1645 /* We make a difference here because the compiler
1646 cannot optimize the `else' case that good and
1647 this reflects all currently used FLOAT types
1648 and GMP implementations. */
1649#if RETURN_LIMB_SIZE <= 2
1650 assert (empty == 1);
1651 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
1652 BITS_PER_MP_LIMB, 0);
1653#else
1654 for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i)
1655 retval[i] = retval[i - empty];
1656 while (i >= 0)
1657 retval[i--] = 0;
1658#endif
1659 }
1660 else
1661 {
1662 used = MANT_DIG - bits;
1663 if (used >= BITS_PER_MP_LIMB)
1664 {
1665 int i;
1666 (void) __mpn_lshift (&retval[used
1667 / BITS_PER_MP_LIMB],
1668 retval,
1669 (RETURN_LIMB_SIZE
1670 - used / BITS_PER_MP_LIMB),
1671 used % BITS_PER_MP_LIMB);
1672 for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i)
1673 retval[i] = 0;
1674 }
1675 else if (used > 0)
1676 __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
1677 }
1678 bits += empty * BITS_PER_MP_LIMB;
1679 }
1680 for (i = numsize; i > 0; --i)
1681 num[i + empty] = num[i - 1];
1682 MPN_ZERO (num, empty + 1);
1683 }
1684 else
1685 {
1686 int i;
1687 assert (numsize == densize);
1688 for (i = numsize; i > 0; --i)
1689 num[i] = num[i - 1];
1690 num[0] = 0;
1691 }
1692
1693 den[densize] = 0;
1694 n0 = num[densize];
1695
1696 while (bits <= MANT_DIG)
1697 {
1698 if (n0 == dX)
1699 /* This might over-estimate QUOT, but it's probably not
1700 worth the extra code here to find out. */
1701 quot = ~(mp_limb_t) 0;
1702 else
1703 {
1704 mp_limb_t r;
1705
1706 udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
1707 umul_ppmm (n1, n0, d1, quot);
1708
1709 while (n1 > r || (n1 == r && n0 > num[densize - 2]))
1710 {
1711 --quot;
1712 r += dX;
1713 if (r < dX) /* I.e. "carry in previous addition?" */
1714 break;
1715 n1 -= n0 < d1;
1716 n0 -= d1;
1717 }
1718 }
1719
1720 /* Possible optimization: We already have (q * n0) and (1 * n1)
1721 after the calculation of QUOT. Taking advantage of this, we
1722 could make this loop make two iterations less. */
1723
1724 cy = __mpn_submul_1 (num, den, densize + 1, quot);
1725
1726 if (num[densize] != cy)
1727 {
1728 cy = __mpn_add_n (num, num, den, densize);
1729 assert (cy != 0);
1730 --quot;
1731 }
1732 n0 = num[densize] = num[densize - 1];
1733 for (i = densize - 1; i > 0; --i)
1734 num[i] = num[i - 1];
1735 num[0] = 0;
1736
1737 got_limb;
1738 }
1739
1740 for (i = densize; i >= 0 && num[i] == 0; --i)
1741 ;
1742 return round_and_return (retval, exponent - 1, negative,
1743 quot, BITS_PER_MP_LIMB - 1 - used,
1744 more_bits || i >= 0);
1745 }
1746 }
1747 }
1748
1749 /* NOTREACHED */
1750}
1751#if defined _LIBC && !defined USE_WIDE_CHAR
1752libc_hidden_def (____STRTOF_INTERNAL)
1753#endif
1754
1755/* External user entry point. */
1756
1757FLOAT
1758#ifdef weak_function
1759weak_function
1760#endif
1761__STRTOF (const STRING_TYPE *nptr, STRING_TYPE **endptr, locale_t loc)
1762{
1763 return ____STRTOF_INTERNAL (nptr, endptr, 0, loc);
1764}
1765#if defined _LIBC
1766libc_hidden_def (__STRTOF)
1767libc_hidden_ver (__STRTOF, STRTOF)
1768#endif
1769weak_alias (__STRTOF, STRTOF)
1770
1771#ifdef LONG_DOUBLE_COMPAT
1772# if LONG_DOUBLE_COMPAT(libc, GLIBC_2_1)
1773# ifdef USE_WIDE_CHAR
1774compat_symbol (libc, __wcstod_l, __wcstold_l, GLIBC_2_1);
1775# else
1776compat_symbol (libc, __strtod_l, __strtold_l, GLIBC_2_1);
1777# endif
1778# endif
1779# if LONG_DOUBLE_COMPAT(libc, GLIBC_2_3)
1780# ifdef USE_WIDE_CHAR
1781compat_symbol (libc, wcstod_l, wcstold_l, GLIBC_2_3);
1782# else
1783compat_symbol (libc, strtod_l, strtold_l, GLIBC_2_3);
1784# endif
1785# endif
1786#endif
1787