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.
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.
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.
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/>. */
38 #include "quadmath-rounding-mode.h"
40 #include "../printf/quadmath-printf.h"
41 #include "../printf/fpioconst.h"
45 # define STRING_TYPE wchar_t
46 # define CHAR_TYPE wint_t
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)
57 # define STRING_TYPE char
58 # define CHAR_TYPE char
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)
70 # define STRTOULL(S, E, B) strtoull (S, E, B)
72 # define STRTOULL(S, E, B) strtoul (S, E, B)
76 __quadmath_strncasecmp_c (const char *s1
, const char *s2
, size_t n
)
78 const unsigned char *p1
= (const unsigned char *) s1
;
79 const unsigned char *p2
= (const unsigned char *) s2
;
81 if (p1
== p2
|| n
== 0)
83 while ((result
= TOLOWER_C (*p1
) - TOLOWER_C (*p2
++)) == 0)
84 if (*p1
++ == '\0' || --n
== 0)
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)
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
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
);
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
118 # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
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
;
125 #define howmany(x,y) (((x)+((y)-1))/(y))
127 #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
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)
133 #define RETURN(val,end) \
134 do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
135 return val; } while (0)
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))
149 /* Set errno and return an overflowing value with sign specified by
152 overflow_value (int negative
)
154 #if defined HAVE_ERRNO_H && defined ERANGE
157 FLOAT result
= (negative
? -MAX_VALUE
: MAX_VALUE
) * MAX_VALUE
;
161 /* Set errno and return an underflowing value with sign specified by
164 underflow_value (int negative
)
166 #if defined HAVE_ERRNO_H && defined ERANGE
169 FLOAT result
= (negative
? -MIN_VALUE
: MIN_VALUE
) * MIN_VALUE
;
173 /* Return a floating point number of the needed type according to the given
174 multi-precision number after possible rounding. */
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
)
180 int mode
= get_rounding_mode ();
183 if (exponent
< MIN_EXP
- 1)
188 if (exponent
< MIN_EXP
- 1 - MANT_DIG
)
189 return underflow_value (negative
);
191 shift
= MIN_EXP
- 1 - exponent
;
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. */
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
);
207 else if (shift
>= BITS_PER_MP_LIMB
)
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))
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
);
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
);
232 if (TININESS_AFTER_ROUNDING
&& shift
== 1)
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,
241 & (((mp_limb_t
) 1) << round_bit
)) != 0,
244 & ((((mp_limb_t
) 1) << round_bit
) - 1))
248 mp_limb_t cy
= mpn_add_1 (retval_normal
, retval
,
249 RETURN_LIMB_SIZE
, 1);
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
)))
260 round_limb
= retval
[0];
261 round_bit
= shift
- 1;
262 (void) mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
264 /* This is a hook for the m68k long double format, where the
265 exponent bias is the same for normalized and denormalized
268 # define DENORM_EXP (MIN_EXP - 2)
270 exponent
= DENORM_EXP
;
272 && ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
274 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
276 #if defined HAVE_ERRNO_H && defined ERANGE
279 volatile FLOAT force_underflow_exception
= MIN_VALUE
* MIN_VALUE
;
280 (void) force_underflow_exception
;
284 if (exponent
>= MAX_EXP
)
288 if (round_away (negative
,
289 (retval
[0] & 1) != 0,
290 (round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0,
292 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0),
295 mp_limb_t cy
= mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
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))
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
);
307 else if (exponent
== DENORM_EXP
308 && (retval
[RETURN_LIMB_SIZE
- 1]
309 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
311 /* The number was denormalized but now normalized. */
312 exponent
= MIN_EXP
- 1;
316 if (exponent
>= MAX_EXP
)
318 return overflow_value (negative
);
320 return MPN2FLOAT (retval
, exponent
, negative
);
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
*
330 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
332 #ifndef USE_WIDE_CHAR
333 , const char *decimal
, size_t decimal_len
, const char *thousands
338 /* Number of digits for actual limb. */
347 if (cnt
== MAX_DIG_PER_LIMB
)
357 cy
= mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
358 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
361 assert (*nsize
< MPNSIZE
);
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. */
375 if (*str
< L
'0' || *str
> L
'9')
378 if (*str
< '0' || *str
> '9')
381 if (thousands
!= NULL
&& *str
== *thousands
382 && ({ for (inner
= 1; thousands
[inner
] != '\0'; ++inner
)
383 if (thousands
[inner
] != str
[inner
])
385 thousands
[inner
] == '\0'; }))
391 low
= low
* 10 + *str
++ - L_('0');
394 while (--digcnt
> 0);
396 if (*exponent
> 0 && *exponent
<= MAX_DIG_PER_LIMB
- cnt
)
398 low
*= _tens_in_limb
[*exponent
];
399 start
= _tens_in_limb
[cnt
+ *exponent
];
403 start
= _tens_in_limb
[cnt
];
413 cy
= mpn_mul_1 (n
, n
, *nsize
, start
);
414 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
417 assert (*nsize
< MPNSIZE
);
426 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
427 with the COUNT most significant bits of LIMB.
429 Implemented as a macro, so that __builtin_constant_p works even at -O0.
431 Tege doesn't like this macro so I have to write it here myself. :)
433 #define mpn_lshift_1(ptr, size, count, limb) \
436 mp_limb_t *__ptr = (ptr); \
437 if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \
440 for (i = (size) - 1; i > 0; --i) \
441 __ptr[i] = __ptr[i - 1]; \
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); \
455 #define INTERNAL(x) INTERNAL1(x)
456 #define INTERNAL1(x) __##x##_internal
457 #ifndef ____STRTOF_INTERNAL
458 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
461 /* This file defines a function to check for correct grouping. */
462 #include "grouping.h"
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. */
471 ____STRTOF_INTERNAL (nptr
, endptr
, group
)
472 const STRING_TYPE
*nptr
;
473 STRING_TYPE
**endptr
;
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. */
480 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
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). */
488 /* Representation for the return value. */
489 mp_limb_t retval
[RETURN_LIMB_SIZE
];
490 /* Number of bits currently in result value. */
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. */
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. */
507 typedef unsigned int wint_t;
509 /* The radix character of the current locale. */
516 /* The thousands character of the current locale. */
518 wchar_t thousands
= L
'\0';
520 const char *thousands
= NULL
;
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. */
528 #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
529 const struct lconv
*lc
= localeconv ();
532 if (__builtin_expect (group
, 0))
534 #ifdef USE_NL_LANGINFO
535 grouping
= nl_langinfo (GROUPING
);
536 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
540 /* Figure out the thousands separator character. */
542 thousands
= nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC
);
543 if (thousands
== L
'\0')
546 thousands
= nl_langinfo (THOUSANDS_SEP
);
547 if (*thousands
== '\0')
554 #elif defined USE_LOCALECONV
555 grouping
= lc
->grouping
;
556 if (grouping
== NULL
|| *grouping
<= 0 || *grouping
== CHAR_MAX
)
560 /* Figure out the thousands separator character. */
561 thousands
= lc
->thousands_sep
;
562 if (thousands
== NULL
|| *thousands
== '\0')
575 /* Find the locale's decimal point character. */
577 decimal
= nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC
);
578 assert (decimal
!= L
'\0');
579 # define decimal_len 1
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')
589 decimal_len
= strlen (decimal
);
596 /* Prepare number representation. */
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. */
604 /* Ignore leading white space. */
609 /* Get sign of the result. */
615 else if (c
== L_('+'))
618 /* Return 0.0 if no legal string is found.
619 No character is used even if a sign was found. */
621 if (c
== (wint_t) decimal
622 && (wint_t) cp
[1] >= L
'0' && (wint_t) cp
[1] <= L
'9')
624 /* We accept it. This funny construct is here only to indent
625 the code correctly. */
628 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
629 if (cp
[cnt
] != decimal
[cnt
])
631 if (decimal
[cnt
] == '\0' && cp
[cnt
] >= '0' && cp
[cnt
] <= '9')
633 /* We accept it. This funny construct is here only to indent
634 the code correctly. */
637 else if (c
< L_('0') || c
> L_('9'))
639 /* Check for `INF' or `INFINITY'. */
640 CHAR_TYPE lowc
= TOLOWER_C (c
);
642 if (lowc
== L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
644 /* Return +/- infinity. */
646 *endptr
= (STRING_TYPE
*)
647 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
650 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
653 if (lowc
== L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
660 /* Match `(n-char-sequence-digit)'. */
663 const STRING_TYPE
*startp
= cp
;
666 while ((*cp
>= L_('0') && *cp
<= L_('9'))
667 || ({ CHAR_TYPE lo
= TOLOWER (*cp
);
668 lo
>= L_('a') && lo
<= L_('z'); })
672 /* The closing brace is missing. Only match the NAN
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
682 unsigned long long int mant
;
684 mant
= STRTOULL (startp
+ 1, &endp
, 0);
686 SET_MANTISSA (retval
, mant
);
688 /* Consume the closing brace. */
694 *endptr
= (STRING_TYPE
*) cp
;
696 return negative
? -retval
: retval
;
699 /* It is really a text we do not recognize. */
703 /* First look whether we are faced with a hexadecimal number. */
704 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
706 /* Okay, it is a hexa-decimal number. Remember this and skip
707 the characters. BTW: hexadecimal numbers must not be
715 /* Record the start of the digits, in case we will check their grouping. */
716 start_of_digits
= startp
= cp
;
718 /* Ignore leading zeroes. This helps us to avoid useless computations. */
720 while (c
== L
'0' || ((wint_t) thousands
!= L
'\0' && c
== (wint_t) thousands
))
723 if (__builtin_expect (thousands
== NULL
, 1))
728 /* We also have the multibyte thousands string. */
733 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
734 if (thousands
[cnt
] != cp
[cnt
])
736 if (thousands
[cnt
] != '\0')
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'))
752 c
== (wint_t) decimal
754 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
755 if (decimal
[cnt
] != cp
[cnt
])
757 decimal
[cnt
] == '\0'; })
759 /* '0x.' alone is not a valid hexadecimal number.
760 '.' alone is not valid either, but that has been checked
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
769 || (base
!= 16 && lowc
== L_('e'))))
772 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
775 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
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
);
784 /* Remember first significant digit and read following characters until the
785 decimal point, exponent character or any non-FP number character. */
790 if ((c
>= L_('0') && c
<= L_('9'))
792 && ({ CHAR_TYPE lo
= TOLOWER (c
);
793 lo
>= L_('a') && lo
<= L_('f'); })))
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. */
803 if (__builtin_expect (thousands
== NULL
, 1))
807 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
808 if (thousands
[cnt
] != cp
[cnt
])
810 if (thousands
[cnt
] != '\0')
819 if (__builtin_expect (grouping
!= NULL
, 0) && cp
> start_of_digits
)
821 /* Check the grouping of the digits. */
823 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
826 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
831 /* Less than the entire string was correctly grouped. */
833 if (tp
== start_of_digits
)
834 /* No valid group of numbers at all: no valid number. */
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
);
842 /* Recompute DIG_NO so we won't read more digits than
843 are properly grouped. */
846 for (tp
= startp
; tp
< cp
; ++tp
)
847 if (*tp
>= L_('0') && *tp
<= L_('9'))
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
861 lead_zero
= int_no
== 0 ? (size_t) -1 : 0;
863 /* Read the fractional digits. A special case are the 'american
864 style' numbers like `16.' i.e. with decimal point but without
868 c
== (wint_t) decimal
870 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
871 if (decimal
[cnt
] != cp
[cnt
])
873 decimal
[cnt
] == '\0'; })
879 while ((c
>= L_('0') && c
<= L_('9')) ||
880 (base
== 16 && ({ CHAR_TYPE lo
= TOLOWER (c
);
881 lo
>= L_('a') && lo
<= L_('f'); })))
883 if (c
!= L_('0') && lead_zero
== (size_t) -1)
884 lead_zero
= dig_no
- int_no
;
889 assert (dig_no
<= (uintmax_t) INTMAX_MAX
);
891 /* Remember start of exponent (if any). */
896 if ((base
== 16 && lowc
== L_('p'))
897 || (base
!= 16 && lowc
== L_('e')))
899 int exp_negative
= 0;
907 else if (c
== L_('+'))
910 if (c
>= L_('0') && c
<= L_('9'))
914 /* Get the exponent limit. */
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
;
927 assert (lead_zero
== 0
928 && int_no
<= (uintmax_t) INTMAX_MAX
/ 4);
929 exp_limit
= MAX_EXP
- 4 * (intmax_t) int_no
+ 3;
931 else if (lead_zero
== (size_t) -1)
933 /* The number is zero and this limit is
935 exp_limit
= MAX_EXP
+ 3;
940 <= (uintmax_t) (INTMAX_MAX
- MAX_EXP
- 3) / 4);
942 + 4 * (intmax_t) lead_zero
952 <= (uintmax_t) (INTMAX_MAX
+ MIN_10_EXP
- MANT_DIG
));
953 exp_limit
= -MIN_10_EXP
+ MANT_DIG
+ (intmax_t) int_no
;
959 assert (lead_zero
== 0
960 && int_no
<= (uintmax_t) INTMAX_MAX
);
961 exp_limit
= MAX_10_EXP
- (intmax_t) int_no
+ 1;
963 else if (lead_zero
== (size_t) -1)
965 /* The number is zero and this limit is
967 exp_limit
= MAX_10_EXP
+ 1;
972 <= (uintmax_t) (INTMAX_MAX
- MAX_10_EXP
- 1));
973 exp_limit
= MAX_10_EXP
+ (intmax_t) lead_zero
+ 1;
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
991 /* We have to take care for special situation: a joker
992 might have written "0.0e100000" which is in fact
994 if (lead_zero
== (size_t) -1)
995 result
= negative
? -0.0 : 0.0;
998 /* Overflow or underflow. */
999 #if defined HAVE_ERRNO_H && defined ERANGE
1002 result
= (exp_negative
? (negative
? -0.0 : 0.0) :
1003 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
1006 /* Accept all following digits as part of the exponent. */
1009 while (*cp
>= L_('0') && *cp
<= L_('9'));
1011 RETURN (result
, cp
);
1016 exponent
+= c
- L_('0');
1020 while (c
>= L_('0') && c
<= L_('9'));
1023 exponent
= -exponent
;
1029 /* We don't want to have to work with trailing zeroes after the radix. */
1030 if (dig_no
> int_no
)
1032 while (expp
[-1] == L_('0'))
1037 assert (dig_no
>= int_no
);
1040 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
1043 while (! (base
== 16 ? ISXDIGIT (expp
[-1]) : ISDIGIT (expp
[-1])))
1046 if (expp
[-1] != L_('0'))
1052 exponent
+= base
== 16 ? 4 : 1;
1054 while (dig_no
> 0 && exponent
< 0);
1058 /* The whole string is parsed. Store the address of the next character. */
1060 *endptr
= (STRING_TYPE
*) cp
;
1063 return negative
? -0.0 : 0.0;
1067 /* Find the decimal point */
1068 #ifdef USE_WIDE_CHAR
1069 while (*startp
!= decimal
)
1074 if (*startp
== decimal
[0])
1076 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
1077 if (decimal
[cnt
] != startp
[cnt
])
1079 if (decimal
[cnt
] == '\0')
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
;
1097 /* If the BASE is 16 we can use a simpler algorithm. */
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
;
1106 while (!ISXDIGIT (*startp
))
1108 while (*startp
== L_('0'))
1110 if (ISDIGIT (*startp
))
1111 val
= *startp
++ - L_('0');
1113 val
= 10 + TOLOWER (*startp
++) - L_('a');
1115 /* We cannot have a leading zero. */
1118 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
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);
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);
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;
1140 while (--dig_no
> 0 && idx
>= 0)
1142 if (!ISXDIGIT (*startp
))
1143 startp
+= decimal_len
;
1144 if (ISDIGIT (*startp
))
1145 val
= *startp
++ - L_('0');
1147 val
= 10 + TOLOWER (*startp
++) - L_('a');
1151 retval
[idx
] |= val
<< (pos
- 4 + 1);
1156 retval
[idx
--] |= val
>> (4 - pos
- 1);
1157 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
1160 int rest_nonzero
= 0;
1161 while (--dig_no
> 0)
1163 if (*startp
!= L_('0'))
1170 return round_and_return (retval
, exponent
, negative
, val
,
1171 BITS_PER_MP_LIMB
- 1, rest_nonzero
);
1175 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
1179 /* We ran out of digits. */
1180 MPN_ZERO (retval
, idx
);
1182 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
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
1190 register intmax_t incr
= (exponent
< 0
1191 ? MAX (-(intmax_t) int_no
, exponent
)
1192 : MIN ((intmax_t) dig_no
- (intmax_t) int_no
,
1198 if (__builtin_expect (exponent
> MAX_10_EXP
+ 1 - (intmax_t) int_no
, 0))
1199 return overflow_value (negative
);
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)
1210 if (__builtin_expect (exponent
< MIN_10_EXP
- (DIG
+ 2), 0))
1211 return underflow_value (negative
);
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
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
;
1228 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1232 if ((exponent
& expbit
) != 0)
1234 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
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
1242 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1243 cy
= mpn_mul (pdest
, psrc
, numsize
,
1244 &__tens
[ttab
->arrayoff
1245 + _FPIO_CONST_OFFSET
],
1248 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1249 + _FPIO_CONST_OFFSET
],
1250 size
, psrc
, numsize
);
1254 (void) SWAP (psrc
, pdest
);
1259 while (exponent
!= 0);
1262 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
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
;
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
);
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
)
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
1284 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1288 memcpy (retval
, &num
[least_idx
],
1289 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1292 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1293 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
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
;
1300 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1301 for (i
= 0; num
[i
] == 0; ++i
)
1304 return round_and_return (retval
, bits
- 1, negative
,
1305 num
[round_idx
], round_bit
,
1306 int_no
< dig_no
|| i
< round_idx
);
1309 else if (dig_no
== int_no
)
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
;
1314 if (target_bit
== is_bit
)
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
);
1322 else if (target_bit
> is_bit
)
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
);
1333 assert (numsize
< RETURN_LIMB_SIZE
);
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);
1343 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
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;
1354 MPN_ZERO (retval
+ numsize
, RETURN_LIMB_SIZE
- numsize
);
1359 /* We have to compute at least some of the fractional digits. */
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. */
1369 int need_frac_digits
;
1371 mp_limb_t
*psrc
= den
;
1372 mp_limb_t
*pdest
= num
;
1373 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1375 assert (dig_no
> int_no
1377 && exponent
>= MIN_10_EXP
- (DIG
+ 2));
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. */
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
;
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 <
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
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
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);
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
;
1420 if ((intmax_t) dig_no
> (intmax_t) int_no
+ need_frac_digits
)
1422 dig_no
= int_no
+ need_frac_digits
;
1428 neg_exp
= (intmax_t) dig_no
- (intmax_t) int_no
- exponent
;
1430 /* Construct the denominator. */
1435 if ((neg_exp
& expbit
) != 0)
1442 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1443 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1444 densize
* sizeof (mp_limb_t
));
1448 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1449 + _FPIO_CONST_OFFSET
],
1450 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1452 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1455 (void) SWAP (psrc
, pdest
);
1461 while (neg_exp
!= 0);
1464 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
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
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.
1478 num: ddddddddddd . 0000000000000000000000
1480 den: ddddddddddd n >= m
1484 count_leading_zeros (cnt
, den
[densize
- 1]);
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
);
1493 num
[numsize
++] = cy
;
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. */
1509 mp_limb_t d
, n
, quot
;
1514 assert (numsize
== 1 && n
< d
);
1518 udiv_qrnnd (quot
, n
, n
, 0, d
);
1525 cnt = BITS_PER_MP_LIMB; \
1527 count_leading_zeros (cnt, quot); \
1529 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1531 used = MANT_DIG + cnt; \
1532 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1533 bits = MANT_DIG + 1; \
1537 /* Note that we only clear the second element. */ \
1538 /* The conditional is determined at compile time. */ \
1539 if (RETURN_LIMB_SIZE > 1) \
1545 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1546 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1550 used = MANT_DIG - bits; \
1552 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1554 bits += BITS_PER_MP_LIMB
1558 while (bits
<= MANT_DIG
);
1560 return round_and_return (retval
, exponent
- 1, negative
,
1561 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1562 more_bits
|| n
!= 0);
1566 mp_limb_t d0
, d1
, n0
, n1
;
1573 if (numsize
< densize
)
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. */
1586 exponent
-= BITS_PER_MP_LIMB
;
1589 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1590 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1591 BITS_PER_MP_LIMB
, 0);
1594 used
= MANT_DIG
- bits
;
1596 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1598 bits
+= BITS_PER_MP_LIMB
;
1610 while (bits
<= MANT_DIG
)
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;
1622 if (r
< d1
) /* Carry in the addition? */
1624 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1627 n1
= d0
- (d0
!= 0);
1632 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1633 umul_ppmm (n1
, n0
, d0
, quot
);
1637 if (n1
> r
|| (n1
== r
&& n0
> 0))
1639 /* The estimated QUOT was too large. */
1642 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1644 if (r
>= d1
) /* If not carry, test QUOT again. */
1647 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1653 return round_and_return (retval
, exponent
- 1, negative
,
1654 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1655 more_bits
|| n1
!= 0 || n0
!= 0);
1660 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1664 dX
= den
[densize
- 1];
1665 d1
= den
[densize
- 2];
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)
1672 if (numsize
< densize
)
1674 mp_size_t empty
= densize
- numsize
;
1678 exponent
-= empty
* BITS_PER_MP_LIMB
;
1681 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
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);
1692 for (i
= RETURN_LIMB_SIZE
- 1; i
>= empty
; --i
)
1693 retval
[i
] = retval
[i
- empty
];
1700 used
= MANT_DIG
- bits
;
1701 if (used
>= BITS_PER_MP_LIMB
)
1704 (void) mpn_lshift (&retval
[used
1705 / BITS_PER_MP_LIMB
],
1708 - used
/ BITS_PER_MP_LIMB
),
1709 used
% BITS_PER_MP_LIMB
);
1710 for (i
= used
/ BITS_PER_MP_LIMB
- 1; i
>= 0; --i
)
1714 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1716 bits
+= empty
* BITS_PER_MP_LIMB
;
1718 for (i
= numsize
; i
> 0; --i
)
1719 num
[i
+ empty
] = num
[i
- 1];
1720 MPN_ZERO (num
, empty
+ 1);
1725 assert (numsize
== densize
);
1726 for (i
= numsize
; i
> 0; --i
)
1727 num
[i
] = num
[i
- 1];
1734 while (bits
<= MANT_DIG
)
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;
1744 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1745 umul_ppmm (n1
, n0
, d1
, quot
);
1747 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1751 if (r
< dX
) /* I.e. "carry in previous addition?" */
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. */
1762 cy
= mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1764 if (num
[densize
] != cy
)
1766 cy
= mpn_add_n (num
, num
, den
, densize
);
1770 n0
= num
[densize
] = num
[densize
- 1];
1771 for (i
= densize
- 1; i
> 0; --i
)
1772 num
[i
] = num
[i
- 1];
1778 for (i
= densize
; i
>= 0 && num
[i
] == 0; --i
)
1780 return round_and_return (retval
, exponent
- 1, negative
,
1781 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1782 more_bits
|| i
>= 0);