]>
Commit | Line | Data |
---|---|---|
0501d603 | 1 | /* Convert string representing a number to float value, using given locale. |
ccadf7b5 | 2 | Copyright (C) 1997,98,2002, 2004 Free Software Foundation, Inc. |
0501d603 UD |
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 | |
41bdb6e2 AJ |
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. | |
0501d603 UD |
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 | |
41bdb6e2 | 14 | Lesser General Public License for more details. |
0501d603 | 15 | |
41bdb6e2 AJ |
16 | You should have received a copy of the GNU Lesser General Public |
17 | License along with the GNU C Library; if not, write to the Free | |
18 | Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA | |
19 | 02111-1307 USA. */ | |
0501d603 | 20 | |
ebbad4cc UD |
21 | #include <xlocale.h> |
22 | ||
0501d603 | 23 | extern double ____strtod_l_internal (const char *, char **, int, __locale_t); |
c3d6c951 UD |
24 | extern unsigned long long int ____strtoull_l_internal (const char *, char **, |
25 | int, int, __locale_t); | |
0501d603 | 26 | |
ccadf7b5 UD |
27 | /* Configuration part. These macros are defined by `strtold.c', |
28 | `strtof.c', `wcstod.c', `wcstold.c', and `wcstof.c' to produce the | |
29 | `long double' and `float' versions of the reader. */ | |
30 | #ifndef FLOAT | |
31 | # define FLOAT double | |
32 | # define FLT DBL | |
33 | # ifdef USE_WIDE_CHAR | |
34 | # define STRTOF wcstod_l | |
35 | # define __STRTOF __wcstod_l | |
36 | # else | |
37 | # define STRTOF strtod_l | |
38 | # define __STRTOF __strtod_l | |
39 | # endif | |
40 | # define MPN2FLOAT __mpn_construct_double | |
41 | # define FLOAT_HUGE_VAL HUGE_VAL | |
42 | # define SET_MANTISSA(flt, mant) \ | |
43 | do { union ieee754_double u; \ | |
44 | u.d = (flt); \ | |
45 | if ((mant & 0xfffffffffffffULL) == 0) \ | |
46 | mant = 0x8000000000000ULL; \ | |
47 | u.ieee.mantissa0 = ((mant) >> 32) & 0xfffff; \ | |
48 | u.ieee.mantissa1 = (mant) & 0xffffffff; \ | |
49 | (flt) = u.d; \ | |
50 | } while (0) | |
51 | #endif | |
52 | /* End of configuration part. */ | |
53 | \f | |
54 | #include <ctype.h> | |
55 | #include <errno.h> | |
56 | #include <float.h> | |
57 | #include <ieee754.h> | |
58 | #include "../locale/localeinfo.h" | |
59 | #include <locale.h> | |
60 | #include <math.h> | |
61 | #include <stdlib.h> | |
62 | #include <string.h> | |
63 | ||
64 | /* The gmp headers need some configuration frobs. */ | |
65 | #define HAVE_ALLOCA 1 | |
66 | ||
67 | /* Include gmp-mparam.h first, such that definitions of _SHORT_LIMB | |
68 | and _LONG_LONG_LIMB in it can take effect into gmp.h. */ | |
69 | #include <gmp-mparam.h> | |
70 | #include <gmp.h> | |
71 | #include <gmp-impl.h> | |
72 | #include <longlong.h> | |
73 | #include "fpioconst.h" | |
74 | ||
75 | #define NDEBUG 1 | |
76 | #include <assert.h> | |
77 | ||
78 | ||
79 | /* We use this code for the extended locale handling where the | |
80 | function gets as an additional argument the locale which has to be | |
81 | used. To access the values we have to redefine the _NL_CURRENT and | |
82 | _NL_CURRENT_WORD macros. */ | |
83 | #undef _NL_CURRENT | |
84 | #define _NL_CURRENT(category, item) \ | |
85 | (current->values[_NL_ITEM_INDEX (item)].string) | |
86 | #undef _NL_CURRENT_WORD | |
87 | #define _NL_CURRENT_WORD(category, item) \ | |
88 | ((uint32_t) current->values[_NL_ITEM_INDEX (item)].word) | |
89 | ||
90 | #if defined _LIBC || defined HAVE_WCHAR_H | |
91 | # include <wchar.h> | |
92 | #endif | |
93 | ||
94 | #ifdef USE_WIDE_CHAR | |
95 | # include <wctype.h> | |
96 | # define STRING_TYPE wchar_t | |
97 | # define CHAR_TYPE wint_t | |
98 | # define L_(Ch) L##Ch | |
99 | # define ISSPACE(Ch) __iswspace_l ((Ch), loc) | |
100 | # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc) | |
101 | # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc) | |
102 | # define TOLOWER(Ch) __towlower_l ((Ch), loc) | |
103 | # define STRNCASECMP(S1, S2, N) __wcsncasecmp_l ((S1), (S2), (N), loc) | |
104 | # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc) | |
105 | #else | |
106 | # define STRING_TYPE char | |
107 | # define CHAR_TYPE char | |
108 | # define L_(Ch) Ch | |
109 | # define ISSPACE(Ch) __isspace_l ((Ch), loc) | |
110 | # define ISDIGIT(Ch) __isdigit_l ((Ch), loc) | |
111 | # define ISXDIGIT(Ch) __isxdigit_l ((Ch), loc) | |
112 | # define TOLOWER(Ch) __tolower_l ((Ch), loc) | |
113 | # define STRNCASECMP(S1, S2, N) __strncasecmp_l ((S1), (S2), (N), loc) | |
114 | # define STRTOULL(S, E, B) ____strtoull_l_internal ((S), (E), (B), 0, loc) | |
115 | #endif | |
116 | ||
117 | ||
118 | /* Constants we need from float.h; select the set for the FLOAT precision. */ | |
119 | #define MANT_DIG PASTE(FLT,_MANT_DIG) | |
120 | #define DIG PASTE(FLT,_DIG) | |
121 | #define MAX_EXP PASTE(FLT,_MAX_EXP) | |
122 | #define MIN_EXP PASTE(FLT,_MIN_EXP) | |
123 | #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP) | |
124 | #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP) | |
125 | ||
126 | /* Extra macros required to get FLT expanded before the pasting. */ | |
127 | #define PASTE(a,b) PASTE1(a,b) | |
128 | #define PASTE1(a,b) a##b | |
129 | ||
130 | /* Function to construct a floating point number from an MP integer | |
131 | containing the fraction bits, a base 2 exponent, and a sign flag. */ | |
132 | extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative); | |
133 | \f | |
134 | /* Definitions according to limb size used. */ | |
135 | #if BITS_PER_MP_LIMB == 32 | |
136 | # define MAX_DIG_PER_LIMB 9 | |
137 | # define MAX_FAC_PER_LIMB 1000000000UL | |
138 | #elif BITS_PER_MP_LIMB == 64 | |
139 | # define MAX_DIG_PER_LIMB 19 | |
140 | # define MAX_FAC_PER_LIMB 10000000000000000000ULL | |
141 | #else | |
142 | # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for" | |
143 | #endif | |
144 | ||
145 | ||
146 | /* Local data structure. */ | |
147 | static const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] = | |
148 | { 0, 10, 100, | |
149 | 1000, 10000, 100000L, | |
150 | 1000000L, 10000000L, 100000000L, | |
151 | 1000000000L | |
152 | #if BITS_PER_MP_LIMB > 32 | |
153 | , 10000000000ULL, 100000000000ULL, | |
154 | 1000000000000ULL, 10000000000000ULL, 100000000000000ULL, | |
155 | 1000000000000000ULL, 10000000000000000ULL, 100000000000000000ULL, | |
156 | 1000000000000000000ULL, 10000000000000000000ULL | |
157 | #endif | |
158 | #if BITS_PER_MP_LIMB > 64 | |
159 | #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB | |
160 | #endif | |
161 | }; | |
162 | \f | |
163 | #ifndef howmany | |
164 | #define howmany(x,y) (((x)+((y)-1))/(y)) | |
165 | #endif | |
166 | #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; }) | |
167 | ||
168 | #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG) | |
169 | #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG) | |
170 | #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB) | |
171 | ||
172 | #define RETURN(val,end) \ | |
173 | do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \ | |
174 | return val; } while (0) | |
175 | ||
176 | /* Maximum size necessary for mpn integers to hold floating point numbers. */ | |
177 | #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \ | |
178 | + 2) | |
179 | /* Declare an mpn integer variable that big. */ | |
180 | #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size | |
181 | /* Copy an mpn integer value. */ | |
182 | #define MPN_ASSIGN(dst, src) \ | |
183 | memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t)) | |
184 | ||
185 | ||
186 | /* Return a floating point number of the needed type according to the given | |
187 | multi-precision number after possible rounding. */ | |
188 | static FLOAT | |
189 | round_and_return (mp_limb_t *retval, int exponent, int negative, | |
190 | mp_limb_t round_limb, mp_size_t round_bit, int more_bits) | |
191 | { | |
192 | if (exponent < MIN_EXP - 1) | |
193 | { | |
194 | mp_size_t shift = MIN_EXP - 1 - exponent; | |
195 | ||
196 | if (shift > MANT_DIG) | |
197 | { | |
198 | __set_errno (EDOM); | |
199 | return 0.0; | |
200 | } | |
201 | ||
202 | more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0; | |
203 | if (shift == MANT_DIG) | |
204 | /* This is a special case to handle the very seldom case where | |
205 | the mantissa will be empty after the shift. */ | |
206 | { | |
207 | int i; | |
208 | ||
209 | round_limb = retval[RETURN_LIMB_SIZE - 1]; | |
210 | round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; | |
211 | for (i = 0; i < RETURN_LIMB_SIZE; ++i) | |
212 | more_bits |= retval[i] != 0; | |
213 | MPN_ZERO (retval, RETURN_LIMB_SIZE); | |
214 | } | |
215 | else if (shift >= BITS_PER_MP_LIMB) | |
216 | { | |
217 | int i; | |
218 | ||
219 | round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB]; | |
220 | round_bit = (shift - 1) % BITS_PER_MP_LIMB; | |
221 | for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i) | |
222 | more_bits |= retval[i] != 0; | |
223 | more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) | |
224 | != 0); | |
225 | ||
226 | (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB], | |
227 | RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB), | |
228 | shift % BITS_PER_MP_LIMB); | |
229 | MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)], | |
230 | shift / BITS_PER_MP_LIMB); | |
231 | } | |
232 | else if (shift > 0) | |
233 | { | |
234 | round_limb = retval[0]; | |
235 | round_bit = shift - 1; | |
236 | (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift); | |
237 | } | |
238 | /* This is a hook for the m68k long double format, where the | |
239 | exponent bias is the same for normalized and denormalized | |
240 | numbers. */ | |
241 | #ifndef DENORM_EXP | |
242 | # define DENORM_EXP (MIN_EXP - 2) | |
243 | #endif | |
244 | exponent = DENORM_EXP; | |
245 | } | |
246 | ||
247 | if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0 | |
248 | && (more_bits || (retval[0] & 1) != 0 | |
249 | || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0)) | |
250 | { | |
251 | mp_limb_t cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1); | |
252 | ||
253 | if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || | |
254 | ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && | |
255 | (retval[RETURN_LIMB_SIZE - 1] | |
256 | & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)) | |
257 | { | |
258 | ++exponent; | |
259 | (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1); | |
260 | retval[RETURN_LIMB_SIZE - 1] | |
261 | |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB); | |
262 | } | |
263 | else if (exponent == DENORM_EXP | |
264 | && (retval[RETURN_LIMB_SIZE - 1] | |
265 | & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) | |
266 | != 0) | |
267 | /* The number was denormalized but now normalized. */ | |
268 | exponent = MIN_EXP - 1; | |
269 | } | |
270 | ||
271 | if (exponent > MAX_EXP) | |
272 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
273 | ||
274 | return MPN2FLOAT (retval, exponent, negative); | |
275 | } | |
276 | ||
277 | ||
278 | /* Read a multi-precision integer starting at STR with exactly DIGCNT digits | |
279 | into N. Return the size of the number limbs in NSIZE at the first | |
280 | character od the string that is not part of the integer as the function | |
281 | value. If the EXPONENT is small enough to be taken as an additional | |
282 | factor for the resulting number (see code) multiply by it. */ | |
283 | static const STRING_TYPE * | |
284 | str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize, | |
285 | int *exponent | |
286 | #ifndef USE_WIDE_CHAR | |
287 | , const char *decimal, size_t decimal_len, const char *thousands | |
288 | #endif | |
289 | ||
290 | ) | |
291 | { | |
292 | /* Number of digits for actual limb. */ | |
293 | int cnt = 0; | |
294 | mp_limb_t low = 0; | |
295 | mp_limb_t start; | |
296 | ||
297 | *nsize = 0; | |
298 | assert (digcnt > 0); | |
299 | do | |
300 | { | |
301 | if (cnt == MAX_DIG_PER_LIMB) | |
302 | { | |
303 | if (*nsize == 0) | |
304 | { | |
305 | n[0] = low; | |
306 | *nsize = 1; | |
307 | } | |
308 | else | |
309 | { | |
310 | mp_limb_t cy; | |
311 | cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB); | |
312 | cy += __mpn_add_1 (n, n, *nsize, low); | |
313 | if (cy != 0) | |
314 | { | |
315 | n[*nsize] = cy; | |
316 | ++(*nsize); | |
317 | } | |
318 | } | |
319 | cnt = 0; | |
320 | low = 0; | |
321 | } | |
322 | ||
323 | /* There might be thousands separators or radix characters in | |
324 | the string. But these all can be ignored because we know the | |
325 | format of the number is correct and we have an exact number | |
326 | of characters to read. */ | |
327 | #ifdef USE_WIDE_CHAR | |
328 | if (*str < L'0' || *str > L'9') | |
329 | ++str; | |
330 | #else | |
331 | if (*str < '0' || *str > '9') | |
332 | { | |
333 | int inner = 0; | |
334 | if (thousands != NULL && *str == *thousands | |
335 | && ({ for (inner = 1; thousands[inner] != '\0'; ++inner) | |
336 | if (thousands[inner] != str[inner]) | |
337 | break; | |
338 | thousands[inner] == '\0'; })) | |
339 | str += inner; | |
340 | else | |
341 | str += decimal_len; | |
342 | } | |
343 | #endif | |
344 | low = low * 10 + *str++ - L_('0'); | |
345 | ++cnt; | |
346 | } | |
347 | while (--digcnt > 0); | |
348 | ||
349 | if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB) | |
350 | { | |
351 | low *= _tens_in_limb[*exponent]; | |
352 | start = _tens_in_limb[cnt + *exponent]; | |
353 | *exponent = 0; | |
354 | } | |
355 | else | |
356 | start = _tens_in_limb[cnt]; | |
357 | ||
358 | if (*nsize == 0) | |
359 | { | |
360 | n[0] = low; | |
361 | *nsize = 1; | |
362 | } | |
363 | else | |
364 | { | |
365 | mp_limb_t cy; | |
366 | cy = __mpn_mul_1 (n, n, *nsize, start); | |
367 | cy += __mpn_add_1 (n, n, *nsize, low); | |
368 | if (cy != 0) | |
369 | n[(*nsize)++] = cy; | |
370 | } | |
371 | ||
372 | return str; | |
373 | } | |
374 | ||
375 | ||
376 | /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits | |
377 | with the COUNT most significant bits of LIMB. | |
378 | ||
379 | Tege doesn't like this function so I have to write it here myself. :) | |
380 | --drepper */ | |
381 | static inline void | |
382 | __attribute ((always_inline)) | |
383 | __mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count, | |
384 | mp_limb_t limb) | |
385 | { | |
386 | if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) | |
387 | { | |
388 | /* Optimize the case of shifting by exactly a word: | |
389 | just copy words, with no actual bit-shifting. */ | |
390 | mp_size_t i; | |
391 | for (i = size - 1; i > 0; --i) | |
392 | ptr[i] = ptr[i - 1]; | |
393 | ptr[0] = limb; | |
394 | } | |
395 | else | |
396 | { | |
397 | (void) __mpn_lshift (ptr, ptr, size, count); | |
398 | ptr[0] |= limb >> (BITS_PER_MP_LIMB - count); | |
399 | } | |
400 | } | |
401 | ||
402 | ||
403 | #define INTERNAL(x) INTERNAL1(x) | |
404 | #define INTERNAL1(x) __##x##_internal | |
405 | ||
406 | /* This file defines a function to check for correct grouping. */ | |
407 | #include "grouping.h" | |
408 | ||
409 | ||
410 | /* Return a floating point number with the value of the given string NPTR. | |
411 | Set *ENDPTR to the character after the last used one. If the number is | |
412 | smaller than the smallest representable number, set `errno' to ERANGE and | |
413 | return 0.0. If the number is too big to be represented, set `errno' to | |
414 | ERANGE and return HUGE_VAL with the appropriate sign. */ | |
415 | FLOAT | |
416 | INTERNAL (__STRTOF) (nptr, endptr, group, loc) | |
417 | const STRING_TYPE *nptr; | |
418 | STRING_TYPE **endptr; | |
419 | int group; | |
420 | __locale_t loc; | |
421 | { | |
422 | int negative; /* The sign of the number. */ | |
423 | MPN_VAR (num); /* MP representation of the number. */ | |
424 | int exponent; /* Exponent of the number. */ | |
425 | ||
426 | /* Numbers starting `0X' or `0x' have to be processed with base 16. */ | |
427 | int base = 10; | |
428 | ||
429 | /* When we have to compute fractional digits we form a fraction with a | |
430 | second multi-precision number (and we sometimes need a second for | |
431 | temporary results). */ | |
432 | MPN_VAR (den); | |
433 | ||
434 | /* Representation for the return value. */ | |
435 | mp_limb_t retval[RETURN_LIMB_SIZE]; | |
436 | /* Number of bits currently in result value. */ | |
437 | int bits; | |
438 | ||
439 | /* Running pointer after the last character processed in the string. */ | |
440 | const STRING_TYPE *cp, *tp; | |
441 | /* Start of significant part of the number. */ | |
442 | const STRING_TYPE *startp, *start_of_digits; | |
443 | /* Points at the character following the integer and fractional digits. */ | |
444 | const STRING_TYPE *expp; | |
445 | /* Total number of digit and number of digits in integer part. */ | |
446 | int dig_no, int_no, lead_zero; | |
447 | /* Contains the last character read. */ | |
448 | CHAR_TYPE c; | |
449 | ||
450 | /* We should get wint_t from <stddef.h>, but not all GCC versions define it | |
451 | there. So define it ourselves if it remains undefined. */ | |
452 | #ifndef _WINT_T | |
453 | typedef unsigned int wint_t; | |
454 | #endif | |
455 | /* The radix character of the current locale. */ | |
456 | #ifdef USE_WIDE_CHAR | |
457 | wchar_t decimal; | |
458 | #else | |
459 | const char *decimal; | |
460 | size_t decimal_len; | |
461 | #endif | |
462 | /* The thousands character of the current locale. */ | |
463 | #ifdef USE_WIDE_CHAR | |
464 | wchar_t thousands = L'\0'; | |
465 | #else | |
466 | const char *thousands = NULL; | |
467 | #endif | |
468 | /* The numeric grouping specification of the current locale, | |
469 | in the format described in <locale.h>. */ | |
470 | const char *grouping; | |
471 | /* Used in several places. */ | |
472 | int cnt; | |
473 | ||
474 | struct locale_data *current = loc->__locales[LC_NUMERIC]; | |
475 | ||
476 | if (group) | |
477 | { | |
478 | grouping = _NL_CURRENT (LC_NUMERIC, GROUPING); | |
479 | if (*grouping <= 0 || *grouping == CHAR_MAX) | |
480 | grouping = NULL; | |
481 | else | |
482 | { | |
483 | /* Figure out the thousands separator character. */ | |
484 | #ifdef USE_WIDE_CHAR | |
485 | thousands = _NL_CURRENT_WORD (LC_NUMERIC, | |
486 | _NL_NUMERIC_THOUSANDS_SEP_WC); | |
487 | if (thousands == L'\0') | |
488 | grouping = NULL; | |
489 | #else | |
490 | thousands = _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP); | |
491 | if (*thousands == '\0') | |
492 | { | |
493 | thousands = NULL; | |
494 | grouping = NULL; | |
495 | } | |
496 | #endif | |
497 | } | |
498 | } | |
499 | else | |
500 | grouping = NULL; | |
501 | ||
502 | /* Find the locale's decimal point character. */ | |
503 | #ifdef USE_WIDE_CHAR | |
504 | decimal = _NL_CURRENT_WORD (LC_NUMERIC, _NL_NUMERIC_DECIMAL_POINT_WC); | |
505 | assert (decimal != L'\0'); | |
506 | # define decimal_len 1 | |
507 | #else | |
508 | decimal = _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT); | |
509 | decimal_len = strlen (decimal); | |
510 | assert (decimal_len > 0); | |
511 | #endif | |
512 | ||
513 | /* Prepare number representation. */ | |
514 | exponent = 0; | |
515 | negative = 0; | |
516 | bits = 0; | |
517 | ||
518 | /* Parse string to get maximal legal prefix. We need the number of | |
519 | characters of the integer part, the fractional part and the exponent. */ | |
520 | cp = nptr - 1; | |
521 | /* Ignore leading white space. */ | |
522 | do | |
523 | c = *++cp; | |
524 | while (ISSPACE (c)); | |
525 | ||
526 | /* Get sign of the result. */ | |
527 | if (c == L_('-')) | |
528 | { | |
529 | negative = 1; | |
530 | c = *++cp; | |
531 | } | |
532 | else if (c == L_('+')) | |
533 | c = *++cp; | |
534 | ||
535 | /* Return 0.0 if no legal string is found. | |
536 | No character is used even if a sign was found. */ | |
537 | #ifdef USE_WIDE_CHAR | |
538 | if (c == (wint_t) decimal | |
539 | && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9') | |
540 | { | |
541 | /* We accept it. This funny construct is here only to indent | |
542 | the code directly. */ | |
543 | } | |
544 | #else | |
545 | for (cnt = 0; decimal[cnt] != '\0'; ++cnt) | |
546 | if (cp[cnt] != decimal[cnt]) | |
547 | break; | |
548 | if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9') | |
549 | { | |
550 | /* We accept it. This funny construct is here only to indent | |
551 | the code directly. */ | |
552 | } | |
553 | #endif | |
554 | else if (c < L_('0') || c > L_('9')) | |
555 | { | |
556 | /* Check for `INF' or `INFINITY'. */ | |
557 | if (TOLOWER (c) == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0) | |
558 | { | |
559 | /* Return +/- infinity. */ | |
560 | if (endptr != NULL) | |
561 | *endptr = (STRING_TYPE *) | |
562 | (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0 | |
563 | ? 8 : 3)); | |
564 | ||
565 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
566 | } | |
567 | ||
568 | if (TOLOWER (c) == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0) | |
569 | { | |
570 | /* Return NaN. */ | |
571 | FLOAT retval = NAN; | |
572 | ||
573 | cp += 3; | |
574 | ||
575 | /* Match `(n-char-sequence-digit)'. */ | |
576 | if (*cp == L_('(')) | |
577 | { | |
578 | const STRING_TYPE *startp = cp; | |
579 | do | |
580 | ++cp; | |
581 | while ((*cp >= L_('0') && *cp <= L_('9')) | |
582 | || (TOLOWER (*cp) >= L_('a') && TOLOWER (*cp) <= L_('z')) | |
583 | || *cp == L_('_')); | |
584 | ||
585 | if (*cp != L_(')')) | |
586 | /* The closing brace is missing. Only match the NAN | |
587 | part. */ | |
588 | cp = startp; | |
589 | else | |
590 | { | |
591 | /* This is a system-dependent way to specify the | |
592 | bitmask used for the NaN. We expect it to be | |
593 | a number which is put in the mantissa of the | |
594 | number. */ | |
595 | STRING_TYPE *endp; | |
596 | unsigned long long int mant; | |
597 | ||
598 | mant = STRTOULL (startp + 1, &endp, 0); | |
599 | if (endp == cp) | |
600 | SET_MANTISSA (retval, mant); | |
601 | } | |
602 | } | |
603 | ||
604 | if (endptr != NULL) | |
605 | *endptr = (STRING_TYPE *) cp; | |
606 | ||
607 | return retval; | |
608 | } | |
609 | ||
610 | /* It is really a text we do not recognize. */ | |
611 | RETURN (0.0, nptr); | |
612 | } | |
613 | ||
614 | /* First look whether we are faced with a hexadecimal number. */ | |
615 | if (c == L_('0') && TOLOWER (cp[1]) == L_('x')) | |
616 | { | |
617 | /* Okay, it is a hexa-decimal number. Remember this and skip | |
618 | the characters. BTW: hexadecimal numbers must not be | |
619 | grouped. */ | |
620 | base = 16; | |
621 | cp += 2; | |
622 | c = *cp; | |
623 | grouping = NULL; | |
624 | } | |
625 | ||
626 | /* Record the start of the digits, in case we will check their grouping. */ | |
627 | start_of_digits = startp = cp; | |
628 | ||
629 | /* Ignore leading zeroes. This helps us to avoid useless computations. */ | |
630 | #ifdef USE_WIDE_CHAR | |
631 | while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands)) | |
632 | c = *++cp; | |
633 | #else | |
634 | if (thousands == NULL) | |
635 | while (c == '0') | |
636 | c = *++cp; | |
637 | else | |
638 | { | |
639 | /* We also have the multibyte thousands string. */ | |
640 | while (1) | |
641 | { | |
642 | if (c != '0') | |
643 | { | |
644 | for (cnt = 0; thousands[cnt] != '\0'; ++cnt) | |
645 | if (c != thousands[cnt]) | |
646 | break; | |
647 | if (thousands[cnt] != '\0') | |
648 | break; | |
649 | } | |
650 | c = *++cp; | |
651 | } | |
652 | } | |
653 | #endif | |
654 | ||
655 | /* If no other digit but a '0' is found the result is 0.0. | |
656 | Return current read pointer. */ | |
657 | if ((c < L_('0') || c > L_('9')) | |
658 | && (base == 16 && (c < (CHAR_TYPE) TOLOWER (L_('a')) | |
659 | || c > (CHAR_TYPE) TOLOWER (L_('f')))) | |
660 | #ifdef USE_WIDE_CHAR | |
661 | && c != (wint_t) decimal | |
662 | #else | |
663 | && ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) | |
664 | if (decimal[cnt] != cp[cnt]) | |
665 | break; | |
666 | decimal[cnt] != '\0'; }) | |
667 | #endif | |
668 | && (base == 16 && (cp == start_of_digits | |
669 | || (CHAR_TYPE) TOLOWER (c) != L_('p'))) | |
670 | && (base != 16 && (CHAR_TYPE) TOLOWER (c) != L_('e'))) | |
671 | { | |
672 | #ifdef USE_WIDE_CHAR | |
673 | tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, | |
674 | grouping); | |
675 | #else | |
676 | tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, | |
677 | grouping); | |
678 | #endif | |
679 | /* If TP is at the start of the digits, there was no correctly | |
680 | grouped prefix of the string; so no number found. */ | |
681 | RETURN (0.0, tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp); | |
682 | } | |
683 | ||
684 | /* Remember first significant digit and read following characters until the | |
685 | decimal point, exponent character or any non-FP number character. */ | |
686 | startp = cp; | |
687 | dig_no = 0; | |
688 | while (1) | |
689 | { | |
690 | if ((c >= L_('0') && c <= L_('9')) | |
691 | || (base == 16 && (wint_t) TOLOWER (c) >= L_('a') | |
692 | && (wint_t) TOLOWER (c) <= L_('f'))) | |
693 | ++dig_no; | |
694 | else | |
695 | { | |
696 | #ifdef USE_WIDE_CHAR | |
697 | if ((wint_t) thousands == L'\0' || c != (wint_t) thousands) | |
698 | /* Not a digit or separator: end of the integer part. */ | |
699 | break; | |
700 | #else | |
701 | if (thousands == NULL) | |
702 | break; | |
703 | else | |
704 | { | |
705 | for (cnt = 0; thousands[cnt] != '\0'; ++cnt) | |
706 | if (thousands[cnt] != cp[cnt]) | |
707 | break; | |
708 | if (thousands[cnt] != '\0') | |
709 | break; | |
710 | } | |
711 | #endif | |
712 | } | |
713 | c = *++cp; | |
714 | } | |
715 | ||
716 | if (grouping && dig_no > 0) | |
717 | { | |
718 | /* Check the grouping of the digits. */ | |
719 | #ifdef USE_WIDE_CHAR | |
720 | tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, | |
721 | grouping); | |
722 | #else | |
723 | tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, | |
724 | grouping); | |
725 | #endif | |
726 | if (cp != tp) | |
727 | { | |
728 | /* Less than the entire string was correctly grouped. */ | |
729 | ||
730 | if (tp == start_of_digits) | |
731 | /* No valid group of numbers at all: no valid number. */ | |
732 | RETURN (0.0, nptr); | |
733 | ||
734 | if (tp < startp) | |
735 | /* The number is validly grouped, but consists | |
736 | only of zeroes. The whole value is zero. */ | |
737 | RETURN (0.0, tp); | |
738 | ||
739 | /* Recompute DIG_NO so we won't read more digits than | |
740 | are properly grouped. */ | |
741 | cp = tp; | |
742 | dig_no = 0; | |
743 | for (tp = startp; tp < cp; ++tp) | |
744 | if (*tp >= L_('0') && *tp <= L_('9')) | |
745 | ++dig_no; | |
746 | ||
747 | int_no = dig_no; | |
748 | lead_zero = 0; | |
749 | ||
750 | goto number_parsed; | |
751 | } | |
752 | } | |
753 | ||
754 | /* We have the number digits in the integer part. Whether these are all or | |
755 | any is really a fractional digit will be decided later. */ | |
756 | int_no = dig_no; | |
757 | lead_zero = int_no == 0 ? -1 : 0; | |
758 | ||
759 | /* Read the fractional digits. A special case are the 'american style' | |
760 | numbers like `16.' i.e. with decimal but without trailing digits. */ | |
761 | if ( | |
762 | #ifdef USE_WIDE_CHAR | |
763 | c == (wint_t) decimal | |
764 | #else | |
765 | ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) | |
766 | if (decimal[cnt] != cp[cnt]) | |
767 | break; | |
768 | decimal[cnt] == '\0'; }) | |
769 | #endif | |
770 | ) | |
771 | { | |
772 | cp += decimal_len; | |
773 | c = *cp; | |
774 | while ((c >= L_('0') && c <= L_('9')) || | |
775 | (base == 16 && TOLOWER (c) >= L_('a') && TOLOWER (c) <= L_('f'))) | |
776 | { | |
777 | if (c != L_('0') && lead_zero == -1) | |
778 | lead_zero = dig_no - int_no; | |
779 | ++dig_no; | |
780 | c = *++cp; | |
781 | } | |
782 | } | |
783 | ||
784 | /* Remember start of exponent (if any). */ | |
785 | expp = cp; | |
786 | ||
787 | /* Read exponent. */ | |
788 | if ((base == 16 && TOLOWER (c) == L_('p')) | |
789 | || (base != 16 && TOLOWER (c) == L_('e'))) | |
790 | { | |
791 | int exp_negative = 0; | |
792 | ||
793 | c = *++cp; | |
794 | if (c == L_('-')) | |
795 | { | |
796 | exp_negative = 1; | |
797 | c = *++cp; | |
798 | } | |
799 | else if (c == L_('+')) | |
800 | c = *++cp; | |
801 | ||
802 | if (c >= L_('0') && c <= L_('9')) | |
803 | { | |
804 | int exp_limit; | |
805 | ||
806 | /* Get the exponent limit. */ | |
807 | if (base == 16) | |
808 | exp_limit = (exp_negative ? | |
809 | -MIN_EXP + MANT_DIG + 4 * int_no : | |
810 | MAX_EXP - 4 * int_no + lead_zero); | |
811 | else | |
812 | exp_limit = (exp_negative ? | |
813 | -MIN_10_EXP + MANT_DIG + int_no : | |
814 | MAX_10_EXP - int_no + lead_zero); | |
815 | ||
816 | do | |
817 | { | |
818 | exponent *= 10; | |
819 | ||
820 | if (exponent > exp_limit) | |
821 | /* The exponent is too large/small to represent a valid | |
822 | number. */ | |
823 | { | |
824 | FLOAT result; | |
825 | ||
826 | /* We have to take care for special situation: a joker | |
827 | might have written "0.0e100000" which is in fact | |
828 | zero. */ | |
829 | if (lead_zero == -1) | |
830 | result = negative ? -0.0 : 0.0; | |
831 | else | |
832 | { | |
833 | /* Overflow or underflow. */ | |
834 | __set_errno (ERANGE); | |
835 | result = (exp_negative ? 0.0 : | |
836 | negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL); | |
837 | } | |
838 | ||
839 | /* Accept all following digits as part of the exponent. */ | |
840 | do | |
841 | ++cp; | |
842 | while (*cp >= L_('0') && *cp <= L_('9')); | |
843 | ||
844 | RETURN (result, cp); | |
845 | /* NOTREACHED */ | |
846 | } | |
847 | ||
848 | exponent += c - L_('0'); | |
849 | c = *++cp; | |
850 | } | |
851 | while (c >= L_('0') && c <= L_('9')); | |
852 | ||
853 | if (exp_negative) | |
854 | exponent = -exponent; | |
855 | } | |
856 | else | |
857 | cp = expp; | |
858 | } | |
859 | ||
860 | /* We don't want to have to work with trailing zeroes after the radix. */ | |
861 | if (dig_no > int_no) | |
862 | { | |
863 | while (expp[-1] == L_('0')) | |
864 | { | |
865 | --expp; | |
866 | --dig_no; | |
867 | } | |
868 | assert (dig_no >= int_no); | |
869 | } | |
870 | ||
871 | if (dig_no == int_no && dig_no > 0 && exponent < 0) | |
872 | do | |
873 | { | |
874 | while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1]))) | |
875 | --expp; | |
876 | ||
877 | if (expp[-1] != L_('0')) | |
878 | break; | |
879 | ||
880 | --expp; | |
881 | --dig_no; | |
882 | --int_no; | |
883 | ++exponent; | |
884 | } | |
885 | while (dig_no > 0 && exponent < 0); | |
886 | ||
887 | number_parsed: | |
888 | ||
889 | /* The whole string is parsed. Store the address of the next character. */ | |
890 | if (endptr) | |
891 | *endptr = (STRING_TYPE *) cp; | |
892 | ||
893 | if (dig_no == 0) | |
894 | return negative ? -0.0 : 0.0; | |
895 | ||
896 | if (lead_zero) | |
897 | { | |
898 | /* Find the decimal point */ | |
899 | #ifdef USE_WIDE_CHAR | |
900 | while (*startp != decimal) | |
901 | ++startp; | |
902 | #else | |
903 | while (1) | |
904 | { | |
905 | if (*startp == decimal[0]) | |
906 | { | |
907 | for (cnt = 1; decimal[cnt] != '\0'; ++cnt) | |
908 | if (decimal[cnt] != startp[cnt]) | |
909 | break; | |
910 | if (decimal[cnt] == '\0') | |
911 | break; | |
912 | } | |
913 | ++startp; | |
914 | } | |
915 | #endif | |
916 | startp += lead_zero + decimal_len; | |
917 | exponent -= base == 16 ? 4 * lead_zero : lead_zero; | |
918 | dig_no -= lead_zero; | |
919 | } | |
920 | ||
921 | /* If the BASE is 16 we can use a simpler algorithm. */ | |
922 | if (base == 16) | |
923 | { | |
924 | static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3, | |
925 | 4, 4, 4, 4, 4, 4, 4, 4 }; | |
926 | int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB; | |
927 | int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB; | |
928 | mp_limb_t val; | |
929 | ||
930 | while (!ISXDIGIT (*startp)) | |
931 | ++startp; | |
932 | while (*startp == L_('0')) | |
933 | ++startp; | |
934 | if (ISDIGIT (*startp)) | |
935 | val = *startp++ - L_('0'); | |
936 | else | |
937 | val = 10 + TOLOWER (*startp++) - L_('a'); | |
938 | bits = nbits[val]; | |
939 | /* We cannot have a leading zero. */ | |
940 | assert (bits != 0); | |
941 | ||
942 | if (pos + 1 >= 4 || pos + 1 >= bits) | |
943 | { | |
944 | /* We don't have to care for wrapping. This is the normal | |
945 | case so we add the first clause in the `if' expression as | |
946 | an optimization. It is a compile-time constant and so does | |
947 | not cost anything. */ | |
948 | retval[idx] = val << (pos - bits + 1); | |
949 | pos -= bits; | |
950 | } | |
951 | else | |
952 | { | |
953 | retval[idx--] = val >> (bits - pos - 1); | |
954 | retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1)); | |
955 | pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1); | |
956 | } | |
957 | ||
958 | /* Adjust the exponent for the bits we are shifting in. */ | |
959 | exponent += bits - 1 + (int_no - 1) * 4; | |
960 | ||
961 | while (--dig_no > 0 && idx >= 0) | |
962 | { | |
963 | if (!ISXDIGIT (*startp)) | |
964 | startp += decimal_len; | |
965 | if (ISDIGIT (*startp)) | |
966 | val = *startp++ - L_('0'); | |
967 | else | |
968 | val = 10 + TOLOWER (*startp++) - L_('a'); | |
969 | ||
970 | if (pos + 1 >= 4) | |
971 | { | |
972 | retval[idx] |= val << (pos - 4 + 1); | |
973 | pos -= 4; | |
974 | } | |
975 | else | |
976 | { | |
977 | retval[idx--] |= val >> (4 - pos - 1); | |
978 | val <<= BITS_PER_MP_LIMB - (4 - pos - 1); | |
979 | if (idx < 0) | |
980 | return round_and_return (retval, exponent, negative, val, | |
981 | BITS_PER_MP_LIMB - 1, dig_no > 0); | |
982 | ||
983 | retval[idx] = val; | |
984 | pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1); | |
985 | } | |
986 | } | |
987 | ||
988 | /* We ran out of digits. */ | |
989 | MPN_ZERO (retval, idx); | |
990 | ||
991 | return round_and_return (retval, exponent, negative, 0, 0, 0); | |
992 | } | |
993 | ||
994 | /* Now we have the number of digits in total and the integer digits as well | |
995 | as the exponent and its sign. We can decide whether the read digits are | |
996 | really integer digits or belong to the fractional part; i.e. we normalize | |
997 | 123e-2 to 1.23. */ | |
998 | { | |
999 | register int incr = (exponent < 0 ? MAX (-int_no, exponent) | |
1000 | : MIN (dig_no - int_no, exponent)); | |
1001 | int_no += incr; | |
1002 | exponent -= incr; | |
1003 | } | |
1004 | ||
1005 | if (int_no + exponent > MAX_10_EXP + 1) | |
1006 | { | |
1007 | __set_errno (ERANGE); | |
1008 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
1009 | } | |
1010 | ||
1011 | if (exponent < MIN_10_EXP - (DIG + 1)) | |
1012 | { | |
1013 | __set_errno (ERANGE); | |
1014 | return 0.0; | |
1015 | } | |
1016 | ||
1017 | if (int_no > 0) | |
1018 | { | |
1019 | /* Read the integer part as a multi-precision number to NUM. */ | |
1020 | startp = str_to_mpn (startp, int_no, num, &numsize, &exponent | |
1021 | #ifndef USE_WIDE_CHAR | |
1022 | , decimal, decimal_len, thousands | |
1023 | #endif | |
1024 | ); | |
1025 | ||
1026 | if (exponent > 0) | |
1027 | { | |
1028 | /* We now multiply the gained number by the given power of ten. */ | |
1029 | mp_limb_t *psrc = num; | |
1030 | mp_limb_t *pdest = den; | |
1031 | int expbit = 1; | |
1032 | const struct mp_power *ttab = &_fpioconst_pow10[0]; | |
1033 | ||
1034 | do | |
1035 | { | |
1036 | if ((exponent & expbit) != 0) | |
1037 | { | |
1038 | size_t size = ttab->arraysize - _FPIO_CONST_OFFSET; | |
1039 | mp_limb_t cy; | |
1040 | exponent ^= expbit; | |
1041 | ||
1042 | /* FIXME: not the whole multiplication has to be | |
1043 | done. If we have the needed number of bits we | |
1044 | only need the information whether more non-zero | |
1045 | bits follow. */ | |
1046 | if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET) | |
1047 | cy = __mpn_mul (pdest, psrc, numsize, | |
1048 | &__tens[ttab->arrayoff | |
1049 | + _FPIO_CONST_OFFSET], | |
1050 | size); | |
1051 | else | |
1052 | cy = __mpn_mul (pdest, &__tens[ttab->arrayoff | |
1053 | + _FPIO_CONST_OFFSET], | |
1054 | size, psrc, numsize); | |
1055 | numsize += size; | |
1056 | if (cy == 0) | |
1057 | --numsize; | |
1058 | (void) SWAP (psrc, pdest); | |
1059 | } | |
1060 | expbit <<= 1; | |
1061 | ++ttab; | |
1062 | } | |
1063 | while (exponent != 0); | |
1064 | ||
1065 | if (psrc == den) | |
1066 | memcpy (num, den, numsize * sizeof (mp_limb_t)); | |
1067 | } | |
1068 | ||
1069 | /* Determine how many bits of the result we already have. */ | |
1070 | count_leading_zeros (bits, num[numsize - 1]); | |
1071 | bits = numsize * BITS_PER_MP_LIMB - bits; | |
1072 | ||
1073 | /* Now we know the exponent of the number in base two. | |
1074 | Check it against the maximum possible exponent. */ | |
1075 | if (bits > MAX_EXP) | |
1076 | { | |
1077 | __set_errno (ERANGE); | |
1078 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
1079 | } | |
1080 | ||
1081 | /* We have already the first BITS bits of the result. Together with | |
1082 | the information whether more non-zero bits follow this is enough | |
1083 | to determine the result. */ | |
1084 | if (bits > MANT_DIG) | |
1085 | { | |
1086 | int i; | |
1087 | const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB; | |
1088 | const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB; | |
1089 | const mp_size_t round_idx = least_bit == 0 ? least_idx - 1 | |
1090 | : least_idx; | |
1091 | const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1 | |
1092 | : least_bit - 1; | |
1093 | ||
1094 | if (least_bit == 0) | |
1095 | memcpy (retval, &num[least_idx], | |
1096 | RETURN_LIMB_SIZE * sizeof (mp_limb_t)); | |
1097 | else | |
1098 | { | |
1099 | for (i = least_idx; i < numsize - 1; ++i) | |
1100 | retval[i - least_idx] = (num[i] >> least_bit) | |
1101 | | (num[i + 1] | |
1102 | << (BITS_PER_MP_LIMB - least_bit)); | |
1103 | if (i - least_idx < RETURN_LIMB_SIZE) | |
1104 | retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit; | |
1105 | } | |
1106 | ||
1107 | /* Check whether any limb beside the ones in RETVAL are non-zero. */ | |
1108 | for (i = 0; num[i] == 0; ++i) | |
1109 | ; | |
1110 | ||
1111 | return round_and_return (retval, bits - 1, negative, | |
1112 | num[round_idx], round_bit, | |
1113 | int_no < dig_no || i < round_idx); | |
1114 | /* NOTREACHED */ | |
1115 | } | |
1116 | else if (dig_no == int_no) | |
1117 | { | |
1118 | const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; | |
1119 | const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB; | |
1120 | ||
1121 | if (target_bit == is_bit) | |
1122 | { | |
1123 | memcpy (&retval[RETURN_LIMB_SIZE - numsize], num, | |
1124 | numsize * sizeof (mp_limb_t)); | |
1125 | /* FIXME: the following loop can be avoided if we assume a | |
1126 | maximal MANT_DIG value. */ | |
1127 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); | |
1128 | } | |
1129 | else if (target_bit > is_bit) | |
1130 | { | |
1131 | (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize], | |
1132 | num, numsize, target_bit - is_bit); | |
1133 | /* FIXME: the following loop can be avoided if we assume a | |
1134 | maximal MANT_DIG value. */ | |
1135 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); | |
1136 | } | |
1137 | else | |
1138 | { | |
1139 | mp_limb_t cy; | |
1140 | assert (numsize < RETURN_LIMB_SIZE); | |
1141 | ||
1142 | cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize], | |
1143 | num, numsize, is_bit - target_bit); | |
1144 | retval[RETURN_LIMB_SIZE - numsize - 1] = cy; | |
1145 | /* FIXME: the following loop can be avoided if we assume a | |
1146 | maximal MANT_DIG value. */ | |
1147 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1); | |
1148 | } | |
1149 | ||
1150 | return round_and_return (retval, bits - 1, negative, 0, 0, 0); | |
1151 | /* NOTREACHED */ | |
1152 | } | |
1153 | ||
1154 | /* Store the bits we already have. */ | |
1155 | memcpy (retval, num, numsize * sizeof (mp_limb_t)); | |
1156 | #if RETURN_LIMB_SIZE > 1 | |
1157 | if (numsize < RETURN_LIMB_SIZE) | |
1158 | retval[numsize] = 0; | |
1159 | #endif | |
1160 | } | |
1161 | ||
1162 | /* We have to compute at least some of the fractional digits. */ | |
1163 | { | |
1164 | /* We construct a fraction and the result of the division gives us | |
1165 | the needed digits. The denominator is 1.0 multiplied by the | |
1166 | exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and | |
1167 | 123e-6 gives 123 / 1000000. */ | |
1168 | ||
1169 | int expbit; | |
1170 | int neg_exp; | |
1171 | int more_bits; | |
1172 | mp_limb_t cy; | |
1173 | mp_limb_t *psrc = den; | |
1174 | mp_limb_t *pdest = num; | |
1175 | const struct mp_power *ttab = &_fpioconst_pow10[0]; | |
1176 | ||
1177 | assert (dig_no > int_no && exponent <= 0); | |
1178 | ||
1179 | ||
1180 | /* For the fractional part we need not process too many digits. One | |
1181 | decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute | |
1182 | ceil(BITS / 3) =: N | |
1183 | digits we should have enough bits for the result. The remaining | |
1184 | decimal digits give us the information that more bits are following. | |
1185 | This can be used while rounding. (Two added as a safety margin.) */ | |
1186 | if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 2) | |
1187 | { | |
1188 | dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 2; | |
1189 | more_bits = 1; | |
1190 | } | |
1191 | else | |
1192 | more_bits = 0; | |
1193 | ||
1194 | neg_exp = dig_no - int_no - exponent; | |
1195 | ||
1196 | /* Construct the denominator. */ | |
1197 | densize = 0; | |
1198 | expbit = 1; | |
1199 | do | |
1200 | { | |
1201 | if ((neg_exp & expbit) != 0) | |
1202 | { | |
1203 | mp_limb_t cy; | |
1204 | neg_exp ^= expbit; | |
1205 | ||
1206 | if (densize == 0) | |
1207 | { | |
1208 | densize = ttab->arraysize - _FPIO_CONST_OFFSET; | |
1209 | memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET], | |
1210 | densize * sizeof (mp_limb_t)); | |
1211 | } | |
1212 | else | |
1213 | { | |
1214 | cy = __mpn_mul (pdest, &__tens[ttab->arrayoff | |
1215 | + _FPIO_CONST_OFFSET], | |
1216 | ttab->arraysize - _FPIO_CONST_OFFSET, | |
1217 | psrc, densize); | |
1218 | densize += ttab->arraysize - _FPIO_CONST_OFFSET; | |
1219 | if (cy == 0) | |
1220 | --densize; | |
1221 | (void) SWAP (psrc, pdest); | |
1222 | } | |
1223 | } | |
1224 | expbit <<= 1; | |
1225 | ++ttab; | |
1226 | } | |
1227 | while (neg_exp != 0); | |
1228 | ||
1229 | if (psrc == num) | |
1230 | memcpy (den, num, densize * sizeof (mp_limb_t)); | |
1231 | ||
1232 | /* Read the fractional digits from the string. */ | |
1233 | (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent | |
1234 | #ifndef USE_WIDE_CHAR | |
1235 | , decimal, decimal_len, thousands | |
1236 | #endif | |
1237 | ); | |
1238 | ||
1239 | /* We now have to shift both numbers so that the highest bit in the | |
1240 | denominator is set. In the same process we copy the numerator to | |
1241 | a high place in the array so that the division constructs the wanted | |
1242 | digits. This is done by a "quasi fix point" number representation. | |
1243 | ||
1244 | num: ddddddddddd . 0000000000000000000000 | |
1245 | |--- m ---| | |
1246 | den: ddddddddddd n >= m | |
1247 | |--- n ---| | |
1248 | */ | |
1249 | ||
1250 | count_leading_zeros (cnt, den[densize - 1]); | |
1251 | ||
1252 | if (cnt > 0) | |
1253 | { | |
1254 | /* Don't call `mpn_shift' with a count of zero since the specification | |
1255 | does not allow this. */ | |
1256 | (void) __mpn_lshift (den, den, densize, cnt); | |
1257 | cy = __mpn_lshift (num, num, numsize, cnt); | |
1258 | if (cy != 0) | |
1259 | num[numsize++] = cy; | |
1260 | } | |
1261 | ||
1262 | /* Now we are ready for the division. But it is not necessary to | |
1263 | do a full multi-precision division because we only need a small | |
1264 | number of bits for the result. So we do not use __mpn_divmod | |
1265 | here but instead do the division here by hand and stop whenever | |
1266 | the needed number of bits is reached. The code itself comes | |
1267 | from the GNU MP Library by Torbj\"orn Granlund. */ | |
1268 | ||
1269 | exponent = bits; | |
1270 | ||
1271 | switch (densize) | |
1272 | { | |
1273 | case 1: | |
1274 | { | |
1275 | mp_limb_t d, n, quot; | |
1276 | int used = 0; | |
1277 | ||
1278 | n = num[0]; | |
1279 | d = den[0]; | |
1280 | assert (numsize == 1 && n < d); | |
1281 | ||
1282 | do | |
1283 | { | |
1284 | udiv_qrnnd (quot, n, n, 0, d); | |
1285 | ||
1286 | #define got_limb \ | |
1287 | if (bits == 0) \ | |
1288 | { \ | |
1289 | register int cnt; \ | |
1290 | if (quot == 0) \ | |
1291 | cnt = BITS_PER_MP_LIMB; \ | |
1292 | else \ | |
1293 | count_leading_zeros (cnt, quot); \ | |
1294 | exponent -= cnt; \ | |
1295 | if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \ | |
1296 | { \ | |
1297 | used = MANT_DIG + cnt; \ | |
1298 | retval[0] = quot >> (BITS_PER_MP_LIMB - used); \ | |
1299 | bits = MANT_DIG + 1; \ | |
1300 | } \ | |
1301 | else \ | |
1302 | { \ | |
1303 | /* Note that we only clear the second element. */ \ | |
1304 | /* The conditional is determined at compile time. */ \ | |
1305 | if (RETURN_LIMB_SIZE > 1) \ | |
1306 | retval[1] = 0; \ | |
1307 | retval[0] = quot; \ | |
1308 | bits = -cnt; \ | |
1309 | } \ | |
1310 | } \ | |
1311 | else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \ | |
1312 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \ | |
1313 | quot); \ | |
1314 | else \ | |
1315 | { \ | |
1316 | used = MANT_DIG - bits; \ | |
1317 | if (used > 0) \ | |
1318 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \ | |
1319 | } \ | |
1320 | bits += BITS_PER_MP_LIMB | |
1321 | ||
1322 | got_limb; | |
1323 | } | |
1324 | while (bits <= MANT_DIG); | |
1325 | ||
1326 | return round_and_return (retval, exponent - 1, negative, | |
1327 | quot, BITS_PER_MP_LIMB - 1 - used, | |
1328 | more_bits || n != 0); | |
1329 | } | |
1330 | case 2: | |
1331 | { | |
1332 | mp_limb_t d0, d1, n0, n1; | |
1333 | mp_limb_t quot = 0; | |
1334 | int used = 0; | |
1335 | ||
1336 | d0 = den[0]; | |
1337 | d1 = den[1]; | |
1338 | ||
1339 | if (numsize < densize) | |
1340 | { | |
1341 | if (num[0] >= d1) | |
1342 | { | |
1343 | /* The numerator of the number occupies fewer bits than | |
1344 | the denominator but the one limb is bigger than the | |
1345 | high limb of the numerator. */ | |
1346 | n1 = 0; | |
1347 | n0 = num[0]; | |
1348 | } | |
1349 | else | |
1350 | { | |
1351 | if (bits <= 0) | |
1352 | exponent -= BITS_PER_MP_LIMB; | |
1353 | else | |
1354 | { | |
1355 | if (bits + BITS_PER_MP_LIMB <= MANT_DIG) | |
1356 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, | |
1357 | BITS_PER_MP_LIMB, 0); | |
1358 | else | |
1359 | { | |
1360 | used = MANT_DIG - bits; | |
1361 | if (used > 0) | |
1362 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); | |
1363 | } | |
1364 | bits += BITS_PER_MP_LIMB; | |
1365 | } | |
1366 | n1 = num[0]; | |
1367 | n0 = 0; | |
1368 | } | |
1369 | } | |
1370 | else | |
1371 | { | |
1372 | n1 = num[1]; | |
1373 | n0 = num[0]; | |
1374 | } | |
1375 | ||
1376 | while (bits <= MANT_DIG) | |
1377 | { | |
1378 | mp_limb_t r; | |
1379 | ||
1380 | if (n1 == d1) | |
1381 | { | |
1382 | /* QUOT should be either 111..111 or 111..110. We need | |
1383 | special treatment of this rare case as normal division | |
1384 | would give overflow. */ | |
1385 | quot = ~(mp_limb_t) 0; | |
1386 | ||
1387 | r = n0 + d1; | |
1388 | if (r < d1) /* Carry in the addition? */ | |
1389 | { | |
1390 | add_ssaaaa (n1, n0, r - d0, 0, 0, d0); | |
1391 | goto have_quot; | |
1392 | } | |
1393 | n1 = d0 - (d0 != 0); | |
1394 | n0 = -d0; | |
1395 | } | |
1396 | else | |
1397 | { | |
1398 | udiv_qrnnd (quot, r, n1, n0, d1); | |
1399 | umul_ppmm (n1, n0, d0, quot); | |
1400 | } | |
1401 | ||
1402 | q_test: | |
1403 | if (n1 > r || (n1 == r && n0 > 0)) | |
1404 | { | |
1405 | /* The estimated QUOT was too large. */ | |
1406 | --quot; | |
1407 | ||
1408 | sub_ddmmss (n1, n0, n1, n0, 0, d0); | |
1409 | r += d1; | |
1410 | if (r >= d1) /* If not carry, test QUOT again. */ | |
1411 | goto q_test; | |
1412 | } | |
1413 | sub_ddmmss (n1, n0, r, 0, n1, n0); | |
1414 | ||
1415 | have_quot: | |
1416 | got_limb; | |
1417 | } | |
1418 | ||
1419 | return round_and_return (retval, exponent - 1, negative, | |
1420 | quot, BITS_PER_MP_LIMB - 1 - used, | |
1421 | more_bits || n1 != 0 || n0 != 0); | |
1422 | } | |
1423 | default: | |
1424 | { | |
1425 | int i; | |
1426 | mp_limb_t cy, dX, d1, n0, n1; | |
1427 | mp_limb_t quot = 0; | |
1428 | int used = 0; | |
1429 | ||
1430 | dX = den[densize - 1]; | |
1431 | d1 = den[densize - 2]; | |
1432 | ||
1433 | /* The division does not work if the upper limb of the two-limb | |
1434 | numerator is greater than the denominator. */ | |
1435 | if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0) | |
1436 | num[numsize++] = 0; | |
1437 | ||
1438 | if (numsize < densize) | |
1439 | { | |
1440 | mp_size_t empty = densize - numsize; | |
1441 | ||
1442 | if (bits <= 0) | |
1443 | { | |
1444 | register int i; | |
1445 | for (i = numsize; i > 0; --i) | |
1446 | num[i + empty] = num[i - 1]; | |
1447 | MPN_ZERO (num, empty + 1); | |
1448 | exponent -= empty * BITS_PER_MP_LIMB; | |
1449 | } | |
1450 | else | |
1451 | { | |
1452 | if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG) | |
1453 | { | |
1454 | /* We make a difference here because the compiler | |
1455 | cannot optimize the `else' case that good and | |
1456 | this reflects all currently used FLOAT types | |
1457 | and GMP implementations. */ | |
1458 | register int i; | |
1459 | #if RETURN_LIMB_SIZE <= 2 | |
1460 | assert (empty == 1); | |
1461 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, | |
1462 | BITS_PER_MP_LIMB, 0); | |
1463 | #else | |
1464 | for (i = RETURN_LIMB_SIZE; i > empty; --i) | |
1465 | retval[i] = retval[i - empty]; | |
1466 | #endif | |
1467 | for (i = numsize; i > 0; --i) | |
1468 | num[i + empty] = num[i - 1]; | |
1469 | MPN_ZERO (num, empty + 1); | |
1470 | } | |
1471 | else | |
1472 | { | |
1473 | used = MANT_DIG - bits; | |
1474 | if (used >= BITS_PER_MP_LIMB) | |
1475 | { | |
1476 | register int i; | |
1477 | (void) __mpn_lshift (&retval[used | |
1478 | / BITS_PER_MP_LIMB], | |
1479 | retval, RETURN_LIMB_SIZE, | |
1480 | used % BITS_PER_MP_LIMB); | |
1481 | for (i = used / BITS_PER_MP_LIMB; i >= 0; --i) | |
1482 | retval[i] = 0; | |
1483 | } | |
1484 | else if (used > 0) | |
1485 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); | |
1486 | } | |
1487 | bits += empty * BITS_PER_MP_LIMB; | |
1488 | } | |
1489 | } | |
1490 | else | |
1491 | { | |
1492 | int i; | |
1493 | assert (numsize == densize); | |
1494 | for (i = numsize; i > 0; --i) | |
1495 | num[i] = num[i - 1]; | |
1496 | } | |
1497 | ||
1498 | den[densize] = 0; | |
1499 | n0 = num[densize]; | |
1500 | ||
1501 | while (bits <= MANT_DIG) | |
1502 | { | |
1503 | if (n0 == dX) | |
1504 | /* This might over-estimate QUOT, but it's probably not | |
1505 | worth the extra code here to find out. */ | |
1506 | quot = ~(mp_limb_t) 0; | |
1507 | else | |
1508 | { | |
1509 | mp_limb_t r; | |
1510 | ||
1511 | udiv_qrnnd (quot, r, n0, num[densize - 1], dX); | |
1512 | umul_ppmm (n1, n0, d1, quot); | |
1513 | ||
1514 | while (n1 > r || (n1 == r && n0 > num[densize - 2])) | |
1515 | { | |
1516 | --quot; | |
1517 | r += dX; | |
1518 | if (r < dX) /* I.e. "carry in previous addition?" */ | |
1519 | break; | |
1520 | n1 -= n0 < d1; | |
1521 | n0 -= d1; | |
1522 | } | |
1523 | } | |
1524 | ||
1525 | /* Possible optimization: We already have (q * n0) and (1 * n1) | |
1526 | after the calculation of QUOT. Taking advantage of this, we | |
1527 | could make this loop make two iterations less. */ | |
1528 | ||
1529 | cy = __mpn_submul_1 (num, den, densize + 1, quot); | |
1530 | ||
1531 | if (num[densize] != cy) | |
1532 | { | |
1533 | cy = __mpn_add_n (num, num, den, densize); | |
1534 | assert (cy != 0); | |
1535 | --quot; | |
1536 | } | |
1537 | n0 = num[densize] = num[densize - 1]; | |
1538 | for (i = densize - 1; i > 0; --i) | |
1539 | num[i] = num[i - 1]; | |
1540 | ||
1541 | got_limb; | |
1542 | } | |
1543 | ||
1544 | for (i = densize; num[i] == 0 && i >= 0; --i) | |
1545 | ; | |
1546 | return round_and_return (retval, exponent - 1, negative, | |
1547 | quot, BITS_PER_MP_LIMB - 1 - used, | |
1548 | more_bits || i >= 0); | |
1549 | } | |
1550 | } | |
1551 | } | |
1552 | ||
1553 | /* NOTREACHED */ | |
1554 | } | |
1555 | #if defined _LIBC && !defined USE_WIDE_CHAR | |
1556 | libc_hidden_def (INTERNAL (__STRTOF)) | |
1557 | #endif | |
1558 | \f | |
1559 | /* External user entry point. */ | |
1ab62b32 | 1560 | |
ccadf7b5 UD |
1561 | FLOAT |
1562 | #ifdef weak_function | |
1563 | weak_function | |
1564 | #endif | |
1565 | __STRTOF (nptr, endptr, loc) | |
1566 | const STRING_TYPE *nptr; | |
1567 | STRING_TYPE **endptr; | |
1568 | __locale_t loc; | |
1569 | { | |
1570 | return INTERNAL (__STRTOF) (nptr, endptr, 0, loc); | |
1571 | } | |
1572 | weak_alias (__STRTOF, STRTOF) |