1 | /* |
---|
2 | FUNCTION |
---|
3 | <<strtod>>, <<strtof>>, <<strtold>>, <<strtod_l>>, <<strtof_l>>, <<strtold_l>>---string to double or float |
---|
4 | |
---|
5 | INDEX |
---|
6 | strtod |
---|
7 | |
---|
8 | INDEX |
---|
9 | strtof |
---|
10 | |
---|
11 | INDEX |
---|
12 | strtold |
---|
13 | |
---|
14 | INDEX |
---|
15 | strtod_l |
---|
16 | |
---|
17 | INDEX |
---|
18 | strtof_l |
---|
19 | |
---|
20 | INDEX |
---|
21 | strtold_l |
---|
22 | |
---|
23 | INDEX |
---|
24 | _strtod_r |
---|
25 | |
---|
26 | SYNOPSIS |
---|
27 | #include <stdlib.h> |
---|
28 | double strtod(const char *restrict <[str]>, char **restrict <[tail]>); |
---|
29 | float strtof(const char *restrict <[str]>, char **restrict <[tail]>); |
---|
30 | long double strtold(const char *restrict <[str]>, |
---|
31 | char **restrict <[tail]>); |
---|
32 | |
---|
33 | #include <stdlib.h> |
---|
34 | double strtod_l(const char *restrict <[str]>, char **restrict <[tail]>, |
---|
35 | locale_t <[locale]>); |
---|
36 | float strtof_l(const char *restrict <[str]>, char **restrict <[tail]>, |
---|
37 | locale_t <[locale]>); |
---|
38 | long double strtold_l(const char *restrict <[str]>, |
---|
39 | char **restrict <[tail]>, |
---|
40 | locale_t <[locale]>); |
---|
41 | |
---|
42 | double _strtod_r(void *<[reent]>, |
---|
43 | const char *restrict <[str]>, char **restrict <[tail]>); |
---|
44 | |
---|
45 | DESCRIPTION |
---|
46 | <<strtod>>, <<strtof>>, <<strtold>> parse the character string |
---|
47 | <[str]>, producing a substring which can be converted to a double, |
---|
48 | float, or long double value, respectively. The substring converted |
---|
49 | is the longest initial subsequence of <[str]>, beginning with the |
---|
50 | first non-whitespace character, that has one of these formats: |
---|
51 | .[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>] |
---|
52 | .[+|-].<[digits]>[(e|E)[+|-]<[digits]>] |
---|
53 | .[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)] |
---|
54 | .[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>] |
---|
55 | .[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>] |
---|
56 | .[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[digits]>] |
---|
57 | The substring contains no characters if <[str]> is empty, consists |
---|
58 | entirely of whitespace, or if the first non-whitespace |
---|
59 | character is something other than <<+>>, <<->>, <<.>>, or a |
---|
60 | digit, and cannot be parsed as infinity or NaN. If the platform |
---|
61 | does not support NaN, then NaN is treated as an empty substring. |
---|
62 | If the substring is empty, no conversion is done, and |
---|
63 | the value of <[str]> is stored in <<*<[tail]>>>. Otherwise, |
---|
64 | the substring is converted, and a pointer to the final string |
---|
65 | (which will contain at least the terminating null character of |
---|
66 | <[str]>) is stored in <<*<[tail]>>>. If you want no |
---|
67 | assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>. |
---|
68 | |
---|
69 | This implementation returns the nearest machine number to the |
---|
70 | input decimal string. Ties are broken by using the IEEE |
---|
71 | round-even rule. However, <<strtof>> is currently subject to |
---|
72 | double rounding errors. |
---|
73 | |
---|
74 | <<strtod_l>>, <<strtof_l>>, <<strtold_l>> are like <<strtod>>, |
---|
75 | <<strtof>>, <<strtold>> but perform the conversion based on the |
---|
76 | locale specified by the locale object locale. If <[locale]> is |
---|
77 | LC_GLOBAL_LOCALE or not a valid locale object, the behaviour is |
---|
78 | undefined. |
---|
79 | |
---|
80 | The alternate function <<_strtod_r>> is a reentrant version. |
---|
81 | The extra argument <[reent]> is a pointer to a reentrancy structure. |
---|
82 | |
---|
83 | RETURNS |
---|
84 | These functions return the converted substring value, if any. If |
---|
85 | no conversion could be performed, 0 is returned. If the correct |
---|
86 | value is out of the range of representable values, plus or minus |
---|
87 | <<HUGE_VAL>> (<<HUGE_VALF>>, <<HUGE_VALL>>) is returned, and |
---|
88 | <<ERANGE>> is stored in errno. If the correct value would cause |
---|
89 | underflow, 0 is returned and <<ERANGE>> is stored in errno. |
---|
90 | |
---|
91 | PORTABILITY |
---|
92 | <<strtod>> is ANSI. |
---|
93 | <<strtof>>, <<strtold>> are C99. |
---|
94 | <<strtod_l>>, <<strtof_l>>, <<strtold_l>> are GNU extensions. |
---|
95 | |
---|
96 | Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>, |
---|
97 | <<lseek>>, <<read>>, <<sbrk>>, <<write>>. |
---|
98 | */ |
---|
99 | |
---|
100 | /**************************************************************** |
---|
101 | |
---|
102 | The author of this software is David M. Gay. |
---|
103 | |
---|
104 | Copyright (C) 1998-2001 by Lucent Technologies |
---|
105 | All Rights Reserved |
---|
106 | |
---|
107 | Permission to use, copy, modify, and distribute this software and |
---|
108 | its documentation for any purpose and without fee is hereby |
---|
109 | granted, provided that the above copyright notice appear in all |
---|
110 | copies and that both that the copyright notice and this |
---|
111 | permission notice and warranty disclaimer appear in supporting |
---|
112 | documentation, and that the name of Lucent or any of its entities |
---|
113 | not be used in advertising or publicity pertaining to |
---|
114 | distribution of the software without specific, written prior |
---|
115 | permission. |
---|
116 | |
---|
117 | LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, |
---|
118 | INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. |
---|
119 | IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY |
---|
120 | SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
---|
121 | WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER |
---|
122 | IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, |
---|
123 | ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF |
---|
124 | THIS SOFTWARE. |
---|
125 | |
---|
126 | ****************************************************************/ |
---|
127 | |
---|
128 | /* Please send bug reports to David M. Gay (dmg at acm dot org, |
---|
129 | * with " at " changed at "@" and " dot " changed to "."). */ |
---|
130 | |
---|
131 | /* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib. */ |
---|
132 | |
---|
133 | #define _GNU_SOURCE |
---|
134 | #include <_ansi.h> |
---|
135 | #include <errno.h> |
---|
136 | #include <stdlib.h> |
---|
137 | #include <string.h> |
---|
138 | #include "mprec.h" |
---|
139 | #include "gdtoa.h" |
---|
140 | #include "gd_qnan.h" |
---|
141 | #include "../locale/setlocale.h" |
---|
142 | |
---|
143 | /* #ifndef NO_FENV_H */ |
---|
144 | /* #include <fenv.h> */ |
---|
145 | /* #endif */ |
---|
146 | |
---|
147 | #include "locale.h" |
---|
148 | |
---|
149 | #ifdef IEEE_Arith |
---|
150 | #ifndef NO_IEEE_Scale |
---|
151 | #define Avoid_Underflow |
---|
152 | #undef tinytens |
---|
153 | /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */ |
---|
154 | /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ |
---|
155 | static const double tinytens[] = { 1e-16, 1e-32, |
---|
156 | #ifdef _DOUBLE_IS_32BITS |
---|
157 | 0.0, 0.0, 0.0 |
---|
158 | #else |
---|
159 | 1e-64, 1e-128, |
---|
160 | 9007199254740992. * 9007199254740992.e-256 |
---|
161 | #endif |
---|
162 | }; |
---|
163 | |
---|
164 | #endif |
---|
165 | #endif |
---|
166 | |
---|
167 | #ifdef Honor_FLT_ROUNDS |
---|
168 | #define Rounding rounding |
---|
169 | #undef Check_FLT_ROUNDS |
---|
170 | #define Check_FLT_ROUNDS |
---|
171 | #else |
---|
172 | #define Rounding Flt_Rounds |
---|
173 | #endif |
---|
174 | |
---|
175 | #ifdef Avoid_Underflow /*{*/ |
---|
176 | static double |
---|
177 | sulp (U x, |
---|
178 | int scale) |
---|
179 | { |
---|
180 | U u; |
---|
181 | double rv; |
---|
182 | int i; |
---|
183 | |
---|
184 | rv = ulp(dval(x)); |
---|
185 | if (!scale || (i = 2*P + 1 - ((dword0(x) & Exp_mask) >> Exp_shift)) <= 0) |
---|
186 | return rv; /* Is there an example where i <= 0 ? */ |
---|
187 | dword0(u) = Exp_1 + ((__int32_t)i << Exp_shift); |
---|
188 | #ifndef _DOUBLE_IS_32BITS |
---|
189 | dword1(u) = 0; |
---|
190 | #endif |
---|
191 | return rv * u.d; |
---|
192 | } |
---|
193 | #endif /*}*/ |
---|
194 | |
---|
195 | |
---|
196 | #ifndef NO_HEX_FP |
---|
197 | |
---|
198 | static void |
---|
199 | ULtod (__ULong *L, |
---|
200 | __ULong *bits, |
---|
201 | Long exp, |
---|
202 | int k) |
---|
203 | { |
---|
204 | switch(k & STRTOG_Retmask) { |
---|
205 | case STRTOG_NoNumber: |
---|
206 | case STRTOG_Zero: |
---|
207 | L[0] = L[1] = 0; |
---|
208 | break; |
---|
209 | |
---|
210 | case STRTOG_Denormal: |
---|
211 | L[_1] = bits[0]; |
---|
212 | L[_0] = bits[1]; |
---|
213 | break; |
---|
214 | |
---|
215 | case STRTOG_Normal: |
---|
216 | case STRTOG_NaNbits: |
---|
217 | L[_1] = bits[0]; |
---|
218 | L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20); |
---|
219 | break; |
---|
220 | |
---|
221 | case STRTOG_Infinite: |
---|
222 | L[_0] = 0x7ff00000; |
---|
223 | L[_1] = 0; |
---|
224 | break; |
---|
225 | |
---|
226 | case STRTOG_NaN: |
---|
227 | L[_0] = 0x7fffffff; |
---|
228 | L[_1] = (__ULong)-1; |
---|
229 | } |
---|
230 | if (k & STRTOG_Neg) |
---|
231 | L[_0] |= 0x80000000L; |
---|
232 | } |
---|
233 | #endif /* !NO_HEX_FP */ |
---|
234 | |
---|
235 | double |
---|
236 | _strtod_l (struct _reent *ptr, const char *__restrict s00, char **__restrict se, |
---|
237 | locale_t loc) |
---|
238 | { |
---|
239 | #ifdef Avoid_Underflow |
---|
240 | int scale; |
---|
241 | #endif |
---|
242 | int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, |
---|
243 | e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; |
---|
244 | const char *s, *s0, *s1; |
---|
245 | double aadj, adj; |
---|
246 | U aadj1, rv, rv0; |
---|
247 | Long L; |
---|
248 | __ULong y, z; |
---|
249 | _Bigint *bb = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL; |
---|
250 | #ifdef Avoid_Underflow |
---|
251 | __ULong Lsb, Lsb1; |
---|
252 | #endif |
---|
253 | #ifdef SET_INEXACT |
---|
254 | int inexact, oldinexact; |
---|
255 | #endif |
---|
256 | #ifdef Honor_FLT_ROUNDS |
---|
257 | int rounding; |
---|
258 | #endif |
---|
259 | struct lconv *lconv = __localeconv_l (loc); |
---|
260 | int dec_len = strlen (lconv->decimal_point); |
---|
261 | |
---|
262 | delta = bs = bd = NULL; |
---|
263 | sign = nz0 = nz = decpt = 0; |
---|
264 | dval(rv) = 0.; |
---|
265 | for(s = s00;;s++) switch(*s) { |
---|
266 | case '-': |
---|
267 | sign = 1; |
---|
268 | /* no break */ |
---|
269 | case '+': |
---|
270 | if (*++s) |
---|
271 | goto break2; |
---|
272 | /* no break */ |
---|
273 | case 0: |
---|
274 | goto ret0; |
---|
275 | case '\t': |
---|
276 | case '\n': |
---|
277 | case '\v': |
---|
278 | case '\f': |
---|
279 | case '\r': |
---|
280 | case ' ': |
---|
281 | continue; |
---|
282 | default: |
---|
283 | goto break2; |
---|
284 | } |
---|
285 | break2: |
---|
286 | if (*s == '0') { |
---|
287 | #ifndef NO_HEX_FP |
---|
288 | { |
---|
289 | static const FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; |
---|
290 | Long exp; |
---|
291 | __ULong bits[2]; |
---|
292 | switch(s[1]) { |
---|
293 | case 'x': |
---|
294 | case 'X': |
---|
295 | /* If the number is not hex, then the parse of |
---|
296 | 0 is still valid. */ |
---|
297 | s00 = s + 1; |
---|
298 | { |
---|
299 | #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) |
---|
300 | FPI fpi1 = fpi; |
---|
301 | switch(fegetround()) { |
---|
302 | case FE_TOWARDZERO: fpi1.rounding = 0; break; |
---|
303 | case FE_UPWARD: fpi1.rounding = 2; break; |
---|
304 | case FE_DOWNWARD: fpi1.rounding = 3; |
---|
305 | } |
---|
306 | #else |
---|
307 | #define fpi1 fpi |
---|
308 | #endif |
---|
309 | switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign, loc)) & STRTOG_Retmask) { |
---|
310 | case STRTOG_NoNumber: |
---|
311 | s = s00; |
---|
312 | sign = 0; |
---|
313 | /* FALLTHROUGH */ |
---|
314 | case STRTOG_Zero: |
---|
315 | break; |
---|
316 | default: |
---|
317 | if (bb) { |
---|
318 | copybits(bits, fpi.nbits, bb); |
---|
319 | Bfree(ptr,bb); |
---|
320 | } |
---|
321 | ULtod(rv.i, bits, exp, i); |
---|
322 | }} |
---|
323 | goto ret; |
---|
324 | } |
---|
325 | } |
---|
326 | #endif |
---|
327 | nz0 = 1; |
---|
328 | while(*++s == '0') ; |
---|
329 | if (!*s) |
---|
330 | goto ret; |
---|
331 | } |
---|
332 | s0 = s; |
---|
333 | y = z = 0; |
---|
334 | for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) |
---|
335 | if (nd < 9) |
---|
336 | y = 10*y + c - '0'; |
---|
337 | else |
---|
338 | z = 10*z + c - '0'; |
---|
339 | nd0 = nd; |
---|
340 | if (strncmp (s, lconv->decimal_point, dec_len) == 0) |
---|
341 | { |
---|
342 | decpt = 1; |
---|
343 | c = *(s += dec_len); |
---|
344 | if (!nd) { |
---|
345 | for(; c == '0'; c = *++s) |
---|
346 | nz++; |
---|
347 | if (c > '0' && c <= '9') { |
---|
348 | s0 = s; |
---|
349 | nf += nz; |
---|
350 | nz = 0; |
---|
351 | goto have_dig; |
---|
352 | } |
---|
353 | goto dig_done; |
---|
354 | } |
---|
355 | for(; c >= '0' && c <= '9'; c = *++s) { |
---|
356 | have_dig: |
---|
357 | nz++; |
---|
358 | if (c -= '0') { |
---|
359 | nf += nz; |
---|
360 | for(i = 1; i < nz; i++) |
---|
361 | if (nd++ < 9) |
---|
362 | y *= 10; |
---|
363 | else if (nd <= DBL_DIG + 1) |
---|
364 | z *= 10; |
---|
365 | if (nd++ < 9) |
---|
366 | y = 10*y + c; |
---|
367 | else if (nd <= DBL_DIG + 1) |
---|
368 | z = 10*z + c; |
---|
369 | nz = 0; |
---|
370 | } |
---|
371 | } |
---|
372 | } |
---|
373 | dig_done: |
---|
374 | e = 0; |
---|
375 | if (c == 'e' || c == 'E') { |
---|
376 | if (!nd && !nz && !nz0) { |
---|
377 | goto ret0; |
---|
378 | } |
---|
379 | s00 = s; |
---|
380 | esign = 0; |
---|
381 | switch(c = *++s) { |
---|
382 | case '-': |
---|
383 | esign = 1; |
---|
384 | case '+': |
---|
385 | c = *++s; |
---|
386 | } |
---|
387 | if (c >= '0' && c <= '9') { |
---|
388 | while(c == '0') |
---|
389 | c = *++s; |
---|
390 | if (c > '0' && c <= '9') { |
---|
391 | L = c - '0'; |
---|
392 | s1 = s; |
---|
393 | while((c = *++s) >= '0' && c <= '9') |
---|
394 | L = 10*L + c - '0'; |
---|
395 | if (s - s1 > 8 || L > 19999) |
---|
396 | /* Avoid confusion from exponents |
---|
397 | * so large that e might overflow. |
---|
398 | */ |
---|
399 | e = 19999; /* safe for 16 bit ints */ |
---|
400 | else |
---|
401 | e = (int)L; |
---|
402 | if (esign) |
---|
403 | e = -e; |
---|
404 | } |
---|
405 | else |
---|
406 | e = 0; |
---|
407 | } |
---|
408 | else |
---|
409 | s = s00; |
---|
410 | } |
---|
411 | if (!nd) { |
---|
412 | if (!nz && !nz0) { |
---|
413 | #ifdef INFNAN_CHECK |
---|
414 | /* Check for Nan and Infinity */ |
---|
415 | __ULong bits[2]; |
---|
416 | static const FPI fpinan = /* only 52 explicit bits */ |
---|
417 | { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; |
---|
418 | if (!decpt) |
---|
419 | switch(c) { |
---|
420 | case 'i': |
---|
421 | case 'I': |
---|
422 | if (match(&s,"nf")) { |
---|
423 | --s; |
---|
424 | if (!match(&s,"inity")) |
---|
425 | ++s; |
---|
426 | dword0(rv) = 0x7ff00000; |
---|
427 | #ifndef _DOUBLE_IS_32BITS |
---|
428 | dword1(rv) = 0; |
---|
429 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
430 | goto ret; |
---|
431 | } |
---|
432 | break; |
---|
433 | case 'n': |
---|
434 | case 'N': |
---|
435 | if (match(&s, "an")) { |
---|
436 | #ifndef No_Hex_NaN |
---|
437 | if (*s == '(' /*)*/ |
---|
438 | && hexnan(&s, &fpinan, bits) |
---|
439 | == STRTOG_NaNbits) { |
---|
440 | dword0(rv) = 0x7ff00000 | bits[1]; |
---|
441 | #ifndef _DOUBLE_IS_32BITS |
---|
442 | dword1(rv) = bits[0]; |
---|
443 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
444 | } |
---|
445 | else { |
---|
446 | #endif |
---|
447 | dword0(rv) = NAN_WORD0; |
---|
448 | #ifndef _DOUBLE_IS_32BITS |
---|
449 | dword1(rv) = NAN_WORD1; |
---|
450 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
451 | #ifndef No_Hex_NaN |
---|
452 | } |
---|
453 | #endif |
---|
454 | goto ret; |
---|
455 | } |
---|
456 | } |
---|
457 | #endif /* INFNAN_CHECK */ |
---|
458 | ret0: |
---|
459 | s = s00; |
---|
460 | sign = 0; |
---|
461 | } |
---|
462 | goto ret; |
---|
463 | } |
---|
464 | e1 = e -= nf; |
---|
465 | |
---|
466 | /* Now we have nd0 digits, starting at s0, followed by a |
---|
467 | * decimal point, followed by nd-nd0 digits. The number we're |
---|
468 | * after is the integer represented by those digits times |
---|
469 | * 10**e */ |
---|
470 | |
---|
471 | if (!nd0) |
---|
472 | nd0 = nd; |
---|
473 | k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; |
---|
474 | dval(rv) = y; |
---|
475 | if (k > 9) { |
---|
476 | #ifdef SET_INEXACT |
---|
477 | if (k > DBL_DIG) |
---|
478 | oldinexact = get_inexact(); |
---|
479 | #endif |
---|
480 | dval(rv) = tens[k - 9] * dval(rv) + z; |
---|
481 | } |
---|
482 | bd0 = 0; |
---|
483 | if (nd <= DBL_DIG |
---|
484 | #ifndef RND_PRODQUOT |
---|
485 | #ifndef Honor_FLT_ROUNDS |
---|
486 | && Flt_Rounds == 1 |
---|
487 | #endif |
---|
488 | #endif |
---|
489 | ) { |
---|
490 | if (!e) |
---|
491 | goto ret; |
---|
492 | if (e > 0) { |
---|
493 | if (e <= Ten_pmax) { |
---|
494 | #ifdef VAX |
---|
495 | goto vax_ovfl_check; |
---|
496 | #else |
---|
497 | #ifdef Honor_FLT_ROUNDS |
---|
498 | /* round correctly FLT_ROUNDS = 2 or 3 */ |
---|
499 | if (sign) { |
---|
500 | dval(rv) = -dval(rv); |
---|
501 | sign = 0; |
---|
502 | } |
---|
503 | #endif |
---|
504 | /* rv = */ rounded_product(dval(rv), tens[e]); |
---|
505 | goto ret; |
---|
506 | #endif |
---|
507 | } |
---|
508 | i = DBL_DIG - nd; |
---|
509 | if (e <= Ten_pmax + i) { |
---|
510 | /* A fancier test would sometimes let us do |
---|
511 | * this for larger i values. |
---|
512 | */ |
---|
513 | #ifdef Honor_FLT_ROUNDS |
---|
514 | /* round correctly FLT_ROUNDS = 2 or 3 */ |
---|
515 | if (sign) { |
---|
516 | dval(rv) = -dval(rv); |
---|
517 | sign = 0; |
---|
518 | } |
---|
519 | #endif |
---|
520 | e -= i; |
---|
521 | dval(rv) *= tens[i]; |
---|
522 | #ifdef VAX |
---|
523 | /* VAX exponent range is so narrow we must |
---|
524 | * worry about overflow here... |
---|
525 | */ |
---|
526 | vax_ovfl_check: |
---|
527 | dword0(rv) -= P*Exp_msk1; |
---|
528 | /* rv = */ rounded_product(dval(rv), tens[e]); |
---|
529 | if ((dword0(rv) & Exp_mask) |
---|
530 | > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) |
---|
531 | goto ovfl; |
---|
532 | dword0(rv) += P*Exp_msk1; |
---|
533 | #else |
---|
534 | /* rv = */ rounded_product(dval(rv), tens[e]); |
---|
535 | #endif |
---|
536 | goto ret; |
---|
537 | } |
---|
538 | } |
---|
539 | #ifndef Inaccurate_Divide |
---|
540 | else if (e >= -Ten_pmax) { |
---|
541 | #ifdef Honor_FLT_ROUNDS |
---|
542 | /* round correctly FLT_ROUNDS = 2 or 3 */ |
---|
543 | if (sign) { |
---|
544 | dval(rv) = -dval(rv); |
---|
545 | sign = 0; |
---|
546 | } |
---|
547 | #endif |
---|
548 | /* rv = */ rounded_quotient(dval(rv), tens[-e]); |
---|
549 | goto ret; |
---|
550 | } |
---|
551 | #endif |
---|
552 | } |
---|
553 | e1 += nd - k; |
---|
554 | |
---|
555 | #ifdef IEEE_Arith |
---|
556 | #ifdef SET_INEXACT |
---|
557 | inexact = 1; |
---|
558 | if (k <= DBL_DIG) |
---|
559 | oldinexact = get_inexact(); |
---|
560 | #endif |
---|
561 | #ifdef Avoid_Underflow |
---|
562 | scale = 0; |
---|
563 | #endif |
---|
564 | #ifdef Honor_FLT_ROUNDS |
---|
565 | if ((rounding = Flt_Rounds) >= 2) { |
---|
566 | if (sign) |
---|
567 | rounding = rounding == 2 ? 0 : 2; |
---|
568 | else |
---|
569 | if (rounding != 2) |
---|
570 | rounding = 0; |
---|
571 | } |
---|
572 | #endif |
---|
573 | #endif /*IEEE_Arith*/ |
---|
574 | |
---|
575 | /* Get starting approximation = rv * 10**e1 */ |
---|
576 | |
---|
577 | if (e1 > 0) { |
---|
578 | if ( (i = e1 & 15) !=0) |
---|
579 | dval(rv) *= tens[i]; |
---|
580 | if (e1 &= ~15) { |
---|
581 | if (e1 > DBL_MAX_10_EXP) { |
---|
582 | ovfl: |
---|
583 | #ifndef NO_ERRNO |
---|
584 | ptr->_errno = ERANGE; |
---|
585 | #endif |
---|
586 | /* Can't trust HUGE_VAL */ |
---|
587 | #ifdef IEEE_Arith |
---|
588 | #ifdef Honor_FLT_ROUNDS |
---|
589 | switch(rounding) { |
---|
590 | case 0: /* toward 0 */ |
---|
591 | case 3: /* toward -infinity */ |
---|
592 | dword0(rv) = Big0; |
---|
593 | #ifndef _DOUBLE_IS_32BITS |
---|
594 | dword1(rv) = Big1; |
---|
595 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
596 | break; |
---|
597 | default: |
---|
598 | dword0(rv) = Exp_mask; |
---|
599 | #ifndef _DOUBLE_IS_32BITS |
---|
600 | dword1(rv) = 0; |
---|
601 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
602 | } |
---|
603 | #else /*Honor_FLT_ROUNDS*/ |
---|
604 | dword0(rv) = Exp_mask; |
---|
605 | #ifndef _DOUBLE_IS_32BITS |
---|
606 | dword1(rv) = 0; |
---|
607 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
608 | #endif /*Honor_FLT_ROUNDS*/ |
---|
609 | #ifdef SET_INEXACT |
---|
610 | /* set overflow bit */ |
---|
611 | dval(rv0) = 1e300; |
---|
612 | dval(rv0) *= dval(rv0); |
---|
613 | #endif |
---|
614 | #else /*IEEE_Arith*/ |
---|
615 | dword0(rv) = Big0; |
---|
616 | #ifndef _DOUBLE_IS_32BITS |
---|
617 | dword1(rv) = Big1; |
---|
618 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
619 | #endif /*IEEE_Arith*/ |
---|
620 | if (bd0) |
---|
621 | goto retfree; |
---|
622 | goto ret; |
---|
623 | } |
---|
624 | e1 >>= 4; |
---|
625 | for(j = 0; e1 > 1; j++, e1 >>= 1) |
---|
626 | if (e1 & 1) |
---|
627 | dval(rv) *= bigtens[j]; |
---|
628 | /* The last multiplication could overflow. */ |
---|
629 | dword0(rv) -= P*Exp_msk1; |
---|
630 | dval(rv) *= bigtens[j]; |
---|
631 | if ((z = dword0(rv) & Exp_mask) |
---|
632 | > Exp_msk1*(DBL_MAX_EXP+Bias-P)) |
---|
633 | goto ovfl; |
---|
634 | if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { |
---|
635 | /* set to largest number */ |
---|
636 | /* (Can't trust DBL_MAX) */ |
---|
637 | dword0(rv) = Big0; |
---|
638 | #ifndef _DOUBLE_IS_32BITS |
---|
639 | dword1(rv) = Big1; |
---|
640 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
641 | } |
---|
642 | else |
---|
643 | dword0(rv) += P*Exp_msk1; |
---|
644 | } |
---|
645 | } |
---|
646 | else if (e1 < 0) { |
---|
647 | e1 = -e1; |
---|
648 | if ( (i = e1 & 15) !=0) |
---|
649 | dval(rv) /= tens[i]; |
---|
650 | if (e1 >>= 4) { |
---|
651 | if (e1 >= 1 << n_bigtens) |
---|
652 | goto undfl; |
---|
653 | #ifdef Avoid_Underflow |
---|
654 | if (e1 & Scale_Bit) |
---|
655 | scale = 2*P; |
---|
656 | for(j = 0; e1 > 0; j++, e1 >>= 1) |
---|
657 | if (e1 & 1) |
---|
658 | dval(rv) *= tinytens[j]; |
---|
659 | if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask) |
---|
660 | >> Exp_shift)) > 0) { |
---|
661 | /* scaled rv is denormal; zap j low bits */ |
---|
662 | if (j >= 32) { |
---|
663 | #ifndef _DOUBLE_IS_32BITS |
---|
664 | dword1(rv) = 0; |
---|
665 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
666 | if (j >= 53) |
---|
667 | dword0(rv) = (P+2)*Exp_msk1; |
---|
668 | else |
---|
669 | dword0(rv) &= 0xffffffff << (j-32); |
---|
670 | } |
---|
671 | #ifndef _DOUBLE_IS_32BITS |
---|
672 | else |
---|
673 | dword1(rv) &= 0xffffffff << j; |
---|
674 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
675 | } |
---|
676 | #else |
---|
677 | for(j = 0; e1 > 1; j++, e1 >>= 1) |
---|
678 | if (e1 & 1) |
---|
679 | dval(rv) *= tinytens[j]; |
---|
680 | /* The last multiplication could underflow. */ |
---|
681 | dval(rv0) = dval(rv); |
---|
682 | dval(rv) *= tinytens[j]; |
---|
683 | if (!dval(rv)) { |
---|
684 | dval(rv) = 2.*dval(rv0); |
---|
685 | dval(rv) *= tinytens[j]; |
---|
686 | #endif |
---|
687 | if (!dval(rv)) { |
---|
688 | undfl: |
---|
689 | dval(rv) = 0.; |
---|
690 | #ifndef NO_ERRNO |
---|
691 | ptr->_errno = ERANGE; |
---|
692 | #endif |
---|
693 | if (bd0) |
---|
694 | goto retfree; |
---|
695 | goto ret; |
---|
696 | } |
---|
697 | #ifndef Avoid_Underflow |
---|
698 | #ifndef _DOUBLE_IS_32BITS |
---|
699 | dword0(rv) = Tiny0; |
---|
700 | dword1(rv) = Tiny1; |
---|
701 | #else |
---|
702 | dword0(rv) = Tiny1; |
---|
703 | #endif /*_DOUBLE_IS_32BITS*/ |
---|
704 | /* The refinement below will clean |
---|
705 | * this approximation up. |
---|
706 | */ |
---|
707 | } |
---|
708 | #endif |
---|
709 | } |
---|
710 | } |
---|
711 | |
---|
712 | /* Now the hard part -- adjusting rv to the correct value.*/ |
---|
713 | |
---|
714 | /* Put digits into bd: true value = bd * 10^e */ |
---|
715 | |
---|
716 | bd0 = s2b(ptr, s0, nd0, nd, y); |
---|
717 | if (bd0 == NULL) |
---|
718 | goto ovfl; |
---|
719 | |
---|
720 | for(;;) { |
---|
721 | bd = Balloc(ptr,bd0->_k); |
---|
722 | if (bd == NULL) |
---|
723 | goto ovfl; |
---|
724 | Bcopy(bd, bd0); |
---|
725 | bb = d2b(ptr,dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */ |
---|
726 | if (bb == NULL) |
---|
727 | goto ovfl; |
---|
728 | bs = i2b(ptr,1); |
---|
729 | if (bs == NULL) |
---|
730 | goto ovfl; |
---|
731 | |
---|
732 | if (e >= 0) { |
---|
733 | bb2 = bb5 = 0; |
---|
734 | bd2 = bd5 = e; |
---|
735 | } |
---|
736 | else { |
---|
737 | bb2 = bb5 = -e; |
---|
738 | bd2 = bd5 = 0; |
---|
739 | } |
---|
740 | if (bbe >= 0) |
---|
741 | bb2 += bbe; |
---|
742 | else |
---|
743 | bd2 -= bbe; |
---|
744 | bs2 = bb2; |
---|
745 | #ifdef Honor_FLT_ROUNDS |
---|
746 | if (rounding != 1) |
---|
747 | bs2++; |
---|
748 | #endif |
---|
749 | #ifdef Avoid_Underflow |
---|
750 | Lsb = LSB; |
---|
751 | Lsb1 = 0; |
---|
752 | j = bbe - scale; |
---|
753 | i = j + bbbits - 1; /* logb(rv) */ |
---|
754 | j = P + 1 - bbbits; |
---|
755 | if (i < Emin) { /* denormal */ |
---|
756 | i = Emin - i; |
---|
757 | j -= i; |
---|
758 | if (i < 32) |
---|
759 | Lsb <<= i; |
---|
760 | else |
---|
761 | Lsb1 = Lsb << (i-32); |
---|
762 | } |
---|
763 | #else /*Avoid_Underflow*/ |
---|
764 | #ifdef Sudden_Underflow |
---|
765 | #ifdef IBM |
---|
766 | j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); |
---|
767 | #else |
---|
768 | j = P + 1 - bbbits; |
---|
769 | #endif |
---|
770 | #else /*Sudden_Underflow*/ |
---|
771 | j = bbe; |
---|
772 | i = j + bbbits - 1; /* logb(rv) */ |
---|
773 | if (i < Emin) /* denormal */ |
---|
774 | j += P - Emin; |
---|
775 | else |
---|
776 | j = P + 1 - bbbits; |
---|
777 | #endif /*Sudden_Underflow*/ |
---|
778 | #endif /*Avoid_Underflow*/ |
---|
779 | bb2 += j; |
---|
780 | bd2 += j; |
---|
781 | #ifdef Avoid_Underflow |
---|
782 | bd2 += scale; |
---|
783 | #endif |
---|
784 | i = bb2 < bd2 ? bb2 : bd2; |
---|
785 | if (i > bs2) |
---|
786 | i = bs2; |
---|
787 | if (i > 0) { |
---|
788 | bb2 -= i; |
---|
789 | bd2 -= i; |
---|
790 | bs2 -= i; |
---|
791 | } |
---|
792 | if (bb5 > 0) { |
---|
793 | bs = pow5mult(ptr, bs, bb5); |
---|
794 | if (bs == NULL) |
---|
795 | goto ovfl; |
---|
796 | bb1 = mult(ptr, bs, bb); |
---|
797 | if (bb1 == NULL) |
---|
798 | goto ovfl; |
---|
799 | Bfree(ptr, bb); |
---|
800 | bb = bb1; |
---|
801 | } |
---|
802 | if (bb2 > 0) { |
---|
803 | bb = lshift(ptr, bb, bb2); |
---|
804 | if (bb == NULL) |
---|
805 | goto ovfl; |
---|
806 | } |
---|
807 | if (bd5 > 0) { |
---|
808 | bd = pow5mult(ptr, bd, bd5); |
---|
809 | if (bd == NULL) |
---|
810 | goto ovfl; |
---|
811 | } |
---|
812 | if (bd2 > 0) { |
---|
813 | bd = lshift(ptr, bd, bd2); |
---|
814 | if (bd == NULL) |
---|
815 | goto ovfl; |
---|
816 | } |
---|
817 | if (bs2 > 0) { |
---|
818 | bs = lshift(ptr, bs, bs2); |
---|
819 | if (bs == NULL) |
---|
820 | goto ovfl; |
---|
821 | } |
---|
822 | delta = diff(ptr, bb, bd); |
---|
823 | if (delta == NULL) |
---|
824 | goto ovfl; |
---|
825 | dsign = delta->_sign; |
---|
826 | delta->_sign = 0; |
---|
827 | i = cmp(delta, bs); |
---|
828 | #ifdef Honor_FLT_ROUNDS |
---|
829 | if (rounding != 1) { |
---|
830 | if (i < 0) { |
---|
831 | /* Error is less than an ulp */ |
---|
832 | if (!delta->_x[0] && delta->_wds <= 1) { |
---|
833 | /* exact */ |
---|
834 | #ifdef SET_INEXACT |
---|
835 | inexact = 0; |
---|
836 | #endif |
---|
837 | break; |
---|
838 | } |
---|
839 | if (rounding) { |
---|
840 | if (dsign) { |
---|
841 | adj = 1.; |
---|
842 | goto apply_adj; |
---|
843 | } |
---|
844 | } |
---|
845 | else if (!dsign) { |
---|
846 | adj = -1.; |
---|
847 | if (!dword1(rv) |
---|
848 | && !(dword0(rv) & Frac_mask)) { |
---|
849 | y = dword0(rv) & Exp_mask; |
---|
850 | #ifdef Avoid_Underflow |
---|
851 | if (!scale || y > 2*P*Exp_msk1) |
---|
852 | #else |
---|
853 | if (y) |
---|
854 | #endif |
---|
855 | { |
---|
856 | delta = lshift(ptr, delta,Log2P); |
---|
857 | if (cmp(delta, bs) <= 0) |
---|
858 | adj = -0.5; |
---|
859 | } |
---|
860 | } |
---|
861 | apply_adj: |
---|
862 | #ifdef Avoid_Underflow |
---|
863 | if (scale && (y = dword0(rv) & Exp_mask) |
---|
864 | <= 2*P*Exp_msk1) |
---|
865 | dword0(adj) += (2*P+1)*Exp_msk1 - y; |
---|
866 | #else |
---|
867 | #ifdef Sudden_Underflow |
---|
868 | if ((dword0(rv) & Exp_mask) <= |
---|
869 | P*Exp_msk1) { |
---|
870 | dword0(rv) += P*Exp_msk1; |
---|
871 | dval(rv) += adj*ulp(dval(rv)); |
---|
872 | dword0(rv) -= P*Exp_msk1; |
---|
873 | } |
---|
874 | else |
---|
875 | #endif /*Sudden_Underflow*/ |
---|
876 | #endif /*Avoid_Underflow*/ |
---|
877 | dval(rv) += adj*ulp(dval(rv)); |
---|
878 | } |
---|
879 | break; |
---|
880 | } |
---|
881 | adj = ratio(delta, bs); |
---|
882 | if (adj < 1.) |
---|
883 | adj = 1.; |
---|
884 | if (adj <= 0x7ffffffe) { |
---|
885 | /* adj = rounding ? ceil(adj) : floor(adj); */ |
---|
886 | y = adj; |
---|
887 | if (y != adj) { |
---|
888 | if (!((rounding>>1) ^ dsign)) |
---|
889 | y++; |
---|
890 | adj = y; |
---|
891 | } |
---|
892 | } |
---|
893 | #ifdef Avoid_Underflow |
---|
894 | if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) |
---|
895 | dword0(adj) += (2*P+1)*Exp_msk1 - y; |
---|
896 | #else |
---|
897 | #ifdef Sudden_Underflow |
---|
898 | if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { |
---|
899 | dword0(rv) += P*Exp_msk1; |
---|
900 | adj *= ulp(dval(rv)); |
---|
901 | if (dsign) |
---|
902 | dval(rv) += adj; |
---|
903 | else |
---|
904 | dval(rv) -= adj; |
---|
905 | dword0(rv) -= P*Exp_msk1; |
---|
906 | goto cont; |
---|
907 | } |
---|
908 | #endif /*Sudden_Underflow*/ |
---|
909 | #endif /*Avoid_Underflow*/ |
---|
910 | adj *= ulp(dval(rv)); |
---|
911 | if (dsign) { |
---|
912 | if (dword0(rv) == Big0 && dword1(rv) == Big1) |
---|
913 | goto ovfl; |
---|
914 | dval(rv) += adj; |
---|
915 | else |
---|
916 | dval(rv) -= adj; |
---|
917 | goto cont; |
---|
918 | } |
---|
919 | #endif /*Honor_FLT_ROUNDS*/ |
---|
920 | |
---|
921 | if (i < 0) { |
---|
922 | /* Error is less than half an ulp -- check for |
---|
923 | * special case of mantissa a power of two. |
---|
924 | */ |
---|
925 | if (dsign || dword1(rv) || dword0(rv) & Bndry_mask |
---|
926 | #ifdef IEEE_Arith |
---|
927 | #ifdef Avoid_Underflow |
---|
928 | || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 |
---|
929 | #else |
---|
930 | || (dword0(rv) & Exp_mask) <= Exp_msk1 |
---|
931 | #endif |
---|
932 | #endif |
---|
933 | ) { |
---|
934 | #ifdef SET_INEXACT |
---|
935 | if (!delta->x[0] && delta->wds <= 1) |
---|
936 | inexact = 0; |
---|
937 | #endif |
---|
938 | break; |
---|
939 | } |
---|
940 | if (!delta->_x[0] && delta->_wds <= 1) { |
---|
941 | /* exact result */ |
---|
942 | #ifdef SET_INEXACT |
---|
943 | inexact = 0; |
---|
944 | #endif |
---|
945 | break; |
---|
946 | } |
---|
947 | delta = lshift(ptr,delta,Log2P); |
---|
948 | if (cmp(delta, bs) > 0) |
---|
949 | goto drop_down; |
---|
950 | break; |
---|
951 | } |
---|
952 | if (i == 0) { |
---|
953 | /* exactly half-way between */ |
---|
954 | if (dsign) { |
---|
955 | if ((dword0(rv) & Bndry_mask1) == Bndry_mask1 |
---|
956 | && dword1(rv) == ( |
---|
957 | #ifdef Avoid_Underflow |
---|
958 | (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) |
---|
959 | ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |
---|
960 | #endif |
---|
961 | 0xffffffff)) { |
---|
962 | /*boundary case -- increment exponent*/ |
---|
963 | if (dword0(rv) == Big0 && dword1(rv) == Big1) |
---|
964 | goto ovfl; |
---|
965 | dword0(rv) = (dword0(rv) & Exp_mask) |
---|
966 | + Exp_msk1 |
---|
967 | #ifdef IBM |
---|
968 | | Exp_msk1 >> 4 |
---|
969 | #endif |
---|
970 | ; |
---|
971 | #ifndef _DOUBLE_IS_32BITS |
---|
972 | dword1(rv) = 0; |
---|
973 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
974 | #ifdef Avoid_Underflow |
---|
975 | dsign = 0; |
---|
976 | #endif |
---|
977 | break; |
---|
978 | } |
---|
979 | } |
---|
980 | else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) { |
---|
981 | drop_down: |
---|
982 | /* boundary case -- decrement exponent */ |
---|
983 | #ifdef Sudden_Underflow /*{{*/ |
---|
984 | L = dword0(rv) & Exp_mask; |
---|
985 | #ifdef IBM |
---|
986 | if (L < Exp_msk1) |
---|
987 | #else |
---|
988 | #ifdef Avoid_Underflow |
---|
989 | if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) |
---|
990 | #else |
---|
991 | if (L <= Exp_msk1) |
---|
992 | #endif /*Avoid_Underflow*/ |
---|
993 | #endif /*IBM*/ |
---|
994 | goto undfl; |
---|
995 | L -= Exp_msk1; |
---|
996 | #else /*Sudden_Underflow}{*/ |
---|
997 | #ifdef Avoid_Underflow |
---|
998 | if (scale) { |
---|
999 | L = dword0(rv) & Exp_mask; |
---|
1000 | if (L <= (2*P+1)*Exp_msk1) { |
---|
1001 | if (L > (P+2)*Exp_msk1) |
---|
1002 | /* round even ==> */ |
---|
1003 | /* accept rv */ |
---|
1004 | break; |
---|
1005 | /* rv = smallest denormal */ |
---|
1006 | goto undfl; |
---|
1007 | } |
---|
1008 | } |
---|
1009 | #endif /*Avoid_Underflow*/ |
---|
1010 | L = (dword0(rv) & Exp_mask) - Exp_msk1; |
---|
1011 | #endif /*Sudden_Underflow}*/ |
---|
1012 | dword0(rv) = L | Bndry_mask1; |
---|
1013 | #ifndef _DOUBLE_IS_32BITS |
---|
1014 | dword1(rv) = 0xffffffff; |
---|
1015 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
1016 | #ifdef IBM |
---|
1017 | goto cont; |
---|
1018 | #else |
---|
1019 | break; |
---|
1020 | #endif |
---|
1021 | } |
---|
1022 | #ifndef ROUND_BIASED |
---|
1023 | #ifdef Avoid_Underflow |
---|
1024 | if (Lsb1) { |
---|
1025 | if (!(dword0(rv) & Lsb1)) |
---|
1026 | break; |
---|
1027 | } |
---|
1028 | else if (!(dword1(rv) & Lsb)) |
---|
1029 | break; |
---|
1030 | #else |
---|
1031 | if (!(dword1(rv) & LSB)) |
---|
1032 | break; |
---|
1033 | #endif |
---|
1034 | #endif |
---|
1035 | if (dsign) |
---|
1036 | #ifdef Avoid_Underflow |
---|
1037 | dval(rv) += sulp(rv, scale); |
---|
1038 | #else |
---|
1039 | dval(rv) += ulp(dval(rv)); |
---|
1040 | #endif |
---|
1041 | #ifndef ROUND_BIASED |
---|
1042 | else { |
---|
1043 | #ifdef Avoid_Underflow |
---|
1044 | dval(rv) -= sulp(rv, scale); |
---|
1045 | #else |
---|
1046 | dval(rv) -= ulp(dval(rv)); |
---|
1047 | #endif |
---|
1048 | #ifndef Sudden_Underflow |
---|
1049 | if (!dval(rv)) |
---|
1050 | goto undfl; |
---|
1051 | #endif |
---|
1052 | } |
---|
1053 | #ifdef Avoid_Underflow |
---|
1054 | dsign = 1 - dsign; |
---|
1055 | #endif |
---|
1056 | #endif |
---|
1057 | break; |
---|
1058 | } |
---|
1059 | if ((aadj = ratio(delta, bs)) <= 2.) { |
---|
1060 | if (dsign) |
---|
1061 | aadj = dval(aadj1) = 1.; |
---|
1062 | else if (dword1(rv) || dword0(rv) & Bndry_mask) { |
---|
1063 | #ifndef Sudden_Underflow |
---|
1064 | if (dword1(rv) == Tiny1 && !dword0(rv)) |
---|
1065 | goto undfl; |
---|
1066 | #endif |
---|
1067 | aadj = 1.; |
---|
1068 | dval(aadj1) = -1.; |
---|
1069 | } |
---|
1070 | else { |
---|
1071 | /* special case -- power of FLT_RADIX to be */ |
---|
1072 | /* rounded down... */ |
---|
1073 | |
---|
1074 | if (aadj < 2./FLT_RADIX) |
---|
1075 | aadj = 1./FLT_RADIX; |
---|
1076 | else |
---|
1077 | aadj *= 0.5; |
---|
1078 | dval(aadj1) = -aadj; |
---|
1079 | } |
---|
1080 | } |
---|
1081 | else { |
---|
1082 | aadj *= 0.5; |
---|
1083 | dval(aadj1) = dsign ? aadj : -aadj; |
---|
1084 | #ifdef Check_FLT_ROUNDS |
---|
1085 | switch(Rounding) { |
---|
1086 | case 2: /* towards +infinity */ |
---|
1087 | dval(aadj1) -= 0.5; |
---|
1088 | break; |
---|
1089 | case 0: /* towards 0 */ |
---|
1090 | case 3: /* towards -infinity */ |
---|
1091 | dval(aadj1) += 0.5; |
---|
1092 | } |
---|
1093 | #else |
---|
1094 | if (Flt_Rounds == 0) |
---|
1095 | dval(aadj1) += 0.5; |
---|
1096 | #endif /*Check_FLT_ROUNDS*/ |
---|
1097 | } |
---|
1098 | y = dword0(rv) & Exp_mask; |
---|
1099 | |
---|
1100 | /* Check for overflow */ |
---|
1101 | |
---|
1102 | if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { |
---|
1103 | dval(rv0) = dval(rv); |
---|
1104 | dword0(rv) -= P*Exp_msk1; |
---|
1105 | adj = dval(aadj1) * ulp(dval(rv)); |
---|
1106 | dval(rv) += adj; |
---|
1107 | if ((dword0(rv) & Exp_mask) >= |
---|
1108 | Exp_msk1*(DBL_MAX_EXP+Bias-P)) { |
---|
1109 | if (dword0(rv0) == Big0 && dword1(rv0) == Big1) |
---|
1110 | goto ovfl; |
---|
1111 | dword0(rv) = Big0; |
---|
1112 | #ifndef _DOUBLE_IS_32BITS |
---|
1113 | dword1(rv) = Big1; |
---|
1114 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
1115 | goto cont; |
---|
1116 | } |
---|
1117 | else |
---|
1118 | dword0(rv) += P*Exp_msk1; |
---|
1119 | } |
---|
1120 | else { |
---|
1121 | #ifdef Avoid_Underflow |
---|
1122 | if (scale && y <= 2*P*Exp_msk1) { |
---|
1123 | if (aadj <= 0x7fffffff) { |
---|
1124 | if ((z = aadj) == 0) |
---|
1125 | z = 1; |
---|
1126 | aadj = z; |
---|
1127 | dval(aadj1) = dsign ? aadj : -aadj; |
---|
1128 | } |
---|
1129 | dword0(aadj1) += (2*P+1)*Exp_msk1 - y; |
---|
1130 | } |
---|
1131 | adj = dval(aadj1) * ulp(dval(rv)); |
---|
1132 | dval(rv) += adj; |
---|
1133 | #else |
---|
1134 | #ifdef Sudden_Underflow |
---|
1135 | if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { |
---|
1136 | dval(rv0) = dval(rv); |
---|
1137 | dword0(rv) += P*Exp_msk1; |
---|
1138 | adj = dval(aadj1) * ulp(dval(rv)); |
---|
1139 | dval(rv) += adj; |
---|
1140 | #ifdef IBM |
---|
1141 | if ((dword0(rv) & Exp_mask) < P*Exp_msk1) |
---|
1142 | #else |
---|
1143 | if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) |
---|
1144 | #endif |
---|
1145 | { |
---|
1146 | if (dword0(rv0) == Tiny0 |
---|
1147 | && dword1(rv0) == Tiny1) |
---|
1148 | goto undfl; |
---|
1149 | #ifndef _DOUBLE_IS_32BITS |
---|
1150 | dword0(rv) = Tiny0; |
---|
1151 | dword1(rv) = Tiny1; |
---|
1152 | #else |
---|
1153 | dword0(rv) = Tiny1; |
---|
1154 | #endif /*_DOUBLE_IS_32BITS*/ |
---|
1155 | goto cont; |
---|
1156 | } |
---|
1157 | else |
---|
1158 | dword0(rv) -= P*Exp_msk1; |
---|
1159 | } |
---|
1160 | else { |
---|
1161 | adj = dval(aadj1) * ulp(dval(rv)); |
---|
1162 | dval(rv) += adj; |
---|
1163 | } |
---|
1164 | #else /*Sudden_Underflow*/ |
---|
1165 | /* Compute adj so that the IEEE rounding rules will |
---|
1166 | * correctly round rv + adj in some half-way cases. |
---|
1167 | * If rv * ulp(rv) is denormalized (i.e., |
---|
1168 | * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid |
---|
1169 | * trouble from bits lost to denormalization; |
---|
1170 | * example: 1.2e-307 . |
---|
1171 | */ |
---|
1172 | if (y <= (P-1)*Exp_msk1 && aadj > 1.) { |
---|
1173 | dval(aadj1) = (double)(int)(aadj + 0.5); |
---|
1174 | if (!dsign) |
---|
1175 | dval(aadj1) = -dval(aadj1); |
---|
1176 | } |
---|
1177 | adj = dval(aadj1) * ulp(dval(rv)); |
---|
1178 | dval(rv) += adj; |
---|
1179 | #endif /*Sudden_Underflow*/ |
---|
1180 | #endif /*Avoid_Underflow*/ |
---|
1181 | } |
---|
1182 | z = dword0(rv) & Exp_mask; |
---|
1183 | #ifndef SET_INEXACT |
---|
1184 | #ifdef Avoid_Underflow |
---|
1185 | if (!scale) |
---|
1186 | #endif |
---|
1187 | if (y == z) { |
---|
1188 | /* Can we stop now? */ |
---|
1189 | L = (Long)aadj; |
---|
1190 | aadj -= L; |
---|
1191 | /* The tolerances below are conservative. */ |
---|
1192 | if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) { |
---|
1193 | if (aadj < .4999999 || aadj > .5000001) |
---|
1194 | break; |
---|
1195 | } |
---|
1196 | else if (aadj < .4999999/FLT_RADIX) |
---|
1197 | break; |
---|
1198 | } |
---|
1199 | #endif |
---|
1200 | cont: |
---|
1201 | Bfree(ptr,bb); |
---|
1202 | Bfree(ptr,bd); |
---|
1203 | Bfree(ptr,bs); |
---|
1204 | Bfree(ptr,delta); |
---|
1205 | } |
---|
1206 | #ifdef SET_INEXACT |
---|
1207 | if (inexact) { |
---|
1208 | if (!oldinexact) { |
---|
1209 | dword0(rv0) = Exp_1 + (70 << Exp_shift); |
---|
1210 | #ifndef _DOUBLE_IS_32BITS |
---|
1211 | dword1(rv0) = 0; |
---|
1212 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
1213 | dval(rv0) += 1.; |
---|
1214 | } |
---|
1215 | } |
---|
1216 | else if (!oldinexact) |
---|
1217 | clear_inexact(); |
---|
1218 | #endif |
---|
1219 | #ifdef Avoid_Underflow |
---|
1220 | if (scale) { |
---|
1221 | dword0(rv0) = Exp_1 - 2*P*Exp_msk1; |
---|
1222 | #ifndef _DOUBLE_IS_32BITS |
---|
1223 | dword1(rv0) = 0; |
---|
1224 | #endif /*!_DOUBLE_IS_32BITS*/ |
---|
1225 | dval(rv) *= dval(rv0); |
---|
1226 | #ifndef NO_ERRNO |
---|
1227 | /* try to avoid the bug of testing an 8087 register value */ |
---|
1228 | if (dword0(rv) == 0 && dword1(rv) == 0) |
---|
1229 | ptr->_errno = ERANGE; |
---|
1230 | #endif |
---|
1231 | } |
---|
1232 | #endif /* Avoid_Underflow */ |
---|
1233 | #ifdef SET_INEXACT |
---|
1234 | if (inexact && !(dword0(rv) & Exp_mask)) { |
---|
1235 | /* set underflow bit */ |
---|
1236 | dval(rv0) = 1e-300; |
---|
1237 | dval(rv0) *= dval(rv0); |
---|
1238 | } |
---|
1239 | #endif |
---|
1240 | retfree: |
---|
1241 | Bfree(ptr,bb); |
---|
1242 | Bfree(ptr,bd); |
---|
1243 | Bfree(ptr,bs); |
---|
1244 | Bfree(ptr,bd0); |
---|
1245 | Bfree(ptr,delta); |
---|
1246 | ret: |
---|
1247 | if (se) |
---|
1248 | *se = (char *)s; |
---|
1249 | return sign ? -dval(rv) : dval(rv); |
---|
1250 | } |
---|
1251 | |
---|
1252 | double |
---|
1253 | _strtod_r (struct _reent *ptr, |
---|
1254 | const char *__restrict s00, |
---|
1255 | char **__restrict se) |
---|
1256 | { |
---|
1257 | return _strtod_l (ptr, s00, se, __get_current_locale ()); |
---|
1258 | } |
---|
1259 | |
---|
1260 | #ifndef _REENT_ONLY |
---|
1261 | |
---|
1262 | double |
---|
1263 | strtod_l (const char *__restrict s00, char **__restrict se, locale_t loc) |
---|
1264 | { |
---|
1265 | return _strtod_l (_REENT, s00, se, loc); |
---|
1266 | } |
---|
1267 | |
---|
1268 | double |
---|
1269 | strtod (const char *__restrict s00, char **__restrict se) |
---|
1270 | { |
---|
1271 | return _strtod_l (_REENT, s00, se, __get_current_locale ()); |
---|
1272 | } |
---|
1273 | |
---|
1274 | float |
---|
1275 | strtof_l (const char *__restrict s00, char **__restrict se, locale_t loc) |
---|
1276 | { |
---|
1277 | double val = _strtod_l (_REENT, s00, se, loc); |
---|
1278 | if (isnan (val)) |
---|
1279 | return nanf (NULL); |
---|
1280 | float retval = (float) val; |
---|
1281 | #ifndef NO_ERRNO |
---|
1282 | if (isinf (retval) && !isinf (val)) |
---|
1283 | _REENT->_errno = ERANGE; |
---|
1284 | #endif |
---|
1285 | return retval; |
---|
1286 | } |
---|
1287 | |
---|
1288 | float |
---|
1289 | strtof (const char *__restrict s00, |
---|
1290 | char **__restrict se) |
---|
1291 | { |
---|
1292 | double val = _strtod_l (_REENT, s00, se, __get_current_locale ()); |
---|
1293 | if (isnan (val)) |
---|
1294 | return nanf (NULL); |
---|
1295 | float retval = (float) val; |
---|
1296 | #ifndef NO_ERRNO |
---|
1297 | if (isinf (retval) && !isinf (val)) |
---|
1298 | _REENT->_errno = ERANGE; |
---|
1299 | #endif |
---|
1300 | return retval; |
---|
1301 | } |
---|
1302 | |
---|
1303 | #endif |
---|