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