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