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