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