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