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
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
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
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	last_digit_odd: (retval[0] & 1) != 0,
271	half_bit: (round_limb
272	& (((mp_limb_t) 1) << round_bit)) != 0,
273	more_bits: (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	last_digit_odd: (retval[0] & 1) != 0,
322	half_bit,
323	more_bits: 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 (mpn: 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 *
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
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 (begin: start_of_digits, end: 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 (begin: start_of_digits, end: 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, round_limb: val,
1154	BITS_PER_MP_LIMB - 1, more_bits: 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, round_limb: 0, round_bit: 0, more_bits: 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 (str: startp, digcnt: int_no, n: num, nsize: &numsize, exponent: &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, exponent: bits - 1, negative,
1287	round_limb: num[round_idx], round_bit,
1288	more_bits: 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, exponent: bits - 1, negative, round_limb: 0, round_bit: 0, more_bits: 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 (str: startp, digcnt: dig_no - int_no, n: num, nsize: &numsize, exponent: &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: exponent - 1, negative,
1543	round_limb: quot, BITS_PER_MP_LIMB - 1 - used,
1544	more_bits: 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: exponent - 1, negative,
1636	round_limb: quot, BITS_PER_MP_LIMB - 1 - used,
1637	more_bits: 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: exponent - 1, negative,
1763	round_limb: quot, BITS_PER_MP_LIMB - 1 - used,
1764	more_bits: more_bits || i >= 0);
1765	}
1766	}
1767	}
1768
1769	/* NOTREACHED */
1770	}
1771	#if defined _LIBC && !defined USE_WIDE_CHAR
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, group: 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