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
; ++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 Implemented as a macro, so that __builtin_constant_p works even at -O0.
426 Tege doesn't like this macro so I have to write it here myself. :)
428 #define mpn_lshift_1(ptr, size, count, limb) \
431 mp_limb_t *__ptr = (ptr); \
432 if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \
435 for (i = (size) - 1; i > 0; --i) \
436 __ptr[i] = __ptr[i - 1]; \
441 /* We assume count > 0 && count < BITS_PER_MP_LIMB here. */ \
442 unsigned int __count = (count); \
443 (void) mpn_lshift (__ptr, __ptr, size, __count); \
444 __ptr[0] |= (limb) >> (BITS_PER_MP_LIMB - __count); \
450 #define INTERNAL(x) INTERNAL1(x)
451 #define INTERNAL1(x) __##x##_internal
452 #ifndef ____STRTOF_INTERNAL
453 # define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
456 /* This file defines a function to check for correct grouping. */
457 #include "grouping.h"
460 /* Return a floating point number with the value of the given string NPTR.
461 Set *ENDPTR to the character after the last used one. If the number is
462 smaller than the smallest representable number, set `errno' to ERANGE and
463 return 0.0. If the number is too big to be represented, set `errno' to
464 ERANGE and return HUGE_VAL with the appropriate sign. */
466 ____STRTOF_INTERNAL (nptr
, endptr
, group
)
467 const STRING_TYPE
*nptr
;
468 STRING_TYPE
**endptr
;
471 int negative
; /* The sign of the number. */
472 MPN_VAR (num
); /* MP representation of the number. */
473 intmax_t exponent
; /* Exponent of the number. */
475 /* Numbers starting `0X' or `0x' have to be processed with base 16. */
478 /* When we have to compute fractional digits we form a fraction with a
479 second multi-precision number (and we sometimes need a second for
480 temporary results). */
483 /* Representation for the return value. */
484 mp_limb_t retval
[RETURN_LIMB_SIZE
];
485 /* Number of bits currently in result value. */
488 /* Running pointer after the last character processed in the string. */
489 const STRING_TYPE
*cp
, *tp
;
490 /* Start of significant part of the number. */
491 const STRING_TYPE
*startp
, *start_of_digits
;
492 /* Points at the character following the integer and fractional digits. */
493 const STRING_TYPE
*expp
;
494 /* Total number of digit and number of digits in integer part. */
495 size_t dig_no
, int_no
, lead_zero
;
496 /* Contains the last character read. */
499 /* We should get wint_t from <stddef.h>, but not all GCC versions define it
500 there. So define it ourselves if it remains undefined. */
502 typedef unsigned int wint_t;
504 /* The radix character of the current locale. */
511 /* The thousands character of the current locale. */
513 wchar_t thousands
= L
'\0';
515 const char *thousands
= NULL
;
517 /* The numeric grouping specification of the current locale,
518 in the format described in <locale.h>. */
519 const char *grouping
;
520 /* Used in several places. */
523 #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
524 const struct lconv
*lc
= localeconv ();
527 if (__builtin_expect (group
, 0))
529 #ifdef USE_NL_LANGINFO
530 grouping
= nl_langinfo (GROUPING
);
531 if (*grouping
<= 0 || *grouping
== CHAR_MAX
)
535 /* Figure out the thousands separator character. */
537 thousands
= nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC
);
538 if (thousands
== L
'\0')
541 thousands
= nl_langinfo (THOUSANDS_SEP
);
542 if (*thousands
== '\0')
549 #elif defined USE_LOCALECONV
550 grouping
= lc
->grouping
;
551 if (grouping
== NULL
|| *grouping
<= 0 || *grouping
== CHAR_MAX
)
555 /* Figure out the thousands separator character. */
556 thousands
= lc
->thousands_sep
;
557 if (thousands
== NULL
|| *thousands
== '\0')
570 /* Find the locale's decimal point character. */
572 decimal
= nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC
);
573 assert (decimal
!= L
'\0');
574 # define decimal_len 1
576 #ifdef USE_NL_LANGINFO
577 decimal
= nl_langinfo (DECIMAL_POINT
);
578 decimal_len
= strlen (decimal
);
579 assert (decimal_len
> 0);
580 #elif defined USE_LOCALECONV
581 decimal
= lc
->decimal_point
;
582 if (decimal
== NULL
|| *decimal
== '\0')
584 decimal_len
= strlen (decimal
);
591 /* Prepare number representation. */
596 /* Parse string to get maximal legal prefix. We need the number of
597 characters of the integer part, the fractional part and the exponent. */
599 /* Ignore leading white space. */
604 /* Get sign of the result. */
610 else if (c
== L_('+'))
613 /* Return 0.0 if no legal string is found.
614 No character is used even if a sign was found. */
616 if (c
== (wint_t) decimal
617 && (wint_t) cp
[1] >= L
'0' && (wint_t) cp
[1] <= L
'9')
619 /* We accept it. This funny construct is here only to indent
620 the code correctly. */
623 for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
624 if (cp
[cnt
] != decimal
[cnt
])
626 if (decimal
[cnt
] == '\0' && cp
[cnt
] >= '0' && cp
[cnt
] <= '9')
628 /* We accept it. This funny construct is here only to indent
629 the code correctly. */
632 else if (c
< L_('0') || c
> L_('9'))
634 /* Check for `INF' or `INFINITY'. */
635 CHAR_TYPE lowc
= TOLOWER_C (c
);
637 if (lowc
== L_('i') && STRNCASECMP (cp
, L_("inf"), 3) == 0)
639 /* Return +/- infinity. */
641 *endptr
= (STRING_TYPE
*)
642 (cp
+ (STRNCASECMP (cp
+ 3, L_("inity"), 5) == 0
645 return negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
;
648 if (lowc
== L_('n') && STRNCASECMP (cp
, L_("nan"), 3) == 0)
655 /* Match `(n-char-sequence-digit)'. */
658 const STRING_TYPE
*startp
= cp
;
661 while ((*cp
>= L_('0') && *cp
<= L_('9'))
662 || ({ CHAR_TYPE lo
= TOLOWER (*cp
);
663 lo
>= L_('a') && lo
<= L_('z'); })
667 /* The closing brace is missing. Only match the NAN
672 /* This is a system-dependent way to specify the
673 bitmask used for the NaN. We expect it to be
674 a number which is put in the mantissa of the
677 unsigned long long int mant
;
679 mant
= STRTOULL (startp
+ 1, &endp
, 0);
681 SET_MANTISSA (retval
, mant
);
683 /* Consume the closing brace. */
689 *endptr
= (STRING_TYPE
*) cp
;
694 /* It is really a text we do not recognize. */
698 /* First look whether we are faced with a hexadecimal number. */
699 if (c
== L_('0') && TOLOWER (cp
[1]) == L_('x'))
701 /* Okay, it is a hexa-decimal number. Remember this and skip
702 the characters. BTW: hexadecimal numbers must not be
710 /* Record the start of the digits, in case we will check their grouping. */
711 start_of_digits
= startp
= cp
;
713 /* Ignore leading zeroes. This helps us to avoid useless computations. */
715 while (c
== L
'0' || ((wint_t) thousands
!= L
'\0' && c
== (wint_t) thousands
))
718 if (__builtin_expect (thousands
== NULL
, 1))
723 /* We also have the multibyte thousands string. */
728 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
729 if (thousands
[cnt
] != cp
[cnt
])
731 if (thousands
[cnt
] != '\0')
740 /* If no other digit but a '0' is found the result is 0.0.
741 Return current read pointer. */
742 CHAR_TYPE lowc
= TOLOWER (c
);
743 if (!((c
>= L_('0') && c
<= L_('9'))
744 || (base
== 16 && lowc
>= L_('a') && lowc
<= L_('f'))
747 c
== (wint_t) decimal
749 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
750 if (decimal
[cnt
] != cp
[cnt
])
752 decimal
[cnt
] == '\0'; })
754 /* '0x.' alone is not a valid hexadecimal number.
755 '.' alone is not valid either, but that has been checked
758 || cp
!= start_of_digits
759 || (cp
[decimal_len
] >= L_('0') && cp
[decimal_len
] <= L_('9'))
760 || ({ CHAR_TYPE lo
= TOLOWER (cp
[decimal_len
]);
761 lo
>= L_('a') && lo
<= L_('f'); })))
762 || (base
== 16 && (cp
!= start_of_digits
764 || (base
!= 16 && lowc
== L_('e'))))
767 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
770 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
773 /* If TP is at the start of the digits, there was no correctly
774 grouped prefix of the string; so no number found. */
775 RETURN (negative
? -0.0 : 0.0,
776 tp
== start_of_digits
? (base
== 16 ? cp
- 1 : nptr
) : tp
);
779 /* Remember first significant digit and read following characters until the
780 decimal point, exponent character or any non-FP number character. */
785 if ((c
>= L_('0') && c
<= L_('9'))
787 && ({ CHAR_TYPE lo
= TOLOWER (c
);
788 lo
>= L_('a') && lo
<= L_('f'); })))
793 if (__builtin_expect ((wint_t) thousands
== L
'\0', 1)
794 || c
!= (wint_t) thousands
)
795 /* Not a digit or separator: end of the integer part. */
798 if (__builtin_expect (thousands
== NULL
, 1))
802 for (cnt
= 0; thousands
[cnt
] != '\0'; ++cnt
)
803 if (thousands
[cnt
] != cp
[cnt
])
805 if (thousands
[cnt
] != '\0')
814 if (__builtin_expect (grouping
!= NULL
, 0) && cp
> start_of_digits
)
816 /* Check the grouping of the digits. */
818 tp
= __correctly_grouped_prefixwc (start_of_digits
, cp
, thousands
,
821 tp
= __correctly_grouped_prefixmb (start_of_digits
, cp
, thousands
,
826 /* Less than the entire string was correctly grouped. */
828 if (tp
== start_of_digits
)
829 /* No valid group of numbers at all: no valid number. */
833 /* The number is validly grouped, but consists
834 only of zeroes. The whole value is zero. */
835 RETURN (negative
? -0.0 : 0.0, tp
);
837 /* Recompute DIG_NO so we won't read more digits than
838 are properly grouped. */
841 for (tp
= startp
; tp
< cp
; ++tp
)
842 if (*tp
>= L_('0') && *tp
<= L_('9'))
852 /* We have the number of digits in the integer part. Whether these
853 are all or any is really a fractional digit will be decided
856 lead_zero
= int_no
== 0 ? (size_t) -1 : 0;
858 /* Read the fractional digits. A special case are the 'american
859 style' numbers like `16.' i.e. with decimal point but without
863 c
== (wint_t) decimal
865 ({ for (cnt
= 0; decimal
[cnt
] != '\0'; ++cnt
)
866 if (decimal
[cnt
] != cp
[cnt
])
868 decimal
[cnt
] == '\0'; })
874 while ((c
>= L_('0') && c
<= L_('9')) ||
875 (base
== 16 && ({ CHAR_TYPE lo
= TOLOWER (c
);
876 lo
>= L_('a') && lo
<= L_('f'); })))
878 if (c
!= L_('0') && lead_zero
== (size_t) -1)
879 lead_zero
= dig_no
- int_no
;
884 assert (dig_no
<= (uintmax_t) INTMAX_MAX
);
886 /* Remember start of exponent (if any). */
891 if ((base
== 16 && lowc
== L_('p'))
892 || (base
!= 16 && lowc
== L_('e')))
894 int exp_negative
= 0;
902 else if (c
== L_('+'))
905 if (c
>= L_('0') && c
<= L_('9'))
909 /* Get the exponent limit. */
914 assert (int_no
<= (uintmax_t) (INTMAX_MAX
915 + MIN_EXP
- MANT_DIG
) / 4);
916 exp_limit
= -MIN_EXP
+ MANT_DIG
+ 4 * (intmax_t) int_no
;
922 assert (lead_zero
== 0
923 && int_no
<= (uintmax_t) INTMAX_MAX
/ 4);
924 exp_limit
= MAX_EXP
- 4 * (intmax_t) int_no
+ 3;
926 else if (lead_zero
== (size_t) -1)
928 /* The number is zero and this limit is
930 exp_limit
= MAX_EXP
+ 3;
935 <= (uintmax_t) (INTMAX_MAX
- MAX_EXP
- 3) / 4);
937 + 4 * (intmax_t) lead_zero
947 <= (uintmax_t) (INTMAX_MAX
+ MIN_10_EXP
- MANT_DIG
));
948 exp_limit
= -MIN_10_EXP
+ MANT_DIG
+ (intmax_t) int_no
;
954 assert (lead_zero
== 0
955 && int_no
<= (uintmax_t) INTMAX_MAX
);
956 exp_limit
= MAX_10_EXP
- (intmax_t) int_no
+ 1;
958 else if (lead_zero
== (size_t) -1)
960 /* The number is zero and this limit is
962 exp_limit
= MAX_10_EXP
+ 1;
967 <= (uintmax_t) (INTMAX_MAX
- MAX_10_EXP
- 1));
968 exp_limit
= MAX_10_EXP
+ (intmax_t) lead_zero
+ 1;
978 if (__builtin_expect ((exponent
> exp_limit
/ 10
979 || (exponent
== exp_limit
/ 10
980 && c
- L_('0') > exp_limit
% 10)), 0))
981 /* The exponent is too large/small to represent a valid
986 /* We have to take care for special situation: a joker
987 might have written "0.0e100000" which is in fact
989 if (lead_zero
== (size_t) -1)
990 result
= negative
? -0.0 : 0.0;
993 /* Overflow or underflow. */
994 #if defined HAVE_ERRNO_H && defined ERANGE
997 result
= (exp_negative
? (negative
? -0.0 : 0.0) :
998 negative
? -FLOAT_HUGE_VAL
: FLOAT_HUGE_VAL
);
1001 /* Accept all following digits as part of the exponent. */
1004 while (*cp
>= L_('0') && *cp
<= L_('9'));
1006 RETURN (result
, cp
);
1011 exponent
+= c
- L_('0');
1015 while (c
>= L_('0') && c
<= L_('9'));
1018 exponent
= -exponent
;
1024 /* We don't want to have to work with trailing zeroes after the radix. */
1025 if (dig_no
> int_no
)
1027 while (expp
[-1] == L_('0'))
1032 assert (dig_no
>= int_no
);
1035 if (dig_no
== int_no
&& dig_no
> 0 && exponent
< 0)
1038 while (! (base
== 16 ? ISXDIGIT (expp
[-1]) : ISDIGIT (expp
[-1])))
1041 if (expp
[-1] != L_('0'))
1047 exponent
+= base
== 16 ? 4 : 1;
1049 while (dig_no
> 0 && exponent
< 0);
1053 /* The whole string is parsed. Store the address of the next character. */
1055 *endptr
= (STRING_TYPE
*) cp
;
1058 return negative
? -0.0 : 0.0;
1062 /* Find the decimal point */
1063 #ifdef USE_WIDE_CHAR
1064 while (*startp
!= decimal
)
1069 if (*startp
== decimal
[0])
1071 for (cnt
= 1; decimal
[cnt
] != '\0'; ++cnt
)
1072 if (decimal
[cnt
] != startp
[cnt
])
1074 if (decimal
[cnt
] == '\0')
1080 startp
+= lead_zero
+ decimal_len
;
1081 assert (lead_zero
<= (base
== 16
1082 ? (uintmax_t) INTMAX_MAX
/ 4
1083 : (uintmax_t) INTMAX_MAX
));
1084 assert (lead_zero
<= (base
== 16
1085 ? ((uintmax_t) exponent
1086 - (uintmax_t) INTMAX_MIN
) / 4
1087 : ((uintmax_t) exponent
- (uintmax_t) INTMAX_MIN
)));
1088 exponent
-= base
== 16 ? 4 * (intmax_t) lead_zero
: (intmax_t) lead_zero
;
1089 dig_no
-= lead_zero
;
1092 /* If the BASE is 16 we can use a simpler algorithm. */
1095 static const int nbits
[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
1096 4, 4, 4, 4, 4, 4, 4, 4 };
1097 int idx
= (MANT_DIG
- 1) / BITS_PER_MP_LIMB
;
1098 int pos
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1101 while (!ISXDIGIT (*startp
))
1103 while (*startp
== L_('0'))
1105 if (ISDIGIT (*startp
))
1106 val
= *startp
++ - L_('0');
1108 val
= 10 + TOLOWER (*startp
++) - L_('a');
1110 /* We cannot have a leading zero. */
1113 if (pos
+ 1 >= 4 || pos
+ 1 >= bits
)
1115 /* We don't have to care for wrapping. This is the normal
1116 case so we add the first clause in the `if' expression as
1117 an optimization. It is a compile-time constant and so does
1118 not cost anything. */
1119 retval
[idx
] = val
<< (pos
- bits
+ 1);
1124 retval
[idx
--] = val
>> (bits
- pos
- 1);
1125 retval
[idx
] = val
<< (BITS_PER_MP_LIMB
- (bits
- pos
- 1));
1126 pos
= BITS_PER_MP_LIMB
- 1 - (bits
- pos
- 1);
1129 /* Adjust the exponent for the bits we are shifting in. */
1130 assert (int_no
<= (uintmax_t) (exponent
< 0
1131 ? (INTMAX_MAX
- bits
+ 1) / 4
1132 : (INTMAX_MAX
- exponent
- bits
+ 1) / 4));
1133 exponent
+= bits
- 1 + ((intmax_t) int_no
- 1) * 4;
1135 while (--dig_no
> 0 && idx
>= 0)
1137 if (!ISXDIGIT (*startp
))
1138 startp
+= decimal_len
;
1139 if (ISDIGIT (*startp
))
1140 val
= *startp
++ - L_('0');
1142 val
= 10 + TOLOWER (*startp
++) - L_('a');
1146 retval
[idx
] |= val
<< (pos
- 4 + 1);
1151 retval
[idx
--] |= val
>> (4 - pos
- 1);
1152 val
<<= BITS_PER_MP_LIMB
- (4 - pos
- 1);
1155 int rest_nonzero
= 0;
1156 while (--dig_no
> 0)
1158 if (*startp
!= L_('0'))
1165 return round_and_return (retval
, exponent
, negative
, val
,
1166 BITS_PER_MP_LIMB
- 1, rest_nonzero
);
1170 pos
= BITS_PER_MP_LIMB
- 1 - (4 - pos
- 1);
1174 /* We ran out of digits. */
1175 MPN_ZERO (retval
, idx
);
1177 return round_and_return (retval
, exponent
, negative
, 0, 0, 0);
1180 /* Now we have the number of digits in total and the integer digits as well
1181 as the exponent and its sign. We can decide whether the read digits are
1182 really integer digits or belong to the fractional part; i.e. we normalize
1185 register intmax_t incr
= (exponent
< 0
1186 ? MAX (-(intmax_t) int_no
, exponent
)
1187 : MIN ((intmax_t) dig_no
- (intmax_t) int_no
,
1193 if (__builtin_expect (exponent
> MAX_10_EXP
+ 1 - (intmax_t) int_no
, 0))
1194 return overflow_value (negative
);
1196 if (__builtin_expect (exponent
< MIN_10_EXP
- (DIG
+ 1), 0))
1197 return underflow_value (negative
);
1201 /* Read the integer part as a multi-precision number to NUM. */
1202 startp
= str_to_mpn (startp
, int_no
, num
, &numsize
, &exponent
1203 #ifndef USE_WIDE_CHAR
1204 , decimal
, decimal_len
, thousands
1210 /* We now multiply the gained number by the given power of ten. */
1211 mp_limb_t
*psrc
= num
;
1212 mp_limb_t
*pdest
= den
;
1214 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1218 if ((exponent
& expbit
) != 0)
1220 size_t size
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1224 /* FIXME: not the whole multiplication has to be
1225 done. If we have the needed number of bits we
1226 only need the information whether more non-zero
1228 if (numsize
>= ttab
->arraysize
- _FPIO_CONST_OFFSET
)
1229 cy
= mpn_mul (pdest
, psrc
, numsize
,
1230 &__tens
[ttab
->arrayoff
1231 + _FPIO_CONST_OFFSET
],
1234 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1235 + _FPIO_CONST_OFFSET
],
1236 size
, psrc
, numsize
);
1240 (void) SWAP (psrc
, pdest
);
1245 while (exponent
!= 0);
1248 memcpy (num
, den
, numsize
* sizeof (mp_limb_t
));
1251 /* Determine how many bits of the result we already have. */
1252 count_leading_zeros (bits
, num
[numsize
- 1]);
1253 bits
= numsize
* BITS_PER_MP_LIMB
- bits
;
1255 /* Now we know the exponent of the number in base two.
1256 Check it against the maximum possible exponent. */
1257 if (__builtin_expect (bits
> MAX_EXP
, 0))
1258 return overflow_value (negative
);
1260 /* We have already the first BITS bits of the result. Together with
1261 the information whether more non-zero bits follow this is enough
1262 to determine the result. */
1263 if (bits
> MANT_DIG
)
1266 const mp_size_t least_idx
= (bits
- MANT_DIG
) / BITS_PER_MP_LIMB
;
1267 const mp_size_t least_bit
= (bits
- MANT_DIG
) % BITS_PER_MP_LIMB
;
1268 const mp_size_t round_idx
= least_bit
== 0 ? least_idx
- 1
1270 const mp_size_t round_bit
= least_bit
== 0 ? BITS_PER_MP_LIMB
- 1
1274 memcpy (retval
, &num
[least_idx
],
1275 RETURN_LIMB_SIZE
* sizeof (mp_limb_t
));
1278 for (i
= least_idx
; i
< numsize
- 1; ++i
)
1279 retval
[i
- least_idx
] = (num
[i
] >> least_bit
)
1281 << (BITS_PER_MP_LIMB
- least_bit
));
1282 if (i
- least_idx
< RETURN_LIMB_SIZE
)
1283 retval
[RETURN_LIMB_SIZE
- 1] = num
[i
] >> least_bit
;
1286 /* Check whether any limb beside the ones in RETVAL are non-zero. */
1287 for (i
= 0; num
[i
] == 0; ++i
)
1290 return round_and_return (retval
, bits
- 1, negative
,
1291 num
[round_idx
], round_bit
,
1292 int_no
< dig_no
|| i
< round_idx
);
1295 else if (dig_no
== int_no
)
1297 const mp_size_t target_bit
= (MANT_DIG
- 1) % BITS_PER_MP_LIMB
;
1298 const mp_size_t is_bit
= (bits
- 1) % BITS_PER_MP_LIMB
;
1300 if (target_bit
== is_bit
)
1302 memcpy (&retval
[RETURN_LIMB_SIZE
- numsize
], num
,
1303 numsize
* sizeof (mp_limb_t
));
1304 /* FIXME: the following loop can be avoided if we assume a
1305 maximal MANT_DIG value. */
1306 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1308 else if (target_bit
> is_bit
)
1310 (void) mpn_lshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1311 num
, numsize
, target_bit
- is_bit
);
1312 /* FIXME: the following loop can be avoided if we assume a
1313 maximal MANT_DIG value. */
1314 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
);
1319 assert (numsize
< RETURN_LIMB_SIZE
);
1321 cy
= mpn_rshift (&retval
[RETURN_LIMB_SIZE
- numsize
],
1322 num
, numsize
, is_bit
- target_bit
);
1323 retval
[RETURN_LIMB_SIZE
- numsize
- 1] = cy
;
1324 /* FIXME: the following loop can be avoided if we assume a
1325 maximal MANT_DIG value. */
1326 MPN_ZERO (retval
, RETURN_LIMB_SIZE
- numsize
- 1);
1329 return round_and_return (retval
, bits
- 1, negative
, 0, 0, 0);
1333 /* Store the bits we already have. */
1334 memcpy (retval
, num
, numsize
* sizeof (mp_limb_t
));
1335 #if RETURN_LIMB_SIZE > 1
1336 if (numsize
< RETURN_LIMB_SIZE
)
1337 # if RETURN_LIMB_SIZE == 2
1338 retval
[numsize
] = 0;
1340 MPN_ZERO (retval
+ numsize
, RETURN_LIMB_SIZE
- numsize
);
1345 /* We have to compute at least some of the fractional digits. */
1347 /* We construct a fraction and the result of the division gives us
1348 the needed digits. The denominator is 1.0 multiplied by the
1349 exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
1350 123e-6 gives 123 / 1000000. */
1355 int need_frac_digits
;
1357 mp_limb_t
*psrc
= den
;
1358 mp_limb_t
*pdest
= num
;
1359 const struct mp_power
*ttab
= &_fpioconst_pow10
[0];
1361 assert (dig_no
> int_no
1363 && exponent
>= MIN_10_EXP
- (DIG
+ 1));
1365 /* We need to compute MANT_DIG - BITS fractional bits that lie
1366 within the mantissa of the result, the following bit for
1367 rounding, and to know whether any subsequent bit is 0.
1368 Computing a bit with value 2^-n means looking at n digits after
1369 the decimal point. */
1372 /* The bits required are those immediately after the point. */
1373 assert (int_no
> 0 && exponent
== 0);
1374 need_frac_digits
= 1 + MANT_DIG
- bits
;
1378 /* The number is in the form .123eEXPONENT. */
1379 assert (int_no
== 0 && *startp
!= L_('0'));
1380 /* The number is at least 10^(EXPONENT-1), and 10^3 <
1382 int neg_exp_2
= ((1 - exponent
) * 10) / 3 + 1;
1383 /* The number is at least 2^-NEG_EXP_2. We need up to
1384 MANT_DIG bits following that bit. */
1385 need_frac_digits
= neg_exp_2
+ MANT_DIG
;
1386 /* However, we never need bits beyond 1/4 ulp of the smallest
1387 representable value. (That 1/4 ulp bit is only needed to
1388 determine tinyness on machines where tinyness is determined
1390 if (need_frac_digits
> MANT_DIG
- MIN_EXP
+ 2)
1391 need_frac_digits
= MANT_DIG
- MIN_EXP
+ 2;
1392 /* At this point, NEED_FRAC_DIGITS is the total number of
1393 digits needed after the point, but some of those may be
1395 need_frac_digits
+= exponent
;
1396 /* Any cases underflowing enough that none of the fractional
1397 digits are needed should have been caught earlier (such
1398 cases are on the order of 10^-n or smaller where 2^-n is
1399 the least subnormal). */
1400 assert (need_frac_digits
> 0);
1403 if (need_frac_digits
> (intmax_t) dig_no
- (intmax_t) int_no
)
1404 need_frac_digits
= (intmax_t) dig_no
- (intmax_t) int_no
;
1406 if ((intmax_t) dig_no
> (intmax_t) int_no
+ need_frac_digits
)
1408 dig_no
= int_no
+ need_frac_digits
;
1414 neg_exp
= (intmax_t) dig_no
- (intmax_t) int_no
- exponent
;
1416 /* Construct the denominator. */
1421 if ((neg_exp
& expbit
) != 0)
1428 densize
= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1429 memcpy (psrc
, &__tens
[ttab
->arrayoff
+ _FPIO_CONST_OFFSET
],
1430 densize
* sizeof (mp_limb_t
));
1434 cy
= mpn_mul (pdest
, &__tens
[ttab
->arrayoff
1435 + _FPIO_CONST_OFFSET
],
1436 ttab
->arraysize
- _FPIO_CONST_OFFSET
,
1438 densize
+= ttab
->arraysize
- _FPIO_CONST_OFFSET
;
1441 (void) SWAP (psrc
, pdest
);
1447 while (neg_exp
!= 0);
1450 memcpy (den
, num
, densize
* sizeof (mp_limb_t
));
1452 /* Read the fractional digits from the string. */
1453 (void) str_to_mpn (startp
, dig_no
- int_no
, num
, &numsize
, &exponent
1454 #ifndef USE_WIDE_CHAR
1455 , decimal
, decimal_len
, thousands
1459 /* We now have to shift both numbers so that the highest bit in the
1460 denominator is set. In the same process we copy the numerator to
1461 a high place in the array so that the division constructs the wanted
1462 digits. This is done by a "quasi fix point" number representation.
1464 num: ddddddddddd . 0000000000000000000000
1466 den: ddddddddddd n >= m
1470 count_leading_zeros (cnt
, den
[densize
- 1]);
1474 /* Don't call `mpn_shift' with a count of zero since the specification
1475 does not allow this. */
1476 (void) mpn_lshift (den
, den
, densize
, cnt
);
1477 cy
= mpn_lshift (num
, num
, numsize
, cnt
);
1479 num
[numsize
++] = cy
;
1482 /* Now we are ready for the division. But it is not necessary to
1483 do a full multi-precision division because we only need a small
1484 number of bits for the result. So we do not use mpn_divmod
1485 here but instead do the division here by hand and stop whenever
1486 the needed number of bits is reached. The code itself comes
1487 from the GNU MP Library by Torbj\"orn Granlund. */
1495 mp_limb_t d
, n
, quot
;
1500 assert (numsize
== 1 && n
< d
);
1504 udiv_qrnnd (quot
, n
, n
, 0, d
);
1511 cnt = BITS_PER_MP_LIMB; \
1513 count_leading_zeros (cnt, quot); \
1515 if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
1517 used = MANT_DIG + cnt; \
1518 retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
1519 bits = MANT_DIG + 1; \
1523 /* Note that we only clear the second element. */ \
1524 /* The conditional is determined at compile time. */ \
1525 if (RETURN_LIMB_SIZE > 1) \
1531 else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
1532 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
1536 used = MANT_DIG - bits; \
1538 mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
1540 bits += BITS_PER_MP_LIMB
1544 while (bits
<= MANT_DIG
);
1546 return round_and_return (retval
, exponent
- 1, negative
,
1547 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1548 more_bits
|| n
!= 0);
1552 mp_limb_t d0
, d1
, n0
, n1
;
1559 if (numsize
< densize
)
1563 /* The numerator of the number occupies fewer bits than
1564 the denominator but the one limb is bigger than the
1565 high limb of the numerator. */
1572 exponent
-= BITS_PER_MP_LIMB
;
1575 if (bits
+ BITS_PER_MP_LIMB
<= MANT_DIG
)
1576 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1577 BITS_PER_MP_LIMB
, 0);
1580 used
= MANT_DIG
- bits
;
1582 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1584 bits
+= BITS_PER_MP_LIMB
;
1596 while (bits
<= MANT_DIG
)
1602 /* QUOT should be either 111..111 or 111..110. We need
1603 special treatment of this rare case as normal division
1604 would give overflow. */
1605 quot
= ~(mp_limb_t
) 0;
1608 if (r
< d1
) /* Carry in the addition? */
1610 add_ssaaaa (n1
, n0
, r
- d0
, 0, 0, d0
);
1613 n1
= d0
- (d0
!= 0);
1618 udiv_qrnnd (quot
, r
, n1
, n0
, d1
);
1619 umul_ppmm (n1
, n0
, d0
, quot
);
1623 if (n1
> r
|| (n1
== r
&& n0
> 0))
1625 /* The estimated QUOT was too large. */
1628 sub_ddmmss (n1
, n0
, n1
, n0
, 0, d0
);
1630 if (r
>= d1
) /* If not carry, test QUOT again. */
1633 sub_ddmmss (n1
, n0
, r
, 0, n1
, n0
);
1639 return round_and_return (retval
, exponent
- 1, negative
,
1640 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1641 more_bits
|| n1
!= 0 || n0
!= 0);
1646 mp_limb_t cy
, dX
, d1
, n0
, n1
;
1650 dX
= den
[densize
- 1];
1651 d1
= den
[densize
- 2];
1653 /* The division does not work if the upper limb of the two-limb
1654 numerator is greater than the denominator. */
1655 if (mpn_cmp (num
, &den
[densize
- numsize
], numsize
) > 0)
1658 if (numsize
< densize
)
1660 mp_size_t empty
= densize
- numsize
;
1664 exponent
-= empty
* BITS_PER_MP_LIMB
;
1667 if (bits
+ empty
* BITS_PER_MP_LIMB
<= MANT_DIG
)
1669 /* We make a difference here because the compiler
1670 cannot optimize the `else' case that good and
1671 this reflects all currently used FLOAT types
1672 and GMP implementations. */
1673 #if RETURN_LIMB_SIZE <= 2
1674 assert (empty
== 1);
1675 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
,
1676 BITS_PER_MP_LIMB
, 0);
1678 for (i
= RETURN_LIMB_SIZE
- 1; i
>= empty
; --i
)
1679 retval
[i
] = retval
[i
- empty
];
1686 used
= MANT_DIG
- bits
;
1687 if (used
>= BITS_PER_MP_LIMB
)
1690 (void) mpn_lshift (&retval
[used
1691 / BITS_PER_MP_LIMB
],
1694 - used
/ BITS_PER_MP_LIMB
),
1695 used
% BITS_PER_MP_LIMB
);
1696 for (i
= used
/ BITS_PER_MP_LIMB
- 1; i
>= 0; --i
)
1700 mpn_lshift_1 (retval
, RETURN_LIMB_SIZE
, used
, 0);
1702 bits
+= empty
* BITS_PER_MP_LIMB
;
1704 for (i
= numsize
; i
> 0; --i
)
1705 num
[i
+ empty
] = num
[i
- 1];
1706 MPN_ZERO (num
, empty
+ 1);
1711 assert (numsize
== densize
);
1712 for (i
= numsize
; i
> 0; --i
)
1713 num
[i
] = num
[i
- 1];
1720 while (bits
<= MANT_DIG
)
1723 /* This might over-estimate QUOT, but it's probably not
1724 worth the extra code here to find out. */
1725 quot
= ~(mp_limb_t
) 0;
1730 udiv_qrnnd (quot
, r
, n0
, num
[densize
- 1], dX
);
1731 umul_ppmm (n1
, n0
, d1
, quot
);
1733 while (n1
> r
|| (n1
== r
&& n0
> num
[densize
- 2]))
1737 if (r
< dX
) /* I.e. "carry in previous addition?" */
1744 /* Possible optimization: We already have (q * n0) and (1 * n1)
1745 after the calculation of QUOT. Taking advantage of this, we
1746 could make this loop make two iterations less. */
1748 cy
= mpn_submul_1 (num
, den
, densize
+ 1, quot
);
1750 if (num
[densize
] != cy
)
1752 cy
= mpn_add_n (num
, num
, den
, densize
);
1756 n0
= num
[densize
] = num
[densize
- 1];
1757 for (i
= densize
- 1; i
> 0; --i
)
1758 num
[i
] = num
[i
- 1];
1764 for (i
= densize
; num
[i
] == 0 && i
>= 0; --i
)
1766 return round_and_return (retval
, exponent
- 1, negative
,
1767 quot
, BITS_PER_MP_LIMB
- 1 - used
,
1768 more_bits
|| i
>= 0);