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