]>
Commit | Line | Data |
---|---|---|
28f540f4 RM |
1 | /* Read decimal floating point numbers. |
2 | Copyright (C) 1995 Free Software Foundation, Inc. | |
3 | Contributed by Ulrich Drepper. | |
4 | ||
5 | This file is part of the GNU C Library. | |
6 | ||
7 | The GNU C Library is free software; you can redistribute it and/or | |
8 | modify it under the terms of the GNU Library General Public License as | |
9 | published by the Free Software Foundation; either version 2 of the | |
10 | License, or (at your option) any later version. | |
11 | ||
12 | The GNU C Library is distributed in the hope that it will be useful, | |
13 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 | Library General Public License for more details. | |
16 | ||
17 | You should have received a copy of the GNU Library General Public | |
18 | License along with the GNU C Library; see the file COPYING.LIB. If | |
19 | not, write to the Free Software Foundation, Inc., 675 Mass Ave, | |
20 | Cambridge, MA 02139, USA. */ | |
21 | ||
22 | /* Configuration part. These macros are defined by `strtold.c' and `strtof.c' | |
23 | to produce the `long double' and `float' versions of the reader. */ | |
24 | #ifndef FLOAT | |
25 | #define FLOAT double | |
26 | #define FLT DBL | |
27 | #define STRTOF strtod | |
28 | #define MPN2FLOAT __mpn_construct_double | |
29 | #define FLOAT_HUGE_VAL HUGE_VAL | |
4933a099 | 30 | #define IMPLICIT_ONE 1 |
28f540f4 RM |
31 | #endif |
32 | /* End of configuration part. */ | |
33 | \f | |
34 | #include <ctype.h> | |
35 | #include <errno.h> | |
36 | #include <float.h> | |
933e73fa | 37 | #include "../locale/localeinfo.h" |
28f540f4 RM |
38 | #include <math.h> |
39 | #include <stdlib.h> | |
933e73fa RM |
40 | #include "../stdio/gmp.h" |
41 | #include "../stdio/gmp-impl.h" | |
28f540f4 | 42 | #include <gmp-mparam.h> |
933e73fa RM |
43 | #include "../stdio/longlong.h" |
44 | #include "../stdio/fpioconst.h" | |
28f540f4 | 45 | |
4933a099 | 46 | #define NDEBUG 1 |
28f540f4 RM |
47 | #include <assert.h> |
48 | ||
49 | ||
50 | /* Constants we need from float.h; select the set for the FLOAT precision. */ | |
0923c7a5 | 51 | #define MANT_DIG PASTE(FLT,_MANT_DIG) |
4933a099 | 52 | #define DIG PASTE(FLT,_DIG) |
0923c7a5 RM |
53 | #define MAX_EXP PASTE(FLT,_MAX_EXP) |
54 | #define MIN_EXP PASTE(FLT,_MIN_EXP) | |
55 | #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP) | |
56 | #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP) | |
57 | #define MAX_10_EXP_LOG PASTE(FLT,_MAX_10_EXP_LOG) | |
28f540f4 | 58 | |
0923c7a5 RM |
59 | /* Extra macros required to get FLT expanded before the pasting. */ |
60 | #define PASTE(a,b) PASTE1(a,b) | |
61 | #define PASTE1(a,b) a##b | |
28f540f4 RM |
62 | |
63 | /* Function to construct a floating point number from an MP integer | |
64 | containing the fraction bits, a base 2 exponent, and a sign flag. */ | |
65 | extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative); | |
66 | \f | |
67 | /* Definitions according to limb size used. */ | |
68 | #if BITS_PER_MP_LIMB == 32 | |
69 | # define MAX_DIG_PER_LIMB 9 | |
70 | # define MAX_FAC_PER_LIMB 1000000000L | |
71 | #elif BITS_PER_MP_LIMB == 64 | |
72 | # define MAX_DIG_PER_LIMB 19 | |
73 | # define MAX_FAC_PER_LIMB 10000000000000000000L | |
74 | #else | |
75 | # error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" | |
76 | #endif | |
77 | ||
78 | ||
79 | /* Local data structure. */ | |
4933a099 RM |
80 | static const mp_limb _tens_in_limb[MAX_DIG_PER_LIMB + 1] = |
81 | { 0, 10, 100, | |
82 | 1000, 10000, 100000, | |
83 | 1000000, 10000000, 100000000, | |
84 | 1000000000 | |
28f540f4 | 85 | #if BITS_PER_MP_LIMB > 32 |
4933a099 RM |
86 | , 10000000000, 100000000000, |
87 | 1000000000000, 10000000000000, 100000000000000, | |
88 | 1000000000000000, 10000000000000000, 100000000000000000, | |
89 | 1000000000000000000, 10000000000000000000 | |
28f540f4 RM |
90 | #endif |
91 | #if BITS_PER_MP_LIMB > 64 | |
92 | #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB | |
93 | #endif | |
94 | }; | |
95 | \f | |
96 | #ifndef howmany | |
97 | #define howmany(x,y) (((x)+((y)-1))/(y)) | |
98 | #endif | |
99 | #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; }) | |
100 | ||
101 | #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG) | |
102 | #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB) | |
103 | ||
104 | #define RETURN(val,end) \ | |
105 | do { if (endptr != 0) *endptr = (char *) end; return val; } while (0) | |
106 | ||
107 | /* Maximum size necessary for mpn integers to hold floating point numbers. */ | |
4933a099 | 108 | #define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) + 2) |
28f540f4 RM |
109 | /* Declare an mpn integer variable that big. */ |
110 | #define MPN_VAR(name) mp_limb name[MPNSIZE]; mp_size_t name##size | |
111 | /* Copy an mpn integer value. */ | |
112 | #define MPN_ASSIGN(dst, src) \ | |
113 | memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb)) | |
114 | ||
115 | ||
116 | /* Return a floating point number of the needed type according to the given | |
117 | multi-precision number after possible rounding. */ | |
118 | static inline FLOAT | |
119 | round_and_return (mp_limb *retval, int exponent, int negative, | |
120 | mp_limb round_limb, mp_size_t round_bit, int more_bits) | |
121 | { | |
4933a099 | 122 | if (exponent < MIN_EXP - 2 + IMPLICIT_ONE) |
28f540f4 | 123 | { |
4933a099 | 124 | mp_size_t shift = MIN_EXP - 2 + IMPLICIT_ONE - exponent; |
28f540f4 RM |
125 | |
126 | if (shift >= MANT_DIG) | |
127 | { | |
128 | errno = EDOM; | |
129 | return 0.0; | |
130 | } | |
131 | ||
132 | more_bits |= (round_limb & ((1 << round_bit) - 1)) != 0; | |
133 | if (shift >= BITS_PER_MP_LIMB) | |
134 | { | |
135 | round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB]; | |
136 | round_bit = (shift - 1) % BITS_PER_MP_LIMB; | |
4933a099 RM |
137 | |
138 | (void) __mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB], | |
139 | RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB), | |
140 | shift % BITS_PER_MP_LIMB); | |
141 | MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)], | |
142 | shift / BITS_PER_MP_LIMB); | |
28f540f4 | 143 | } |
4933a099 | 144 | else if (shift > 0) |
28f540f4 RM |
145 | { |
146 | round_limb = retval[0]; | |
147 | round_bit = shift - 1; | |
4933a099 | 148 | (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift); |
28f540f4 | 149 | } |
28f540f4 RM |
150 | exponent = MIN_EXP - 2; |
151 | } | |
152 | ||
4933a099 RM |
153 | if ((round_limb & (1 << round_bit)) != 0 |
154 | && (more_bits || (retval[0] & 1) != 0 | |
155 | || (round_limb & ((1 << round_bit) - 1)) != 0)) | |
28f540f4 RM |
156 | { |
157 | mp_limb cy = __mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1); | |
4933a099 RM |
158 | |
159 | if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || | |
160 | ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && | |
161 | (retval[RETURN_LIMB_SIZE - 1] | |
162 | & (1 << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)) | |
28f540f4 RM |
163 | { |
164 | ++exponent; | |
165 | (void) __mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1); | |
4933a099 RM |
166 | retval[RETURN_LIMB_SIZE - 1] |= 1 << ((MANT_DIG - 1) |
167 | % BITS_PER_MP_LIMB); | |
28f540f4 | 168 | } |
4933a099 RM |
169 | else if (IMPLICIT_ONE && exponent == MIN_EXP - 2 |
170 | && (retval[RETURN_LIMB_SIZE - 1] | |
171 | & (1 << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) != 0) | |
172 | /* The number was denormalized but now normalized. */ | |
173 | exponent = MIN_EXP - 1; | |
28f540f4 RM |
174 | } |
175 | ||
176 | if (exponent > MAX_EXP) | |
177 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
178 | ||
179 | return MPN2FLOAT (retval, exponent, negative); | |
180 | } | |
181 | ||
182 | ||
183 | /* Read a multi-precision integer starting at STR with exactly DIGCNT digits | |
184 | into N. Return the size of the number limbs in NSIZE at the first | |
185 | character od the string that is not part of the integer as the function | |
186 | value. If the EXPONENT is small enough to be taken as an additional | |
187 | factor for the resulting number (see code) multiply by it. */ | |
188 | static inline const char * | |
189 | str_to_mpn (const char *str, int digcnt, mp_limb *n, mp_size_t *nsize, | |
190 | int *exponent) | |
191 | { | |
192 | /* Number of digits for actual limb. */ | |
193 | int cnt = 0; | |
194 | mp_limb low = 0; | |
195 | mp_limb base; | |
196 | ||
197 | *nsize = 0; | |
198 | assert (digcnt > 0); | |
199 | do | |
200 | { | |
201 | if (cnt == MAX_DIG_PER_LIMB) | |
202 | { | |
203 | if (*nsize == 0) | |
204 | n[0] = low; | |
205 | else | |
206 | { | |
207 | mp_limb cy; | |
208 | cy = __mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB); | |
209 | cy += __mpn_add_1 (n, n, *nsize, low); | |
210 | if (cy != 0) | |
211 | n[*nsize] = cy; | |
212 | } | |
213 | ++(*nsize); | |
214 | cnt = 0; | |
215 | low = 0; | |
216 | } | |
217 | ||
218 | /* There might be thousands separators or radix characters in the string. | |
219 | But these all can be ignored because we know the format of the number | |
220 | is correct and we have an exact number of characters to read. */ | |
221 | while (!isdigit (*str)) | |
222 | ++str; | |
223 | low = low * 10 + *str++ - '0'; | |
224 | ++cnt; | |
225 | } | |
226 | while (--digcnt > 0); | |
227 | ||
228 | if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB) | |
229 | { | |
230 | low *= _tens_in_limb[*exponent]; | |
231 | base = _tens_in_limb[cnt + *exponent]; | |
232 | *exponent = 0; | |
233 | } | |
234 | else | |
235 | base = _tens_in_limb[cnt]; | |
236 | ||
237 | if (*nsize == 0) | |
238 | { | |
239 | n[0] = low; | |
240 | *nsize = 1; | |
241 | } | |
242 | else | |
243 | { | |
244 | mp_limb cy; | |
245 | cy = __mpn_mul_1 (n, n, *nsize, base); | |
246 | cy += __mpn_add_1 (n, n, *nsize, low); | |
247 | if (cy != 0) | |
248 | n[(*nsize)++] = cy; | |
249 | } | |
250 | return str; | |
251 | } | |
252 | ||
253 | ||
254 | /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits | |
255 | with the COUNT most significant bits of LIMB. | |
256 | ||
257 | Tege doesn't like this function so I have to write it here myself. :) | |
258 | --drepper */ | |
259 | static inline void | |
260 | __mpn_lshift_1 (mp_limb *ptr, mp_size_t size, unsigned int count, mp_limb limb) | |
261 | { | |
262 | if (count == BITS_PER_MP_LIMB) | |
263 | { | |
264 | /* Optimize the case of shifting by exactly a word: | |
265 | just copy words, with no actual bit-shifting. */ | |
266 | mp_size_t i; | |
267 | for (i = size - 1; i > 0; --i) | |
268 | ptr[i] = ptr[i - 1]; | |
269 | ptr[0] = limb; | |
270 | } | |
271 | else | |
272 | { | |
273 | (void) __mpn_lshift (ptr, ptr, size, count); | |
274 | ptr[0] |= limb >> (BITS_PER_MP_LIMB - count); | |
275 | } | |
276 | } | |
277 | ||
278 | ||
279 | /* Return a floating point number with the value of the given string NPTR. | |
280 | Set *ENDPTR to the character after the last used one. If the number is | |
281 | smaller than the smallest representable number, set `errno' to ERANGE and | |
282 | return 0.0. If the number is too big to be represented, set `errno' to | |
283 | ERANGE and return HUGE_VAL with the approriate sign. */ | |
284 | FLOAT | |
285 | STRTOF (nptr, endptr) | |
286 | const char *nptr; | |
287 | char **endptr; | |
288 | { | |
289 | int negative; /* The sign of the number. */ | |
290 | MPN_VAR (num); /* MP representation of the number. */ | |
291 | int exponent; /* Exponent of the number. */ | |
292 | ||
293 | /* When we have to compute fractional digits we form a fraction with a | |
294 | second multi-precision number (and we sometimes need a second for | |
295 | temporary results). */ | |
296 | MPN_VAR (den); | |
297 | ||
298 | /* Representation for the return value. */ | |
299 | mp_limb retval[RETURN_LIMB_SIZE]; | |
300 | /* Number of bits currently in result value. */ | |
301 | int bits; | |
302 | ||
303 | /* Running pointer after the last character processed in the string. */ | |
304 | const char *cp; | |
305 | /* Start of significant part of the number. */ | |
306 | const char *startp; | |
307 | /* Points at the character following the integer and fractional digits. */ | |
308 | const char *expp; | |
309 | /* Total number of digit and number of digits in integer part. */ | |
4933a099 | 310 | int dig_no, int_no, lead_zero; |
28f540f4 RM |
311 | /* Contains the last character read. */ |
312 | char c; | |
313 | ||
314 | /* The radix character of the current locale. */ | |
315 | wchar_t decimal; | |
316 | #ifdef USE_GROUPING | |
317 | /* The thousands character of the current locale. */ | |
318 | wchar_t thousands; | |
319 | /* The numeric grouping specification of the current locale, | |
320 | in the format described in <locale.h>. */ | |
321 | const char *grouping; | |
322 | ||
323 | /* Check the grouping of the integer part at [BEGIN,END). | |
324 | Return zero iff a separator is found out of place. */ | |
325 | int grouping_ok (const char *begin, const char *end) | |
326 | { | |
327 | if (grouping) | |
328 | while (end > begin) | |
329 | { | |
330 | const char *p = end; | |
331 | do | |
332 | --p; | |
333 | while (*p != thousands && p > begin); | |
334 | if (end - 1 - p != *grouping++) | |
335 | return 0; /* Wrong number of digits in this group. */ | |
336 | end = p; /* Correct group; trim it off the end. */ | |
337 | ||
338 | if (*grouping == 0) | |
339 | --grouping; /* Same grouping repeats in next iteration. */ | |
340 | else if (*grouping == CHAR_MAX || *grouping < 0) | |
341 | { | |
342 | /* No further grouping allowed. */ | |
343 | while (end > begin) | |
344 | if (*--end == thousands) | |
345 | return 0; | |
346 | } | |
347 | } | |
348 | return 1; | |
349 | } | |
350 | /* Return with no conversion if the grouping of [STARTP,CP) is bad. */ | |
351 | #define CHECK_GROUPING if (! grouping_ok (startp, cp)) RETURN (0.0, nptr); else | |
352 | ||
933e73fa | 353 | grouping = _NL_CURRENT (LC_NUMERIC, GROUPING); |
28f540f4 RM |
354 | if (*grouping <= 0 || *grouping == CHAR_MAX) |
355 | grouping = NULL; | |
356 | else | |
357 | { | |
358 | /* Figure out the thousands seperator character. */ | |
933e73fa RM |
359 | if (mbtowc (&thousands_sep, _NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP), |
360 | strlen (_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP))) <= 0) | |
361 | thousands = (wchar_t) *_NL_CURRENT (LC_NUMERIC, THOUSANDS_SEP); | |
28f540f4 RM |
362 | if (thousands == L'\0') |
363 | grouping = NULL; | |
364 | } | |
365 | #else | |
366 | #define grouping NULL | |
367 | #define thousands L'\0' | |
368 | #define CHECK_GROUPING ((void) 0) | |
369 | #endif | |
370 | ||
371 | /* Find the locale's decimal point character. */ | |
933e73fa RM |
372 | if (mbtowc (&decimal, _NL_CURRENT (LC_NUMERIC, DECIMAL_POINT), |
373 | strlen (_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT))) <= 0) | |
374 | decimal = (wchar_t) *_NL_CURRENT (LC_NUMERIC, DECIMAL_POINT); | |
28f540f4 RM |
375 | |
376 | ||
377 | /* Prepare number representation. */ | |
378 | exponent = 0; | |
379 | negative = 0; | |
380 | bits = 0; | |
381 | ||
382 | /* Parse string to get maximal legal prefix. We need the number of | |
383 | characters of the interger part, the fractional part and the exponent. */ | |
384 | cp = nptr - 1; | |
385 | /* Ignore leading white space. */ | |
386 | do | |
387 | c = *++cp; | |
388 | while (isspace (c)); | |
389 | ||
390 | /* Get sign of the result. */ | |
391 | if (c == '-') | |
392 | { | |
393 | negative = 1; | |
394 | c = *++cp; | |
395 | } | |
396 | else if (c == '+') | |
397 | c = *++cp; | |
398 | ||
399 | /* Return 0.0 if no legal string is found. | |
400 | No character is used even if a sign was found. */ | |
37f91d33 | 401 | if (!isdigit (c) && (c != decimal || !isdigit (cp[1]))) |
28f540f4 RM |
402 | RETURN (0.0, nptr); |
403 | ||
404 | /* Record the start of the digits, in case we will check their grouping. */ | |
405 | startp = cp; | |
406 | ||
407 | /* Ignore leading zeroes. This helps us to avoid useless computations. */ | |
408 | while (c == '0' || (thousands != L'\0' && c == thousands)) | |
409 | c = *++cp; | |
410 | ||
411 | CHECK_GROUPING; | |
412 | ||
413 | /* If no other digit but a '0' is found the result is 0.0. | |
414 | Return current read pointer. */ | |
415 | if (!isdigit (c) && c != decimal) | |
416 | RETURN (0.0, cp); | |
417 | ||
418 | /* Remember first significant digit and read following characters until the | |
419 | decimal point, exponent character or any non-FP number character. */ | |
420 | startp = cp; | |
421 | dig_no = 0; | |
422 | while (dig_no < NDIG || | |
423 | /* If parsing grouping info, keep going past useful digits | |
424 | so we can check all the grouping separators. */ | |
425 | grouping) | |
426 | { | |
427 | if (isdigit (c)) | |
428 | ++dig_no; | |
429 | else if (thousands == L'\0' || c != thousands) | |
430 | /* Not a digit or separator: end of the integer part. */ | |
431 | break; | |
432 | c = *++cp; | |
433 | } | |
434 | ||
435 | CHECK_GROUPING; | |
436 | ||
437 | if (dig_no >= NDIG) | |
438 | /* Too many digits to be representable. Assigning this to EXPONENT | |
439 | allows us to read the full number but return HUGE_VAL after parsing. */ | |
440 | exponent = MAX_10_EXP; | |
441 | ||
442 | /* We have the number digits in the integer part. Whether these are all or | |
443 | any is really a fractional digit will be decided later. */ | |
444 | int_no = dig_no; | |
4933a099 | 445 | lead_zero = int_no == 0 ? -1 : 0; |
28f540f4 RM |
446 | |
447 | /* Read the fractional digits. */ | |
448 | if (c == decimal) | |
449 | { | |
450 | if (isdigit (cp[1])) | |
451 | { | |
4933a099 | 452 | c = *++cp; |
28f540f4 RM |
453 | do |
454 | { | |
4933a099 RM |
455 | if (c != '0' && lead_zero == -1) |
456 | lead_zero = dig_no - int_no; | |
28f540f4 RM |
457 | ++dig_no; |
458 | c = *++cp; | |
459 | } | |
460 | while (isdigit (c)); | |
461 | } | |
462 | } | |
463 | ||
464 | /* Remember start of exponent (if any). */ | |
465 | expp = cp; | |
466 | ||
467 | /* Read exponent. */ | |
468 | if (tolower (c) == 'e') | |
469 | { | |
470 | int exp_negative = 0; | |
471 | ||
472 | c = *++cp; | |
473 | if (c == '-') | |
474 | { | |
475 | exp_negative = 1; | |
476 | c = *++cp; | |
477 | } | |
478 | else if (c == '+') | |
479 | c = *++cp; | |
480 | ||
481 | if (isdigit (c)) | |
482 | { | |
483 | do | |
484 | { | |
485 | if ((!exp_negative && exponent * 10 + int_no > MAX_10_EXP) | |
486 | || (exp_negative | |
487 | && exponent * 10 + int_no > -MIN_10_EXP + MANT_DIG)) | |
488 | /* The exponent is too large/small to represent a valid | |
489 | number. */ | |
490 | { | |
491 | FLOAT retval; | |
492 | ||
493 | /* Overflow or underflow. */ | |
494 | errno = ERANGE; | |
495 | retval = (exp_negative ? 0.0 : | |
496 | negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL); | |
497 | ||
498 | /* Accept all following digits as part of the exponent. */ | |
499 | do | |
500 | ++cp; | |
501 | while (isdigit (*cp)); | |
502 | ||
503 | RETURN (retval, cp); | |
504 | /* NOTREACHED */ | |
505 | } | |
506 | ||
507 | exponent *= 10; | |
508 | exponent += c - '0'; | |
509 | c = *++cp; | |
510 | } | |
511 | while (isdigit (c)); | |
512 | } | |
513 | else | |
514 | cp = expp; | |
515 | ||
516 | if (exp_negative) | |
517 | exponent = -exponent; | |
518 | } | |
519 | ||
520 | /* We don't want to have to work with trailing zeroes after the radix. */ | |
521 | if (dig_no > int_no) | |
522 | { | |
523 | while (expp[-1] == '0') | |
524 | { | |
525 | --expp; | |
526 | --dig_no; | |
527 | } | |
528 | assert (dig_no >= int_no); | |
529 | } | |
530 | ||
531 | /* The whole string is parsed. Store the address of the next character. */ | |
532 | if (endptr) | |
533 | *endptr = (char *) cp; | |
534 | ||
535 | if (dig_no == 0) | |
536 | return 0.0; | |
537 | ||
538 | /* Now we have the number of digits in total and the integer digits as well | |
539 | as the exponent and its sign. We can decide whether the read digits are | |
540 | really integer digits or belong to the fractional part; i.e. we normalize | |
541 | 123e-2 to 1.23. */ | |
542 | { | |
543 | register int incr = exponent < 0 ? MAX (-int_no, exponent) | |
544 | : MIN (dig_no - int_no, exponent); | |
545 | int_no += incr; | |
546 | exponent -= incr; | |
547 | } | |
548 | ||
549 | if (int_no + exponent > MAX_10_EXP) | |
550 | { | |
551 | errno = ERANGE; | |
552 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; | |
553 | } | |
554 | ||
4933a099 | 555 | if (exponent - MAX(0, lead_zero) < MIN_10_EXP - (DIG + 1)) |
28f540f4 RM |
556 | { |
557 | errno = ERANGE; | |
558 | return 0.0; | |
559 | } | |
560 | ||
561 | if (int_no > 0) | |
562 | { | |
563 | /* Read the integer part as a multi-precision number to NUM. */ | |
564 | startp = str_to_mpn (startp, int_no, num, &numsize, &exponent); | |
565 | ||
566 | if (exponent > 0) | |
567 | { | |
568 | /* We now multiply the gained number by the given power of ten. */ | |
569 | mp_limb *psrc = num; | |
570 | mp_limb *pdest = den; | |
571 | int expbit = 1; | |
572 | const struct mp_power *ttab = &_fpioconst_pow10[0]; | |
573 | ||
574 | assert (exponent < (1 << (MAX_10_EXP_LOG + 1))); | |
575 | do | |
576 | { | |
577 | if ((exponent & expbit) != 0) | |
578 | { | |
579 | mp_limb cy; | |
580 | exponent ^= expbit; | |
581 | ||
582 | /* FIXME: not the whole multiplication has to be done. | |
583 | If we have the needed number of bits we only need the | |
584 | information whether more non-zero bits follow. */ | |
585 | if (numsize >= ttab->arraysize - 2) | |
586 | cy = __mpn_mul (pdest, psrc, numsize, | |
587 | &ttab->array[2], ttab->arraysize - 2); | |
588 | else | |
589 | cy = __mpn_mul (pdest, &ttab->array[2], | |
590 | ttab->arraysize - 2, | |
591 | psrc, numsize); | |
592 | numsize += ttab->arraysize - 2; | |
593 | if (cy == 0) | |
594 | --numsize; | |
595 | SWAP (psrc, pdest); | |
596 | } | |
597 | expbit <<= 1; | |
598 | ++ttab; | |
599 | } | |
600 | while (exponent != 0); | |
601 | ||
602 | if (psrc == den) | |
603 | memcpy (num, den, numsize * sizeof (mp_limb)); | |
604 | } | |
605 | ||
606 | /* Determine how many bits of the result we already have. */ | |
607 | count_leading_zeros (bits, num[numsize - 1]); | |
608 | bits = numsize * BITS_PER_MP_LIMB - bits; | |
609 | ||
610 | /* We have already the first BITS bits of the result. Together with | |
611 | the information whether more non-zero bits follow this is enough | |
612 | to determine the result. */ | |
613 | if (bits > MANT_DIG) | |
614 | { | |
4933a099 | 615 | int i; |
28f540f4 RM |
616 | const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB; |
617 | const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB; | |
618 | const mp_size_t round_idx = least_bit == 0 ? least_idx - 1 | |
619 | : least_idx; | |
620 | const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1 | |
4933a099 | 621 | : least_bit - 1; |
28f540f4 RM |
622 | |
623 | if (least_bit == 0) | |
624 | memcpy (retval, &num[least_idx], | |
625 | RETURN_LIMB_SIZE * sizeof (mp_limb)); | |
626 | else | |
4933a099 RM |
627 | { |
628 | for (i = least_idx; i < numsize - 1; ++i) | |
629 | retval[i - least_idx] = (num[i] >> least_bit) | |
630 | | (num[i + 1] | |
631 | << (BITS_PER_MP_LIMB - least_bit)); | |
632 | if (i - least_idx < RETURN_LIMB_SIZE) | |
633 | retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit; | |
634 | } | |
28f540f4 RM |
635 | |
636 | /* Check whether any limb beside the ones in RETVAL are non-zero. */ | |
637 | for (i = 0; num[i] == 0; ++i) | |
638 | ; | |
639 | ||
640 | return round_and_return (retval, bits - 1, negative, | |
641 | num[round_idx], round_bit, | |
642 | int_no < dig_no || i < round_idx); | |
643 | /* NOTREACHED */ | |
644 | } | |
645 | else if (dig_no == int_no) | |
646 | { | |
647 | const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; | |
648 | const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB; | |
649 | ||
650 | if (target_bit == is_bit) | |
651 | { | |
652 | memcpy (&retval[RETURN_LIMB_SIZE - numsize], num, | |
653 | numsize * sizeof (mp_limb)); | |
654 | /* FIXME: the following loop can be avoided if we assume a | |
655 | maximal MANT_DIG value. */ | |
656 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); | |
657 | } | |
658 | else if (target_bit > is_bit) | |
659 | { | |
660 | (void) __mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize], | |
661 | num, numsize, target_bit - is_bit); | |
662 | /* FIXME: the following loop can be avoided if we assume a | |
663 | maximal MANT_DIG value. */ | |
664 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); | |
665 | } | |
666 | else | |
667 | { | |
668 | mp_limb cy; | |
669 | assert (numsize < RETURN_LIMB_SIZE); | |
670 | ||
671 | cy = __mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize], | |
672 | num, numsize, is_bit - target_bit); | |
673 | retval[RETURN_LIMB_SIZE - numsize - 1] = cy; | |
674 | /* FIXME: the following loop can be avoided if we assume a | |
675 | maximal MANT_DIG value. */ | |
676 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1); | |
677 | } | |
678 | ||
679 | return round_and_return (retval, bits - 1, negative, 0, 0, 0); | |
680 | /* NOTREACHED */ | |
681 | } | |
682 | ||
683 | /* Store the bits we already have. */ | |
684 | memcpy (retval, num, numsize * sizeof (mp_limb)); | |
685 | #if RETURN_LIMB_SIZE > 1 | |
686 | if (numsize < RETURN_LIMB_SIZE) | |
687 | retval[numsize] = 0; | |
688 | #endif | |
689 | } | |
690 | ||
691 | /* We have to compute at least some of the fractional digits. */ | |
692 | { | |
693 | /* We construct a fraction and the result of the division gives us | |
694 | the needed digits. The denominator is 1.0 multiplied by the | |
695 | exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and | |
696 | 123e6 gives 123 / 1000000. */ | |
697 | ||
698 | int expbit; | |
699 | int cnt; | |
4933a099 RM |
700 | int neg_exp; |
701 | int more_bits; | |
28f540f4 RM |
702 | mp_limb cy; |
703 | mp_limb *psrc = den; | |
704 | mp_limb *pdest = num; | |
28f540f4 RM |
705 | const struct mp_power *ttab = &_fpioconst_pow10[0]; |
706 | ||
707 | assert (dig_no > int_no && exponent <= 0); | |
708 | ||
4933a099 RM |
709 | |
710 | /* For the fractional part we need not process too much digits. One | |
711 | decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute | |
712 | ceil(BITS / 3) =: N | |
713 | digits we should have enough bits for the result. The remaining | |
714 | decimal digits give us the information that more bits are following. | |
715 | This can be used while rounding. (One added as a safety margin.) */ | |
716 | if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 1) | |
717 | { | |
718 | dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 1; | |
719 | more_bits = 1; | |
720 | } | |
721 | else | |
722 | more_bits = 0; | |
723 | ||
724 | neg_exp = dig_no - int_no - exponent; | |
725 | ||
28f540f4 RM |
726 | /* Construct the denominator. */ |
727 | densize = 0; | |
728 | expbit = 1; | |
729 | do | |
730 | { | |
731 | if ((neg_exp & expbit) != 0) | |
732 | { | |
733 | mp_limb cy; | |
734 | neg_exp ^= expbit; | |
735 | ||
736 | if (densize == 0) | |
737 | memcpy (psrc, &ttab->array[2], | |
738 | (densize = ttab->arraysize - 2) * sizeof (mp_limb)); | |
739 | else | |
740 | { | |
741 | cy = __mpn_mul (pdest, &ttab->array[2], ttab->arraysize - 2, | |
742 | psrc, densize); | |
743 | densize += ttab->arraysize - 2; | |
744 | if (cy == 0) | |
745 | --densize; | |
746 | SWAP (psrc, pdest); | |
747 | } | |
748 | } | |
749 | expbit <<= 1; | |
750 | ++ttab; | |
751 | } | |
752 | while (neg_exp != 0); | |
753 | ||
754 | if (psrc == num) | |
755 | memcpy (den, num, densize * sizeof (mp_limb)); | |
756 | ||
757 | /* Read the fractional digits from the string. */ | |
758 | (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent); | |
759 | ||
760 | ||
761 | /* We now have to shift both numbers so that the highest bit in the | |
762 | denominator is set. In the same process we copy the numerator to | |
763 | a high place in the array so that the division constructs the wanted | |
764 | digits. This is done by a "quasi fix point" number representation. | |
765 | ||
766 | num: ddddddddddd . 0000000000000000000000 | |
767 | |--- m ---| | |
768 | den: ddddddddddd n >= m | |
769 | |--- n ---| | |
770 | */ | |
771 | ||
772 | count_leading_zeros (cnt, den[densize - 1]); | |
773 | ||
774 | (void) __mpn_lshift (den, den, densize, cnt); | |
775 | cy = __mpn_lshift (num, num, numsize, cnt); | |
776 | if (cy != 0) | |
777 | num[numsize++] = cy; | |
778 | ||
779 | /* Now we are ready for the division. But it is not necessary to | |
780 | do a full multi-precision division because we only need a small | |
781 | number of bits for the result. So we do not use __mpn_divmod | |
782 | here but instead do the division here by hand and stop whenever | |
783 | the needed number of bits is reached. The code itself comes | |
784 | from the GNU MP Library by Torbj\"orn Granlund. */ | |
785 | ||
786 | exponent = bits; | |
787 | ||
788 | switch (densize) | |
789 | { | |
790 | case 1: | |
791 | { | |
792 | mp_limb d, n, quot; | |
793 | int used = 0; | |
794 | ||
795 | n = num[0]; | |
796 | d = den[0]; | |
797 | assert (numsize == 1 && n < d); | |
798 | ||
799 | do | |
800 | { | |
801 | udiv_qrnnd (quot, n, n, 0, d); | |
802 | ||
803 | #define got_limb \ | |
804 | if (bits == 0) \ | |
805 | { \ | |
806 | register int cnt; \ | |
807 | if (quot == 0) \ | |
808 | cnt = BITS_PER_MP_LIMB; \ | |
809 | else \ | |
810 | count_leading_zeros (cnt, quot); \ | |
811 | exponent -= cnt; \ | |
812 | if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \ | |
813 | { \ | |
4933a099 | 814 | used = MANT_DIG + cnt; \ |
28f540f4 | 815 | retval[0] = quot >> (BITS_PER_MP_LIMB - used); \ |
4933a099 | 816 | bits = MANT_DIG + 1; \ |
28f540f4 RM |
817 | } \ |
818 | else \ | |
819 | { \ | |
820 | /* Note that we only clear the second element. */ \ | |
4933a099 RM |
821 | /* The conditional is determined at compile time. */ \ |
822 | if (RETURN_LIMB_SIZE > 1) \ | |
823 | retval[1] = 0; \ | |
28f540f4 | 824 | retval[0] = quot; \ |
4933a099 | 825 | bits = -cnt; \ |
28f540f4 RM |
826 | } \ |
827 | } \ | |
828 | else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \ | |
4933a099 | 829 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \ |
28f540f4 RM |
830 | quot); \ |
831 | else \ | |
832 | { \ | |
833 | used = MANT_DIG - bits; \ | |
834 | if (used > 0) \ | |
4933a099 | 835 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \ |
28f540f4 RM |
836 | } \ |
837 | bits += BITS_PER_MP_LIMB | |
838 | ||
4933a099 | 839 | got_limb; |
28f540f4 RM |
840 | } |
841 | while (bits <= MANT_DIG); | |
842 | ||
843 | return round_and_return (retval, exponent - 1, negative, | |
844 | quot, BITS_PER_MP_LIMB - 1 - used, | |
4933a099 | 845 | more_bits || n != 0); |
28f540f4 RM |
846 | } |
847 | case 2: | |
848 | { | |
849 | mp_limb d0, d1, n0, n1; | |
850 | mp_limb quot = 0; | |
851 | int used = 0; | |
852 | ||
853 | d0 = den[0]; | |
854 | d1 = den[1]; | |
855 | ||
856 | if (numsize < densize) | |
857 | { | |
858 | if (bits <= 0) | |
859 | exponent -= BITS_PER_MP_LIMB; | |
860 | else | |
861 | { | |
862 | if (bits + BITS_PER_MP_LIMB <= MANT_DIG) | |
863 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, | |
864 | BITS_PER_MP_LIMB, 0); | |
865 | else | |
866 | { | |
867 | used = MANT_DIG - bits; | |
868 | if (used > 0) | |
869 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); | |
870 | } | |
871 | bits += BITS_PER_MP_LIMB; | |
872 | } | |
873 | n1 = num[0]; | |
874 | n0 = 0; | |
875 | } | |
876 | else | |
877 | { | |
878 | n1 = num[1]; | |
879 | n0 = num[0]; | |
880 | } | |
881 | ||
882 | while (bits <= MANT_DIG) | |
883 | { | |
884 | mp_limb r; | |
885 | ||
886 | if (n1 == d1) | |
887 | { | |
888 | /* QUOT should be either 111..111 or 111..110. We need | |
889 | special treatment of this rare case as normal division | |
890 | would give overflow. */ | |
891 | quot = ~(mp_limb) 0; | |
892 | ||
893 | r = n0 + d1; | |
894 | if (r < d1) /* Carry in the addition? */ | |
895 | { | |
896 | add_ssaaaa (n1, n0, r - d0, 0, 0, d0); | |
897 | goto have_quot; | |
898 | } | |
899 | n1 = d0 - (d0 != 0); | |
900 | n0 = -d0; | |
901 | } | |
902 | else | |
903 | { | |
904 | udiv_qrnnd (quot, r, n1, n0, d1); | |
905 | umul_ppmm (n1, n0, d0, quot); | |
906 | } | |
907 | ||
908 | q_test: | |
909 | if (n1 > r || (n1 == r && n0 > 0)) | |
910 | { | |
911 | /* The estimated QUOT was too large. */ | |
912 | --quot; | |
913 | ||
914 | sub_ddmmss (n1, n0, n1, n0, 0, d0); | |
915 | r += d1; | |
916 | if (r >= d1) /* If not carry, test QUOT again. */ | |
917 | goto q_test; | |
918 | } | |
919 | sub_ddmmss (n1, n0, r, 0, n1, n0); | |
920 | ||
921 | have_quot: | |
922 | got_limb; | |
923 | } | |
924 | ||
925 | return round_and_return (retval, exponent - 1, negative, | |
926 | quot, BITS_PER_MP_LIMB - 1 - used, | |
4933a099 | 927 | more_bits || n1 != 0 || n0 != 0); |
28f540f4 RM |
928 | } |
929 | default: | |
930 | { | |
931 | int i; | |
932 | mp_limb cy, dX, d1, n0, n1; | |
933 | mp_limb quot = 0; | |
934 | int used = 0; | |
935 | ||
936 | dX = den[densize - 1]; | |
937 | d1 = den[densize - 2]; | |
938 | ||
939 | /* The division does not work if the upper limb of the two-limb | |
940 | numerator is greater than the denominator. */ | |
4933a099 | 941 | if (__mpn_cmp (num, &den[densize - numsize], numsize) > 0) |
28f540f4 RM |
942 | num[numsize++] = 0; |
943 | ||
944 | if (numsize < densize) | |
945 | { | |
946 | mp_size_t empty = densize - numsize; | |
947 | ||
948 | if (bits <= 0) | |
949 | { | |
950 | register int i; | |
951 | for (i = numsize; i > 0; --i) | |
952 | num[i + empty] = num[i - 1]; | |
953 | MPN_ZERO (num, empty + 1); | |
954 | exponent -= empty * BITS_PER_MP_LIMB; | |
955 | } | |
956 | else | |
957 | { | |
958 | if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG) | |
959 | { | |
960 | /* We make a difference here because the compiler | |
961 | cannot optimize the `else' case that good and | |
962 | this reflects all currently used FLOAT types | |
963 | and GMP implementations. */ | |
964 | register int i; | |
965 | #if RETURN_LIMB_SIZE <= 2 | |
966 | assert (empty == 1); | |
967 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, | |
968 | BITS_PER_MP_LIMB, 0); | |
969 | #else | |
970 | for (i = RETURN_LIMB_SIZE; i > empty; --i) | |
971 | retval[i] = retval[i - empty]; | |
972 | #endif | |
973 | retval[1] = 0; | |
974 | for (i = numsize; i > 0; --i) | |
975 | num[i + empty] = num[i - 1]; | |
976 | MPN_ZERO (num, empty + 1); | |
977 | } | |
978 | else | |
979 | { | |
980 | used = MANT_DIG - bits; | |
981 | if (used >= BITS_PER_MP_LIMB) | |
982 | { | |
983 | register int i; | |
984 | (void) __mpn_lshift (&retval[used | |
985 | / BITS_PER_MP_LIMB], | |
986 | retval, RETURN_LIMB_SIZE, | |
987 | used % BITS_PER_MP_LIMB); | |
988 | for (i = used / BITS_PER_MP_LIMB; i >= 0; --i) | |
989 | retval[i] = 0; | |
990 | } | |
991 | else if (used > 0) | |
992 | __mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); | |
993 | } | |
994 | bits += empty * BITS_PER_MP_LIMB; | |
995 | } | |
996 | } | |
997 | else | |
998 | { | |
999 | int i; | |
1000 | assert (numsize == densize); | |
1001 | for (i = numsize; i > 0; --i) | |
1002 | num[i] = num[i - 1]; | |
1003 | } | |
1004 | ||
1005 | den[densize] = 0; | |
1006 | n0 = num[densize]; | |
1007 | ||
1008 | while (bits <= MANT_DIG) | |
1009 | { | |
1010 | if (n0 == dX) | |
1011 | /* This might over-estimate QUOT, but it's probably not | |
1012 | worth the extra code here to find out. */ | |
1013 | quot = ~(mp_limb) 0; | |
1014 | else | |
1015 | { | |
1016 | mp_limb r; | |
1017 | ||
1018 | udiv_qrnnd (quot, r, n0, num[densize - 1], dX); | |
1019 | umul_ppmm (n1, n0, d1, quot); | |
1020 | ||
1021 | while (n1 > r || (n1 == r && n0 > num[densize - 2])) | |
1022 | { | |
1023 | --quot; | |
1024 | r += dX; | |
1025 | if (r < dX) /* I.e. "carry in previous addition?" */ | |
1026 | break; | |
1027 | n1 -= n0 < d1; | |
1028 | n0 -= d1; | |
1029 | } | |
1030 | } | |
1031 | ||
1032 | /* Possible optimization: We already have (q * n0) and (1 * n1) | |
1033 | after the calculation of QUOT. Taking advantage of this, we | |
1034 | could make this loop make two iterations less. */ | |
1035 | ||
1036 | cy = __mpn_submul_1 (num, den, densize + 1, quot); | |
1037 | ||
1038 | if (num[densize] != cy) | |
1039 | { | |
1040 | cy = __mpn_add_n (num, num, den, densize); | |
1041 | assert (cy != 0); | |
1042 | --quot; | |
1043 | } | |
1044 | n0 = num[densize] = num[densize - 1]; | |
1045 | for (i = densize - 1; i > 0; --i) | |
1046 | num[i] = num[i - 1]; | |
1047 | ||
1048 | got_limb; | |
1049 | } | |
1050 | ||
4933a099 | 1051 | for (i = densize; num[i] == 0 && i >= 0; --i) |
28f540f4 RM |
1052 | ; |
1053 | return round_and_return (retval, exponent - 1, negative, | |
1054 | quot, BITS_PER_MP_LIMB - 1 - used, | |
4933a099 | 1055 | more_bits || i >= 0); |
28f540f4 RM |
1056 | } |
1057 | } | |
1058 | } | |
1059 | ||
1060 | /* NOTREACHED */ | |
1061 | } |