1 /* Convert string representing a number to float value, using given locale.
2 Copyright (C) 1997-2013 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
; ++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 (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
);
227 if (TININESS_AFTER_ROUNDING
&& shift
== 1)
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,
236 & (((mp_limb_t
) 1) << round_bit
)) != 0,
239 & ((((mp_limb_t
) 1) << round_bit
) - 1))
243 mp_limb_t cy
= mpn_add_1 (retval_normal
, retval
,
244 RETURN_LIMB_SIZE
, 1);
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
)))
255 round_limb
= retval
[0];
256 round_bit
= shift
- 1;
257 (void) mpn_rshift (retval
, retval
, RETURN_LIMB_SIZE
, shift
);
259 /* This is a hook for the m68k long double format, where the
260 exponent bias is the same for normalized and denormalized
263 # define DENORM_EXP (MIN_EXP - 2)
265 exponent
= DENORM_EXP
;
267 && ((round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0
269 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0))
271 #if defined HAVE_ERRNO_H && defined ERANGE
274 volatile FLOAT force_underflow_exception
= MIN_VALUE
* MIN_VALUE
;
275 (void) force_underflow_exception
;
279 if (exponent
> MAX_EXP
)
283 if (round_away (negative
,
284 (retval
[0] & 1) != 0,
285 (round_limb
& (((mp_limb_t
) 1) << round_bit
)) != 0,
287 || (round_limb
& ((((mp_limb_t
) 1) << round_bit
) - 1)) != 0),
290 mp_limb_t cy
= mpn_add_1 (retval
, retval
, RETURN_LIMB_SIZE
, 1);
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))
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
);
302 else if (exponent
== DENORM_EXP
303 && (retval
[RETURN_LIMB_SIZE
- 1]
304 & (((mp_limb_t
) 1) << ((MANT_DIG
- 1) % BITS_PER_MP_LIMB
)))
306 /* The number was denormalized but now normalized. */
307 exponent
= MIN_EXP
- 1;
311 if (exponent
> MAX_EXP
)
313 return overflow_value (negative
);
315 return MPN2FLOAT (retval
, exponent
, negative
);
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. */
324 static const STRING_TYPE
*
325 str_to_mpn (const STRING_TYPE
*str
, int digcnt
, mp_limb_t
*n
, mp_size_t
*nsize
,
327 #ifndef USE_WIDE_CHAR
328 , const char *decimal
, size_t decimal_len
, const char *thousands
333 /* Number of digits for actual limb. */
342 if (cnt
== MAX_DIG_PER_LIMB
)
352 cy
= mpn_mul_1 (n
, n
, *nsize
, MAX_FAC_PER_LIMB
);
353 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
356 assert (*nsize
< MPNSIZE
);
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. */
370 if (*str
< L
'0' || *str
> L
'9')
373 if (*str
< '0' || *str
> '9')
376 if (thousands
!= NULL
&& *str
== *thousands
377 && ({ for (inner
= 1; thousands
[inner
] != '\0'; ++inner
)
378 if (thousands
[inner
] != str
[inner
])
380 thousands
[inner
] == '\0'; }))
386 low
= low
* 10 + *str
++ - L_('0');
389 while (--digcnt
> 0);
391 if (*exponent
> 0 && *exponent
<= MAX_DIG_PER_LIMB
- cnt
)
393 low
*= _tens_in_limb
[*exponent
];
394 start
= _tens_in_limb
[cnt
+ *exponent
];
398 start
= _tens_in_limb
[cnt
];
408 cy
= mpn_mul_1 (n
, n
, *nsize
, start
);
409 cy
+= mpn_add_1 (n
, n
, *nsize
, low
);
412 assert (*nsize
< MPNSIZE
);
421 /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
422 with the COUNT most significant bits of LIMB.
424 Tege doesn't like this function so I have to write it here myself. :)
427 __attribute ((always_inline
))
428 mpn_lshift_1 (mp_limb_t
*ptr
, mp_size_t size
, unsigned int count
,
431 if (__builtin_constant_p (count
) && count
== BITS_PER_MP_LIMB
)
433 /* Optimize the case of shifting by exactly a word:
434 just copy words, with no actual bit-shifting. */
436 for (i
= size
- 1; i
> 0; --i
)
442 (void) mpn_lshift (ptr
, ptr
, size
, count
);
443 ptr
[0] |= limb
>> (BITS_PER_MP_LIMB
- count
);
448 #define INTERNAL(x) INTERNAL1(x)
449 #define INTERNAL1(x) __##x##_internal
450 #ifndef ____STRTOF_INTERNAL
451 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
454 /* This file defines a function to check for correct grouping. */
455 #include "grouping.h"
458 /* Return a floating point number with the value of the given string NPTR.
459 Set *ENDPTR to the character after the last used one. If the number is
460 smaller than the smallest representable number, set `errno' to ERANGE and
461 return 0.0. If the number is too big to be represented, set `errno' to
462 ERANGE and return HUGE_VAL with the appropriate sign. */
464 ____STRTOF_INTERNAL (nptr
, endptr
, group
)
465 const STRING_TYPE
*nptr
;
466 STRING_TYPE
**endptr
;
469 int negative
; /* The sign of the number. */
470 MPN_VAR (num
); /* MP representation of the number. */
471 intmax_t exponent
; /* Exponent of the number. */
473 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
476 /* When we have to compute fractional digits we form a fraction with a
477 second multi-precision number (and we sometimes need a second for
478 temporary results). */
481 /* Representation for the return value. */
482 mp_limb_t retval
[RETURN_LIMB_SIZE
];
483 /* Number of bits currently in result value. */
486 /* Running pointer after the last character processed in the string. */
487 const STRING_TYPE
*cp
, *tp
;
488 /* Start of significant part of the number. */
489 const STRING_TYPE
*startp
, *start_of_digits
;
490 /* Points at the character following the integer and fractional digits. */
491 const STRING_TYPE
*expp
;
492 /* Total number of digit and number of digits in integer part. */
493 size_t dig_no
, int_no
, lead_zero
;
494 /* Contains the last character read. */
497 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
498 there. So define it ourselves if it remains undefined. */
500 typedef unsigned int wint_t;
502 /* The radix character of the current locale. */
509 /* The thousands character of the current locale. */
511 wchar_t thousands
= L
'\0';
513 const char *thousands
= NULL
;
515 /* The numeric grouping specification of the current locale,
516 in the format described in <locale.h>. */
517 const char *grouping
;
518 /* Used in several places. */
521 #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
522 const struct lconv
*lc
= localeconv ();
525 if (__builtin_expect (group
, 0))
527 #ifdef USE_NL_LANGINFO
528 grouping
= nl_langinfo (GROUPING
);
529 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
533 /* Figure out the thousands separator character. */
535 thousands
= nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC
);
536 if (thousands
== L
'\0')
539 thousands
= nl_langinfo (THOUSANDS_SEP
);
540 if (*thousands
== '\0')
547 #elif defined USE_LOCALECONV
548 grouping
= lc
->grouping
;
549 if (grouping
== NULL
|| *grouping
<= 0 || *grouping
== CHAR_MAX
)
553 /* Figure out the thousands separator character. */
554 thousands
= lc
->thousands_sep
;
555 if (thousands
== NULL
|| *thousands
== '\0')
568 /* Find the locale's decimal point character. */
570 decimal
= nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC
);
571 assert (decimal
!= L
'\0');
572 # define decimal_len 1
574 #ifdef USE_NL_LANGINFO
575 decimal
= nl_langinfo (DECIMAL_POINT
);
576 decimal_len
= strlen (decimal
);
577 assert (decimal_len
> 0);
578 #elif defined USE_LOCALECONV
579 decimal
= lc
->decimal_point
;
580 if (decimal
== NULL
|| *decimal
== '\0')
582 decimal_len
= strlen (decimal
);
589 /* Prepare number representation. */
594 /* Parse string to get maximal legal prefix. We need the number of
595 characters of the integer part, the fractional part and the exponent. */
597 /* Ignore leading white space. */
602 /* Get sign of the result. */
608 else if (c
== L_('+'))
611 /* Return 0.0 if no legal string is found.
612 No character is used even if a sign was found. */
614 if (c
== (wint_t) decimal
615 && (wint_t) cp
[1] >= L
'0' && (wint_t) cp
[1] <= L
'9')
617 /* We accept it. This funny construct is here only to indent
618 the code correctly. */
621 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
622 if (cp
[cnt
] != decimal
[cnt
])
624 if (decimal
[cnt
] == '\0' && cp
[cnt
] >= '0' && cp
[cnt
] <= '9')
626 /* We accept it. This funny construct is here only to indent
627 the code correctly. */
630 else if (c
< L_('0') || c
> L_('9'))
632 /* Check for `INF' or `INFINITY'. */
633 CHAR_TYPE lowc
= TOLOWER_C (c
);
635 if (lowc
== L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
637 /* Return +/- infinity. */
639 *endptr
= (STRING_TYPE
*)
640 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
643 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
646 if (lowc
== L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
653 /* Match `(n-char-sequence-digit)'. */
656 const STRING_TYPE
*startp
= cp
;
659 while ((*cp
>= L_('0') && *cp
<= L_('9'))
660 || ({ CHAR_TYPE lo
= TOLOWER (*cp
);
661 lo
>= L_('a') && lo
<= L_('z'); })
665 /* The closing brace is missing. Only match the NAN
670 /* This is a system-dependent way to specify the
671 bitmask used for the NaN. We expect it to be
672 a number which is put in the mantissa of the
675 unsigned long long int mant
;
677 mant
= STRTOULL (startp
+ 1, &endp
, 0);
679 SET_MANTISSA (retval
, mant
);
681 /* Consume the closing brace. */
687 *endptr
= (STRING_TYPE
*) cp
;
692 /* It is really a text we do not recognize. */
696 /* First look whether we are faced with a hexadecimal number. */
697 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
699 /* Okay, it is a hexa-decimal number. Remember this and skip
700 the characters. BTW: hexadecimal numbers must not be
708 /* Record the start of the digits, in case we will check their grouping. */
709 start_of_digits
= startp
= cp
;
711 /* Ignore leading zeroes. This helps us to avoid useless computations. */
713 while (c
== L
'0' || ((wint_t) thousands
!= L
'\0' && c
== (wint_t) thousands
))
716 if (__builtin_expect (thousands
== NULL
, 1))
721 /* We also have the multibyte thousands string. */
726 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
727 if (thousands
[cnt
] != cp
[cnt
])
729 if (thousands
[cnt
] != '\0')
738 /* If no other digit but a '0' is found the result is 0.0.
739 Return current read pointer. */
740 CHAR_TYPE lowc
= TOLOWER (c
);
741 if (!((c
>= L_('0') && c
<= L_('9'))
742 || (base
== 16 && lowc
>= L_('a') && lowc
<= L_('f'))
745 c
== (wint_t) decimal
747 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
748 if (decimal
[cnt
] != cp
[cnt
])
750 decimal
[cnt
] == '\0'; })
752 /* '0x.' alone is not a valid hexadecimal number.
753 '.' alone is not valid either, but that has been checked
756 || cp
!= start_of_digits
757 || (cp
[decimal_len
] >= L_('0') && cp
[decimal_len
] <= L_('9'))
758 || ({ CHAR_TYPE lo
= TOLOWER (cp
[decimal_len
]);
759 lo
>= L_('a') && lo
<= L_('f'); })))
760 || (base
== 16 && (cp
!= start_of_digits
762 || (base
!= 16 && lowc
== L_('e'))))
765 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
768 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
771 /* If TP is at the start of the digits, there was no correctly
772 grouped prefix of the string; so no number found. */
773 RETURN (negative
? -0.0 : 0.0,
774 tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
777 /* Remember first significant digit and read following characters until the
778 decimal point, exponent character or any non-FP number character. */
783 if ((c
>= L_('0') && c
<= L_('9'))
785 && ({ CHAR_TYPE lo
= TOLOWER (c
);
786 lo
>= L_('a') && lo
<= L_('f'); })))
791 if (__builtin_expect ((wint_t) thousands
== L
'\0', 1)
792 || c
!= (wint_t) thousands
)
793 /* Not a digit or separator: end of the integer part. */
796 if (__builtin_expect (thousands
== NULL
, 1))
800 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
801 if (thousands
[cnt
] != cp
[cnt
])
803 if (thousands
[cnt
] != '\0')
812 if (__builtin_expect (grouping
!= NULL
, 0) && cp
> start_of_digits
)
814 /* Check the grouping of the digits. */
816 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
819 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
824 /* Less than the entire string was correctly grouped. */
826 if (tp
== start_of_digits
)
827 /* No valid group of numbers at all: no valid number. */
831 /* The number is validly grouped, but consists
832 only of zeroes. The whole value is zero. */
833 RETURN (negative
? -0.0 : 0.0, tp
);
835 /* Recompute DIG_NO so we won't read more digits than
836 are properly grouped. */
839 for (tp
= startp
; tp
< cp
; ++tp
)
840 if (*tp
>= L_('0') && *tp
<= L_('9'))
850 /* We have the number of digits in the integer part. Whether these
851 are all or any is really a fractional digit will be decided
854 lead_zero
= int_no
== 0 ? (size_t) -1 : 0;
856 /* Read the fractional digits. A special case are the 'american
857 style' numbers like `16.' i.e. with decimal point but without
861 c
== (wint_t) decimal
863 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
864 if (decimal
[cnt
] != cp
[cnt
])
866 decimal
[cnt
] == '\0'; })
872 while ((c
>= L_('0') && c
<= L_('9')) ||
873 (base
== 16 && ({ CHAR_TYPE lo
= TOLOWER (c
);
874 lo
>= L_('a') && lo
<= L_('f'); })))
876 if (c
!= L_('0') && lead_zero
== (size_t) -1)
877 lead_zero
= dig_no
- int_no
;
882 assert (dig_no
<= (uintmax_t) INTMAX_MAX
);
884 /* Remember start of exponent (if any). */
889 if ((base
== 16 && lowc
== L_('p'))
890 || (base
!= 16 && lowc
== L_('e')))
892 int exp_negative
= 0;
900 else if (c
== L_('+'))
903 if (c
>= L_('0') && c
<= L_('9'))
907 /* Get the exponent limit. */
912 assert (int_no
<= (uintmax_t) (INTMAX_MAX
913 + MIN_EXP
- MANT_DIG
) / 4);
914 exp_limit
= -MIN_EXP
+ MANT_DIG
+ 4 * (intmax_t) int_no
;
920 assert (lead_zero
== 0
921 && int_no
<= (uintmax_t) INTMAX_MAX
/ 4);
922 exp_limit
= MAX_EXP
- 4 * (intmax_t) int_no
+ 3;
924 else if (lead_zero
== (size_t) -1)
926 /* The number is zero and this limit is
928 exp_limit
= MAX_EXP
+ 3;
933 <= (uintmax_t) (INTMAX_MAX
- MAX_EXP
- 3) / 4);
935 + 4 * (intmax_t) lead_zero
945 <= (uintmax_t) (INTMAX_MAX
+ MIN_10_EXP
- MANT_DIG
));
946 exp_limit
= -MIN_10_EXP
+ MANT_DIG
+ (intmax_t) int_no
;
952 assert (lead_zero
== 0
953 && int_no
<= (uintmax_t) INTMAX_MAX
);
954 exp_limit
= MAX_10_EXP
- (intmax_t) int_no
+ 1;
956 else if (lead_zero
== (size_t) -1)
958 /* The number is zero and this limit is
960 exp_limit
= MAX_10_EXP
+ 1;
965 <= (uintmax_t) (INTMAX_MAX
- MAX_10_EXP
- 1));
966 exp_limit
= MAX_10_EXP
+ (intmax_t) lead_zero
+ 1;
976 if (__builtin_expect ((exponent
> exp_limit
/ 10
977 || (exponent
== exp_limit
/ 10
978 && c
- L_('0') > exp_limit
% 10)), 0))
979 /* The exponent is too large/small to represent a valid
984 /* We have to take care for special situation: a joker
985 might have written "0.0e100000" which is in fact
987 if (lead_zero
== (size_t) -1)
988 result
= negative
? -0.0 : 0.0;
991 /* Overflow or underflow. */
992 #if defined HAVE_ERRNO_H && defined ERANGE
995 result
= (exp_negative
? (negative
? -0.0 : 0.0) :
996 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
999 /* Accept all following digits as part of the exponent. */
1002 while (*cp
>= L_('0') && *cp
<= L_('9'));
1004 RETURN (result
, cp
);
1009 exponent
+= c
- L_('0');
1013 while (c
>= L_('0') && c
<= L_('9'));
1016 exponent
= -exponent
;
1022 /* We don't want to have to work with trailing zeroes after the radix. */
1023 if (dig_no
> int_no
)
1025 while (expp
[-1] == L_('0'))
1030 assert (dig_no
>= int_no
);
1033 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
1036 while (! (base
== 16 ? ISXDIGIT (expp
[-1]) : ISDIGIT (expp
[-1])))
1039 if (expp
[-1] != L_('0'))
1045 exponent
+= base
== 16 ? 4 : 1;
1047 while (dig_no
> 0 && exponent
< 0);
1051 /* The whole string is parsed. Store the address of the next character. */
1053 *endptr
= (STRING_TYPE
*) cp
;
1056 return negative
? -0.0 : 0.0;
1060 /* Find the decimal point */
1061 #ifdef USE_WIDE_CHAR
1062 while (*startp
!= decimal
)
1067 if (*startp
== decimal
[0])
1069 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
1070 if (decimal
[cnt
] != startp
[cnt
])
1072 if (decimal
[cnt
] == '\0')
1078 startp
+= lead_zero
+ decimal_len
;
1079 assert (lead_zero
<= (base
== 16
1080 ? (uintmax_t) INTMAX_MAX
/ 4
1081 : (uintmax_t) INTMAX_MAX
));
1082 assert (lead_zero
<= (base
== 16
1083 ? ((uintmax_t) exponent
1084 - (uintmax_t) INTMAX_MIN
) / 4
1085 : ((uintmax_t) exponent
- (uintmax_t) INTMAX_MIN
)));
1086 exponent
-= base
== 16 ? 4 * (intmax_t) lead_zero
: (intmax_t) lead_zero
;
1087 dig_no
-= lead_zero
;
1090 /* If the BASE is 16 we can use a simpler algorithm. */
1093 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
1094 4, 4, 4, 4, 4, 4, 4, 4 };
1095 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
1096 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1099 while (!ISXDIGIT (*startp
))
1101 while (*startp
== L_('0'))
1103 if (ISDIGIT (*startp
))
1104 val
= *startp
++ - L_('0');
1106 val
= 10 + TOLOWER (*startp
++) - L_('a');
1108 /* We cannot have a leading zero. */
1111 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
1113 /* We don't have to care for wrapping. This is the normal
1114 case so we add the first clause in the `if' expression as
1115 an optimization. It is a compile-time constant and so does
1116 not cost anything. */
1117 retval
[idx
] = val
<< (pos
- bits
+ 1);
1122 retval
[idx
--] = val
>> (bits
- pos
- 1);
1123 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
1124 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
1127 /* Adjust the exponent for the bits we are shifting in. */
1128 assert (int_no
<= (uintmax_t) (exponent
< 0
1129 ? (INTMAX_MAX
- bits
+ 1) / 4
1130 : (INTMAX_MAX
- exponent
- bits
+ 1) / 4));
1131 exponent
+= bits
- 1 + ((intmax_t) int_no
- 1) * 4;
1133 while (--dig_no
> 0 && idx
>= 0)
1135 if (!ISXDIGIT (*startp
))
1136 startp
+= decimal_len
;
1137 if (ISDIGIT (*startp
))
1138 val
= *startp
++ - L_('0');
1140 val
= 10 + TOLOWER (*startp
++) - L_('a');
1144 retval
[idx
] |= val
<< (pos
- 4 + 1);
1149 retval
[idx
--] |= val
>> (4 - pos
- 1);
1150 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
1153 int rest_nonzero
= 0;
1154 while (--dig_no
> 0)
1156 if (*startp
!= L_('0'))
1163 return round_and_return (retval
, exponent
, negative
, val
,
1164 BITS_PER_MP_LIMB
- 1, rest_nonzero
);
1168 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
1172 /* We ran out of digits. */
1173 MPN_ZERO (retval
, idx
);
1175 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
1178 /* Now we have the number of digits in total and the integer digits as well
1179 as the exponent and its sign. We can decide whether the read digits are
1180 really integer digits or belong to the fractional part; i.e. we normalize
1183 register intmax_t incr
= (exponent
< 0
1184 ? MAX (-(intmax_t) int_no
, exponent
)
1185 : MIN ((intmax_t) dig_no
- (intmax_t) int_no
,
1191 if (__builtin_expect (exponent
> MAX_10_EXP
+ 1 - (intmax_t) int_no
, 0))
1192 return overflow_value (negative
);
1194 if (__builtin_expect (exponent
< MIN_10_EXP
- (DIG
+ 1), 0))
1195 return underflow_value (negative
);
1199 /* Read the integer part as a multi-precision number to NUM. */
1200 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1201 #ifndef USE_WIDE_CHAR
1202 , decimal
, decimal_len
, thousands
1208 /* We now multiply the gained number by the given power of ten. */
1209 mp_limb_t
*psrc
= num
;
1210 mp_limb_t
*pdest
= den
;
1212 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1216 if ((exponent
& expbit
) != 0)
1218 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1222 /* FIXME: not the whole multiplication has to be
1223 done. If we have the needed number of bits we
1224 only need the information whether more non-zero
1226 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1227 cy
= mpn_mul (pdest
, psrc
, numsize
,
1228 &__tens
[ttab
->arrayoff
1229 + _FPIO_CONST_OFFSET
],
1232 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1233 + _FPIO_CONST_OFFSET
],
1234 size
, psrc
, numsize
);
1238 (void) SWAP (psrc
, pdest
);
1243 while (exponent
!= 0);
1246 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1249 /* Determine how many bits of the result we already have. */
1250 count_leading_zeros (bits
, num
[numsize
- 1]);
1251 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1253 /* Now we know the exponent of the number in base two.
1254 Check it against the maximum possible exponent. */
1255 if (__builtin_expect (bits
> MAX_EXP
, 0))
1256 return overflow_value (negative
);
1258 /* We have already the first BITS bits of the result. Together with
1259 the information whether more non-zero bits follow this is enough
1260 to determine the result. */
1261 if (bits
> MANT_DIG
)
1264 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1265 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1266 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1268 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1272 memcpy (retval
, &num
[least_idx
],
1273 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1276 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1277 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1279 << (BITS_PER_MP_LIMB
- least_bit
));
1280 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1281 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1284 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1285 for (i
= 0; num
[i
] == 0; ++i
)
1288 return round_and_return (retval
, bits
- 1, negative
,
1289 num
[round_idx
], round_bit
,
1290 int_no
< dig_no
|| i
< round_idx
);
1293 else if (dig_no
== int_no
)
1295 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1296 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1298 if (target_bit
== is_bit
)
1300 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1301 numsize
* sizeof (mp_limb_t
));
1302 /* FIXME: the following loop can be avoided if we assume a
1303 maximal MANT_DIG value. */
1304 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1306 else if (target_bit
> is_bit
)
1308 (void) mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1309 num
, numsize
, target_bit
- is_bit
);
1310 /* FIXME: the following loop can be avoided if we assume a
1311 maximal MANT_DIG value. */
1312 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1317 assert (numsize
< RETURN_LIMB_SIZE
);
1319 cy
= mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1320 num
, numsize
, is_bit
- target_bit
);
1321 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1322 /* FIXME: the following loop can be avoided if we assume a
1323 maximal MANT_DIG value. */
1324 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1327 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1331 /* Store the bits we already have. */
1332 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1333 #if RETURN_LIMB_SIZE > 1
1334 if (numsize
< RETURN_LIMB_SIZE
)
1335 # if RETURN_LIMB_SIZE == 2
1336 retval
[numsize
] = 0;
1338 MPN_ZERO (retval
+ numsize
, RETURN_LIMB_SIZE
- numsize
);
1343 /* We have to compute at least some of the fractional digits. */
1345 /* We construct a fraction and the result of the division gives us
1346 the needed digits. The denominator is 1.0 multiplied by the
1347 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1348 123e-6 gives 123 / 1000000. */
1353 int need_frac_digits
;
1355 mp_limb_t
*psrc
= den
;
1356 mp_limb_t
*pdest
= num
;
1357 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1359 assert (dig_no
> int_no
1361 && exponent
>= MIN_10_EXP
- (DIG
+ 1));
1363 /* We need to compute MANT_DIG - BITS fractional bits that lie
1364 within the mantissa of the result, the following bit for
1365 rounding, and to know whether any subsequent bit is 0.
1366 Computing a bit with value 2^-n means looking at n digits after
1367 the decimal point. */
1370 /* The bits required are those immediately after the point. */
1371 assert (int_no
> 0 && exponent
== 0);
1372 need_frac_digits
= 1 + MANT_DIG
- bits
;
1376 /* The number is in the form .123eEXPONENT. */
1377 assert (int_no
== 0 && *startp
!= L_('0'));
1378 /* The number is at least 10^(EXPONENT-1), and 10^3 <
1380 int neg_exp_2
= ((1 - exponent
) * 10) / 3 + 1;
1381 /* The number is at least 2^-NEG_EXP_2. We need up to
1382 MANT_DIG bits following that bit. */
1383 need_frac_digits
= neg_exp_2
+ MANT_DIG
;
1384 /* However, we never need bits beyond 1/4 ulp of the smallest
1385 representable value. (That 1/4 ulp bit is only needed to
1386 determine tinyness on machines where tinyness is determined
1388 if (need_frac_digits
> MANT_DIG
- MIN_EXP
+ 2)
1389 need_frac_digits
= MANT_DIG
- MIN_EXP
+ 2;
1390 /* At this point, NEED_FRAC_DIGITS is the total number of
1391 digits needed after the point, but some of those may be
1393 need_frac_digits
+= exponent
;
1394 /* Any cases underflowing enough that none of the fractional
1395 digits are needed should have been caught earlier (such
1396 cases are on the order of 10^-n or smaller where 2^-n is
1397 the least subnormal). */
1398 assert (need_frac_digits
> 0);
1401 if (need_frac_digits
> (intmax_t) dig_no
- (intmax_t) int_no
)
1402 need_frac_digits
= (intmax_t) dig_no
- (intmax_t) int_no
;
1404 if ((intmax_t) dig_no
> (intmax_t) int_no
+ need_frac_digits
)
1406 dig_no
= int_no
+ need_frac_digits
;
1412 neg_exp
= (intmax_t) dig_no
- (intmax_t) int_no
- exponent
;
1414 /* Construct the denominator. */
1419 if ((neg_exp
& expbit
) != 0)
1426 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1427 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1428 densize
* sizeof (mp_limb_t
));
1432 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1433 + _FPIO_CONST_OFFSET
],
1434 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1436 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1439 (void) SWAP (psrc
, pdest
);
1445 while (neg_exp
!= 0);
1448 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1450 /* Read the fractional digits from the string. */
1451 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1452 #ifndef USE_WIDE_CHAR
1453 , decimal
, decimal_len
, thousands
1457 /* We now have to shift both numbers so that the highest bit in the
1458 denominator is set. In the same process we copy the numerator to
1459 a high place in the array so that the division constructs the wanted
1460 digits. This is done by a "quasi fix point" number representation.
1462 num: ddddddddddd . 0000000000000000000000
1464 den: ddddddddddd n >= m
1468 count_leading_zeros (cnt
, den
[densize
- 1]);
1472 /* Don't call `mpn_shift' with a count of zero since the specification
1473 does not allow this. */
1474 (void) mpn_lshift (den
, den
, densize
, cnt
);
1475 cy
= mpn_lshift (num
, num
, numsize
, cnt
);
1477 num
[numsize
++] = cy
;
1480 /* Now we are ready for the division. But it is not necessary to
1481 do a full multi-precision division because we only need a small
1482 number of bits for the result. So we do not use mpn_divmod
1483 here but instead do the division here by hand and stop whenever
1484 the needed number of bits is reached. The code itself comes
1485 from the GNU MP Library by Torbj\"orn Granlund. */
1493 mp_limb_t d
, n
, quot
;
1498 assert (numsize
== 1 && n
< d
);
1502 udiv_qrnnd (quot
, n
, n
, 0, d
);
1509 cnt = BITS_PER_MP_LIMB; \
1511 count_leading_zeros (cnt, quot); \
1513 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1515 used = MANT_DIG + cnt; \
1516 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1517 bits = MANT_DIG + 1; \
1521 /* Note that we only clear the second element. */ \
1522 /* The conditional is determined at compile time. */ \
1523 if (RETURN_LIMB_SIZE > 1) \
1529 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1530 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1534 used = MANT_DIG - bits; \
1536 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1538 bits += BITS_PER_MP_LIMB
1542 while (bits
<= MANT_DIG
);
1544 return round_and_return (retval
, exponent
- 1, negative
,
1545 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1546 more_bits
|| n
!= 0);
1550 mp_limb_t d0
, d1
, n0
, n1
;
1557 if (numsize
< densize
)
1561 /* The numerator of the number occupies fewer bits than
1562 the denominator but the one limb is bigger than the
1563 high limb of the numerator. */
1570 exponent
-= BITS_PER_MP_LIMB
;
1573 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1574 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1575 BITS_PER_MP_LIMB
, 0);
1578 used
= MANT_DIG
- bits
;
1580 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1582 bits
+= BITS_PER_MP_LIMB
;
1594 while (bits
<= MANT_DIG
)
1600 /* QUOT should be either 111..111 or 111..110. We need
1601 special treatment of this rare case as normal division
1602 would give overflow. */
1603 quot
= ~(mp_limb_t
) 0;
1606 if (r
< d1
) /* Carry in the addition? */
1608 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1611 n1
= d0
- (d0
!= 0);
1616 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1617 umul_ppmm (n1
, n0
, d0
, quot
);
1621 if (n1
> r
|| (n1
== r
&& n0
> 0))
1623 /* The estimated QUOT was too large. */
1626 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1628 if (r
>= d1
) /* If not carry, test QUOT again. */
1631 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1637 return round_and_return (retval
, exponent
- 1, negative
,
1638 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1639 more_bits
|| n1
!= 0 || n0
!= 0);
1644 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1648 dX
= den
[densize
- 1];
1649 d1
= den
[densize
- 2];
1651 /* The division does not work if the upper limb of the two-limb
1652 numerator is greater than the denominator. */
1653 if (mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1656 if (numsize
< densize
)
1658 mp_size_t empty
= densize
- numsize
;
1662 exponent
-= empty
* BITS_PER_MP_LIMB
;
1665 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1667 /* We make a difference here because the compiler
1668 cannot optimize the `else' case that good and
1669 this reflects all currently used FLOAT types
1670 and GMP implementations. */
1671 #if RETURN_LIMB_SIZE <= 2
1672 assert (empty
== 1);
1673 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1674 BITS_PER_MP_LIMB
, 0);
1676 for (i
= RETURN_LIMB_SIZE
- 1; i
>= empty
; --i
)
1677 retval
[i
] = retval
[i
- empty
];
1684 used
= MANT_DIG
- bits
;
1685 if (used
>= BITS_PER_MP_LIMB
)
1688 (void) mpn_lshift (&retval
[used
1689 / BITS_PER_MP_LIMB
],
1692 - used
/ BITS_PER_MP_LIMB
),
1693 used
% BITS_PER_MP_LIMB
);
1694 for (i
= used
/ BITS_PER_MP_LIMB
- 1; i
>= 0; --i
)
1698 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1700 bits
+= empty
* BITS_PER_MP_LIMB
;
1702 for (i
= numsize
; i
> 0; --i
)
1703 num
[i
+ empty
] = num
[i
- 1];
1704 MPN_ZERO (num
, empty
+ 1);
1709 assert (numsize
== densize
);
1710 for (i
= numsize
; i
> 0; --i
)
1711 num
[i
] = num
[i
- 1];
1718 while (bits
<= MANT_DIG
)
1721 /* This might over-estimate QUOT, but it's probably not
1722 worth the extra code here to find out. */
1723 quot
= ~(mp_limb_t
) 0;
1728 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1729 umul_ppmm (n1
, n0
, d1
, quot
);
1731 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1735 if (r
< dX
) /* I.e. "carry in previous addition?" */
1742 /* Possible optimization: We already have (q * n0) and (1 * n1)
1743 after the calculation of QUOT. Taking advantage of this, we
1744 could make this loop make two iterations less. */
1746 cy
= mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1748 if (num
[densize
] != cy
)
1750 cy
= mpn_add_n (num
, num
, den
, densize
);
1754 n0
= num
[densize
] = num
[densize
- 1];
1755 for (i
= densize
- 1; i
> 0; --i
)
1756 num
[i
] = num
[i
- 1];
1762 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1764 return round_and_return (retval
, exponent
- 1, negative
,
1765 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1766 more_bits
|| i
>= 0);