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