1 /* $NetBSD: strtod.c,v 1.14 2013/05/17 12:55:57 joerg Exp $ */
2
3 /****************************************************************
4
5 The author of this software is David M. Gay.
6
7 Copyright (C) 1998-2001 by Lucent Technologies
8 All Rights Reserved
9
10 Permission to use, copy, modify, and distribute this software and
11 its documentation for any purpose and without fee is hereby
12 granted, provided that the above copyright notice appear in all
13 copies and that both that the copyright notice and this
14 permission notice and warranty disclaimer appear in supporting
15 documentation, and that the name of Lucent or any of its entities
16 not be used in advertising or publicity pertaining to
17 distribution of the software without specific, written prior
18 permission.
19
20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 THIS SOFTWARE.
28
29 ****************************************************************/
30
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
33
34 #include "namespace.h"
35 #include "gdtoaimp.h"
36 #ifndef NO_FENV_H
37 #include <fenv.h>
38 #endif
39
40 #ifdef USE_LOCALE
41 #include <locale.h>
42 #include "setlocale_local.h"
43 #endif
44
45 #ifdef IEEE_Arith
46 #ifndef NO_IEEE_Scale
47 #define Avoid_Underflow
48 #undef tinytens
49 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
50 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
51 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
52 9007199254740992.*9007199254740992.e-256
53 };
54 #endif
55 #endif
56
57 #ifdef Honor_FLT_ROUNDS
58 #undef Check_FLT_ROUNDS
59 #define Check_FLT_ROUNDS
60 #else
61 #define Rounding Flt_Rounds
62 #endif
63
64 #ifndef __HAVE_LONG_DOUBLE
65 __strong_alias(_strtold, strtod)
66 __weak_alias(strtold, _strtold)
67 __strong_alias(_strtold_l, strtod_l)
68 __weak_alias(strtold_l, _strtold_l)
69 #endif
70
71 #ifdef Avoid_Underflow /*{*/
72 static double
73 sulp
74 #ifdef KR_headers
75 (x, scale) U *x; int scale;
76 #else
77 (U *x, int scale)
78 #endif
79 {
80 U u;
81 double rv;
82 int i;
83
84 rv = ulp(x);
85 if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
86 return rv; /* Is there an example where i <= 0 ? */
87 word0(&u) = Exp_1 + (i << Exp_shift);
88 word1(&u) = 0;
89 return rv * u.d;
90 }
91 #endif /*}*/
92
93 static double
_int_strtod_l(CONST char * s00,char ** se,locale_t loc)94 _int_strtod_l(CONST char *s00, char **se, locale_t loc)
95 {
96 #ifdef Avoid_Underflow
97 int scale;
98 #endif
99 #ifdef INFNAN_CHECK
100 int decpt;
101 #endif
102 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
103 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
104 CONST char *s, *s0, *s1;
105 double aadj;
106 Long L;
107 U adj, aadj1, rv, rv0;
108 ULong y, z;
109 Bigint *bb = NULL, *bb1, *bd0;
110 Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
111 #ifdef Avoid_Underflow
112 ULong Lsb, Lsb1;
113 #endif
114 #ifdef SET_INEXACT
115 int inexact, oldinexact;
116 #endif
117 #ifdef USE_LOCALE /*{{*/
118 char *decimalpoint = localeconv_l(loc)->decimal_point;
119 size_t dplen = strlen(decimalpoint);
120 #endif /*USE_LOCALE}}*/
121
122 #ifdef Honor_FLT_ROUNDS /*{*/
123 int Rounding;
124 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
125 Rounding = Flt_Rounds;
126 #else /*}{*/
127 Rounding = 1;
128 switch(fegetround()) {
129 case FE_TOWARDZERO: Rounding = 0; break;
130 case FE_UPWARD: Rounding = 2; break;
131 case FE_DOWNWARD: Rounding = 3;
132 }
133 #endif /*}}*/
134 #endif /*}*/
135
136 #ifdef INFNAN_CHECK
137 decpt = 0;
138 #endif
139 sign = nz0 = nz = 0;
140 dval(&rv) = 0.;
141 for(s = s00;;s++) switch(*s) {
142 case '-':
143 sign = 1;
144 /* FALLTHROUGH */
145 case '+':
146 if (*++s)
147 goto break2;
148 /* FALLTHROUGH */
149 case 0:
150 goto ret0;
151 case '\t':
152 case '\n':
153 case '\v':
154 case '\f':
155 case '\r':
156 case ' ':
157 continue;
158 default:
159 goto break2;
160 }
161 break2:
162 if (*s == '0') {
163 #ifndef NO_HEX_FP /*{*/
164 {
165 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
166 Long expt;
167 ULong bits[2];
168 switch(s[1]) {
169 case 'x':
170 case 'X':
171 {
172 #ifdef Honor_FLT_ROUNDS
173 FPI fpi1 = fpi;
174 fpi1.rounding = Rounding;
175 #else
176 #define fpi1 fpi
177 #endif
178 switch((i = gethex(&s, &fpi1, &expt, &bb, sign, loc)) & STRTOG_Retmask) {
179 case STRTOG_NoNumber:
180 s = s00;
181 sign = 0;
182 /* FALLTHROUGH */
183 case STRTOG_Zero:
184 break;
185 default:
186 if (bb) {
187 copybits(bits, fpi.nbits, bb);
188 Bfree(bb);
189 }
190 ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
191 }}
192 goto ret;
193 }
194 }
195 #endif /*}*/
196 nz0 = 1;
197 while(*++s == '0') ;
198 if (!*s)
199 goto ret;
200 }
201 s0 = s;
202 y = z = 0;
203 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
204 if (nd < 9)
205 y = 10*y + c - '0';
206 else if (nd < 16)
207 z = 10*z + c - '0';
208 nd0 = nd;
209 #ifdef USE_LOCALE
210 if (c == *decimalpoint) {
211 for(i = 1; decimalpoint[i]; ++i)
212 if (s[i] != decimalpoint[i])
213 goto dig_done;
214 s += i;
215 c = *s;
216 #else
217 if (c == '.') {
218 c = *++s;
219 #endif
220 #ifdef INFNAN_CHECK
221 decpt = 1;
222 #endif
223 if (!nd) {
224 for(; c == '0'; c = *++s)
225 nz++;
226 if (c > '0' && c <= '9') {
227 s0 = s;
228 nf += nz;
229 nz = 0;
230 goto have_dig;
231 }
232 goto dig_done;
233 }
234 for(; c >= '0' && c <= '9'; c = *++s) {
235 have_dig:
236 nz++;
237 if (c -= '0') {
238 nf += nz;
239 for(i = 1; i < nz; i++)
240 if (nd++ < 9)
241 y *= 10;
242 else if (nd <= DBL_DIG + 1)
243 z *= 10;
244 if (nd++ < 9)
245 y = 10*y + c;
246 else if (nd <= DBL_DIG + 1)
247 z = 10*z + c;
248 nz = 0;
249 }
250 }
251 }/*}*/
252 dig_done:
253 e = 0;
254 if (c == 'e' || c == 'E') {
255 if (!nd && !nz && !nz0) {
256 goto ret0;
257 }
258 s00 = s;
259 esign = 0;
260 switch(c = *++s) {
261 case '-':
262 esign = 1;
263 /* FALLTHROUGH */
264 case '+':
265 c = *++s;
266 }
267 if (c >= '0' && c <= '9') {
268 while(c == '0')
269 c = *++s;
270 if (c > '0' && c <= '9') {
271 L = c - '0';
272 s1 = s;
273 while((c = *++s) >= '0' && c <= '9')
274 L = 10*L + c - '0';
275 if (s - s1 > 8 || L > 19999)
276 /* Avoid confusion from exponents
277 * so large that e might overflow.
278 */
279 e = 19999; /* safe for 16 bit ints */
280 else
281 e = (int)L;
282 if (esign)
283 e = -e;
284 }
285 else
286 e = 0;
287 }
288 else
289 s = s00;
290 }
291 if (!nd) {
292 if (!nz && !nz0) {
293 #ifdef INFNAN_CHECK
294 /* Check for Nan and Infinity */
295 ULong bits[2];
296 static FPI fpinan = /* only 52 explicit bits */
297 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
298 if (!decpt)
299 switch(c) {
300 case 'i':
301 case 'I':
302 if (match(&s,"nf")) {
303 --s;
304 if (!match(&s,"inity"))
305 ++s;
306 word0(&rv) = 0x7ff00000;
307 word1(&rv) = 0;
308 goto ret;
309 }
310 break;
311 case 'n':
312 case 'N':
313 if (match(&s, "an")) {
314 #ifndef No_Hex_NaN
315 if (*s == '(' /*)*/
316 && hexnan(&s, &fpinan, bits)
317 == STRTOG_NaNbits) {
318 word0(&rv) = 0x7ff00000 | bits[1];
319 word1(&rv) = bits[0];
320 }
321 else {
322 #endif
323 word0(&rv) = NAN_WORD0;
324 word1(&rv) = NAN_WORD1;
325 #ifndef No_Hex_NaN
326 }
327 #endif
328 goto ret;
329 }
330 }
331 #endif /* INFNAN_CHECK */
332 ret0:
333 s = s00;
334 sign = 0;
335 }
336 goto ret;
337 }
338 e1 = e -= nf;
339
340 /* Now we have nd0 digits, starting at s0, followed by a
341 * decimal point, followed by nd-nd0 digits. The number we're
342 * after is the integer represented by those digits times
343 * 10**e */
344
345 if (!nd0)
346 nd0 = nd;
347 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
348 dval(&rv) = y;
349 if (k > 9) {
350 #ifdef SET_INEXACT
351 if (k > DBL_DIG)
352 oldinexact = get_inexact();
353 #endif
354 dval(&rv) = tens[k - 9] * dval(&rv) + z;
355 }
356 bd0 = 0;
357 if (nd <= DBL_DIG
358 #ifndef RND_PRODQUOT
359 #ifndef Honor_FLT_ROUNDS
360 && Flt_Rounds == 1
361 #endif
362 #endif
363 ) {
364 if (!e)
365 goto ret;
366 #ifndef ROUND_BIASED_without_Round_Up
367 if (e > 0) {
368 if (e <= Ten_pmax) {
369 #ifdef VAX
370 goto vax_ovfl_check;
371 #else
372 #ifdef Honor_FLT_ROUNDS
373 /* round correctly FLT_ROUNDS = 2 or 3 */
374 if (sign) {
375 rv.d = -rv.d;
376 sign = 0;
377 }
378 #endif
379 /* rv = */ rounded_product(dval(&rv), tens[e]);
380 goto ret;
381 #endif
382 }
383 i = DBL_DIG - nd;
384 if (e <= Ten_pmax + i) {
385 /* A fancier test would sometimes let us do
386 * this for larger i values.
387 */
388 #ifdef Honor_FLT_ROUNDS
389 /* round correctly FLT_ROUNDS = 2 or 3 */
390 if (sign) {
391 rv.d = -rv.d;
392 sign = 0;
393 }
394 #endif
395 e -= i;
396 dval(&rv) *= tens[i];
397 #ifdef VAX
398 /* VAX exponent range is so narrow we must
399 * worry about overflow here...
400 */
401 vax_ovfl_check:
402 word0(&rv) -= P*Exp_msk1;
403 /* rv = */ rounded_product(dval(&rv), tens[e]);
404 if ((word0(&rv) & Exp_mask)
405 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
406 goto ovfl;
407 word0(&rv) += P*Exp_msk1;
408 #else
409 /* rv = */ rounded_product(dval(&rv), tens[e]);
410 #endif
411 goto ret;
412 }
413 }
414 #ifndef Inaccurate_Divide
415 else if (e >= -Ten_pmax) {
416 #ifdef Honor_FLT_ROUNDS
417 /* round correctly FLT_ROUNDS = 2 or 3 */
418 if (sign) {
419 rv.d = -rv.d;
420 sign = 0;
421 }
422 #endif
423 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
424 goto ret;
425 }
426 #endif
427 #endif /* ROUND_BIASED_without_Round_Up */
428 }
429 e1 += nd - k;
430
431 #ifdef IEEE_Arith
432 #ifdef SET_INEXACT
433 inexact = 1;
434 if (k <= DBL_DIG)
435 oldinexact = get_inexact();
436 #endif
437 #ifdef Avoid_Underflow
438 scale = 0;
439 #endif
440 #ifdef Honor_FLT_ROUNDS
441 if (Rounding >= 2) {
442 if (sign)
443 Rounding = Rounding == 2 ? 0 : 2;
444 else
445 if (Rounding != 2)
446 Rounding = 0;
447 }
448 #endif
449 #endif /*IEEE_Arith*/
450
451 /* Get starting approximation = rv * 10**e1 */
452
453 if (e1 > 0) {
454 if ( (i = e1 & 15) !=0)
455 dval(&rv) *= tens[i];
456 if (e1 &= ~15) {
457 if (e1 > DBL_MAX_10_EXP) {
458 ovfl:
459 /* Can't trust HUGE_VAL */
460 #ifdef IEEE_Arith
461 #ifdef Honor_FLT_ROUNDS
462 switch(Rounding) {
463 case 0: /* toward 0 */
464 case 3: /* toward -infinity */
465 word0(&rv) = Big0;
466 word1(&rv) = Big1;
467 break;
468 default:
469 word0(&rv) = Exp_mask;
470 word1(&rv) = 0;
471 }
472 #else /*Honor_FLT_ROUNDS*/
473 word0(&rv) = Exp_mask;
474 word1(&rv) = 0;
475 #endif /*Honor_FLT_ROUNDS*/
476 #ifdef SET_INEXACT
477 /* set overflow bit */
478 dval(&rv0) = 1e300;
479 dval(&rv0) *= dval(&rv0);
480 #endif
481 #else /*IEEE_Arith*/
482 word0(&rv) = Big0;
483 word1(&rv) = Big1;
484 #endif /*IEEE_Arith*/
485 range_err:
486 if (bd0) {
487 Bfree(bb);
488 Bfree(bd);
489 Bfree(bs);
490 Bfree(bd0);
491 Bfree(delta);
492 }
493 #ifndef NO_ERRNO
494 errno = ERANGE;
495 #endif
496 goto ret;
497 }
498 e1 = (unsigned int)e1 >> 4;
499 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
500 if (e1 & 1)
501 dval(&rv) *= bigtens[j];
502 /* The last multiplication could overflow. */
503 word0(&rv) -= P*Exp_msk1;
504 dval(&rv) *= bigtens[j];
505 if ((z = word0(&rv) & Exp_mask)
506 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
507 goto ovfl;
508 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
509 /* set to largest number */
510 /* (Can't trust DBL_MAX) */
511 word0(&rv) = Big0;
512 word1(&rv) = Big1;
513 }
514 else
515 word0(&rv) += P*Exp_msk1;
516 }
517 }
518 else if (e1 < 0) {
519 e1 = -e1;
520 if ( (i = e1 & 15) !=0)
521 dval(&rv) /= tens[i];
522 if (e1 >>= 4) {
523 if (e1 >= 1 << n_bigtens)
524 goto undfl;
525 #ifdef Avoid_Underflow
526 if (e1 & Scale_Bit)
527 scale = 2*P;
528 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
529 if (e1 & 1)
530 dval(&rv) *= tinytens[j];
531 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
532 >> Exp_shift)) > 0) {
533 /* scaled rv is denormal; zap j low bits */
534 if (j >= 32) {
535 word1(&rv) = 0;
536 if (j >= 53)
537 word0(&rv) = (P+2)*Exp_msk1;
538 else
539 word0(&rv) &= 0xffffffffU << (j-32);
540 }
541 else
542 word1(&rv) &= 0xffffffffU << j;
543 }
544 #else
545 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
546 if (e1 & 1)
547 dval(&rv) *= tinytens[j];
548 /* The last multiplication could underflow. */
549 dval(&rv0) = dval(&rv);
550 dval(&rv) *= tinytens[j];
551 if (!dval(&rv)) {
552 dval(&rv) = 2.*dval(&rv0);
553 dval(&rv) *= tinytens[j];
554 #endif
555 if (!dval(&rv)) {
556 undfl:
557 dval(&rv) = 0.;
558 goto range_err;
559 }
560 #ifndef Avoid_Underflow
561 word0(&rv) = Tiny0;
562 word1(&rv) = Tiny1;
563 /* The refinement below will clean
564 * this approximation up.
565 */
566 }
567 #endif
568 }
569 }
570
571 /* Now the hard part -- adjusting rv to the correct value.*/
572
573 /* Put digits into bd: true value = bd * 10^e */
574
575 bd0 = s2b(s0, nd0, nd, y, dplen);
576 if (bd0 == NULL)
577 goto ovfl;
578
579 for(;;) {
580 bd = Balloc(bd0->k);
581 if (bd == NULL)
582 goto ovfl;
583 Bcopy(bd, bd0);
584 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
585 if (bb == NULL)
586 goto ovfl;
587 bs = i2b(1);
588 if (bs == NULL)
589 goto ovfl;
590
591 if (e >= 0) {
592 bb2 = bb5 = 0;
593 bd2 = bd5 = e;
594 }
595 else {
596 bb2 = bb5 = -e;
597 bd2 = bd5 = 0;
598 }
599 if (bbe >= 0)
600 bb2 += bbe;
601 else
602 bd2 -= bbe;
603 bs2 = bb2;
604 #ifdef Honor_FLT_ROUNDS
605 if (Rounding != 1)
606 bs2++;
607 #endif
608 #ifdef Avoid_Underflow
609 Lsb = LSB;
610 Lsb1 = 0;
611 j = bbe - scale;
612 i = j + bbbits - 1; /* logb(rv) */
613 j = P + 1 - bbbits;
614 if (i < Emin) { /* denormal */
615 i = Emin - i;
616 j -= i;
617 if (i < 32)
618 Lsb <<= i;
619 else
620 Lsb1 = Lsb << (i-32);
621 }
622 #else /*Avoid_Underflow*/
623 #ifdef Sudden_Underflow
624 #ifdef IBM
625 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
626 #else
627 j = P + 1 - bbbits;
628 #endif
629 #else /*Sudden_Underflow*/
630 j = bbe;
631 i = j + bbbits - 1; /* logb(&rv) */
632 if (i < Emin) /* denormal */
633 j += P - Emin;
634 else
635 j = P + 1 - bbbits;
636 #endif /*Sudden_Underflow*/
637 #endif /*Avoid_Underflow*/
638 bb2 += j;
639 bd2 += j;
640 #ifdef Avoid_Underflow
641 bd2 += scale;
642 #endif
643 i = bb2 < bd2 ? bb2 : bd2;
644 if (i > bs2)
645 i = bs2;
646 if (i > 0) {
647 bb2 -= i;
648 bd2 -= i;
649 bs2 -= i;
650 }
651 if (bb5 > 0) {
652 bs = pow5mult(bs, bb5);
653 if (bs == NULL)
654 goto ovfl;
655 bb1 = mult(bs, bb);
656 if (bb1 == NULL)
657 goto ovfl;
658 Bfree(bb);
659 bb = bb1;
660 }
661 if (bb2 > 0) {
662 bb = lshift(bb, bb2);
663 if (bb == NULL)
664 goto ovfl;
665 }
666 if (bd5 > 0) {
667 bd = pow5mult(bd, bd5);
668 if (bd == NULL)
669 goto ovfl;
670 }
671 if (bd2 > 0) {
672 bd = lshift(bd, bd2);
673 if (bd == NULL)
674 goto ovfl;
675 }
676 if (bs2 > 0) {
677 bs = lshift(bs, bs2);
678 if (bs == NULL)
679 goto ovfl;
680 }
681 delta = diff(bb, bd);
682 if (delta == NULL)
683 goto ovfl;
684 dsign = delta->sign;
685 delta->sign = 0;
686 i = cmp(delta, bs);
687 #ifdef Honor_FLT_ROUNDS
688 if (Rounding != 1) {
689 if (i < 0) {
690 /* Error is less than an ulp */
691 if (!delta->x[0] && delta->wds <= 1) {
692 /* exact */
693 #ifdef SET_INEXACT
694 inexact = 0;
695 #endif
696 break;
697 }
698 if (Rounding) {
699 if (dsign) {
700 dval(&adj) = 1.;
701 goto apply_adj;
702 }
703 }
704 else if (!dsign) {
705 dval(&adj) = -1.;
706 if (!word1(&rv)
707 && !(word0(&rv) & Frac_mask)) {
708 y = word0(&rv) & Exp_mask;
709 #ifdef Avoid_Underflow
710 if (!scale || y > 2*P*Exp_msk1)
711 #else
712 if (y)
713 #endif
714 {
715 delta = lshift(delta,Log2P);
716 if (cmp(delta, bs) <= 0)
717 dval(&adj) = -0.5;
718 }
719 }
720 apply_adj:
721 #ifdef Avoid_Underflow
722 if (scale && (y = word0(&rv) & Exp_mask)
723 <= 2*P*Exp_msk1)
724 word0(&adj) += (2*P+1)*Exp_msk1 - y;
725 #else
726 #ifdef Sudden_Underflow
727 if ((word0(&rv) & Exp_mask) <=
728 P*Exp_msk1) {
729 word0(&rv) += P*Exp_msk1;
730 dval(&rv) += adj*ulp(&rv);
731 word0(&rv) -= P*Exp_msk1;
732 }
733 else
734 #endif /*Sudden_Underflow*/
735 #endif /*Avoid_Underflow*/
736 dval(&rv) += adj.d*ulp(&rv);
737 }
738 break;
739 }
740 dval(&adj) = ratio(delta, bs);
741 if (adj.d < 1.)
742 dval(&adj) = 1.;
743 if (adj.d <= 0x7ffffffe) {
744 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
745 y = adj.d;
746 if (y != adj.d) {
747 if (!((Rounding>>1) ^ dsign))
748 y++;
749 dval(&adj) = y;
750 }
751 }
752 #ifdef Avoid_Underflow
753 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
754 word0(&adj) += (2*P+1)*Exp_msk1 - y;
755 #else
756 #ifdef Sudden_Underflow
757 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
758 word0(&rv) += P*Exp_msk1;
759 dval(&adj) *= ulp(&rv);
760 if (dsign)
761 dval(&rv) += adj;
762 else
763 dval(&rv) -= adj;
764 word0(&rv) -= P*Exp_msk1;
765 goto cont;
766 }
767 #endif /*Sudden_Underflow*/
768 #endif /*Avoid_Underflow*/
769 dval(&adj) *= ulp(&rv);
770 if (dsign) {
771 if (word0(&rv) == Big0 && word1(&rv) == Big1)
772 goto ovfl;
773 dval(&rv) += adj.d;
774 }
775 else
776 dval(&rv) -= adj.d;
777 goto cont;
778 }
779 #endif /*Honor_FLT_ROUNDS*/
780
781 if (i < 0) {
782 /* Error is less than half an ulp -- check for
783 * special case of mantissa a power of two.
784 */
785 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
786 #ifdef IEEE_Arith
787 #ifdef Avoid_Underflow
788 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
789 #else
790 || (word0(&rv) & Exp_mask) <= Exp_msk1
791 #endif
792 #endif
793 ) {
794 #ifdef SET_INEXACT
795 if (!delta->x[0] && delta->wds <= 1)
796 inexact = 0;
797 #endif
798 break;
799 }
800 if (!delta->x[0] && delta->wds <= 1) {
801 /* exact result */
802 #ifdef SET_INEXACT
803 inexact = 0;
804 #endif
805 break;
806 }
807 delta = lshift(delta,Log2P);
808 if (cmp(delta, bs) > 0)
809 goto drop_down;
810 break;
811 }
812 if (i == 0) {
813 /* exactly half-way between */
814 if (dsign) {
815 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
816 && word1(&rv) == (
817 #ifdef Avoid_Underflow
818 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
819 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
820 #endif
821 0xffffffff)) {
822 /*boundary case -- increment exponent*/
823 if (word0(&rv) == Big0 && word1(&rv) == Big1)
824 goto ovfl;
825 word0(&rv) = (word0(&rv) & Exp_mask)
826 + Exp_msk1
827 #ifdef IBM
828 | Exp_msk1 >> 4
829 #endif
830 ;
831 word1(&rv) = 0;
832 #ifdef Avoid_Underflow
833 dsign = 0;
834 #endif
835 break;
836 }
837 }
838 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
839 drop_down:
840 /* boundary case -- decrement exponent */
841 #ifdef Sudden_Underflow /*{{*/
842 L = word0(&rv) & Exp_mask;
843 #ifdef IBM
844 if (L < Exp_msk1)
845 #else
846 #ifdef Avoid_Underflow
847 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
848 #else
849 if (L <= Exp_msk1)
850 #endif /*Avoid_Underflow*/
851 #endif /*IBM*/
852 goto undfl;
853 L -= Exp_msk1;
854 #else /*Sudden_Underflow}{*/
855 #ifdef Avoid_Underflow
856 if (scale) {
857 L = word0(&rv) & Exp_mask;
858 if (L <= (2*P+1)*Exp_msk1) {
859 if (L > (P+2)*Exp_msk1)
860 /* round even ==> */
861 /* accept rv */
862 break;
863 /* rv = smallest denormal */
864 goto undfl;
865 }
866 }
867 #endif /*Avoid_Underflow*/
868 L = (word0(&rv) & Exp_mask) - Exp_msk1;
869 #endif /*Sudden_Underflow}}*/
870 word0(&rv) = L | Bndry_mask1;
871 word1(&rv) = 0xffffffff;
872 #ifdef IBM
873 goto cont;
874 #else
875 break;
876 #endif
877 }
878 #ifndef ROUND_BIASED
879 #ifdef Avoid_Underflow
880 if (Lsb1) {
881 if (!(word0(&rv) & Lsb1))
882 break;
883 }
884 else if (!(word1(&rv) & Lsb))
885 break;
886 #else
887 if (!(word1(&rv) & LSB))
888 break;
889 #endif
890 #endif
891 if (dsign)
892 #ifdef Avoid_Underflow
893 dval(&rv) += sulp(&rv, scale);
894 #else
895 dval(&rv) += ulp(&rv);
896 #endif
897 #ifndef ROUND_BIASED
898 else {
899 #ifdef Avoid_Underflow
900 dval(&rv) -= sulp(&rv, scale);
901 #else
902 dval(&rv) -= ulp(&rv);
903 #endif
904 #ifndef Sudden_Underflow
905 if (!dval(&rv))
906 goto undfl;
907 #endif
908 }
909 #ifdef Avoid_Underflow
910 dsign = 1 - dsign;
911 #endif
912 #endif
913 break;
914 }
915 if ((aadj = ratio(delta, bs)) <= 2.) {
916 if (dsign)
917 aadj = dval(&aadj1) = 1.;
918 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
919 #ifndef Sudden_Underflow
920 if (word1(&rv) == Tiny1 && !word0(&rv))
921 goto undfl;
922 #endif
923 aadj = 1.;
924 dval(&aadj1) = -1.;
925 }
926 else {
927 /* special case -- power of FLT_RADIX to be */
928 /* rounded down... */
929
930 if (aadj < 2./FLT_RADIX)
931 aadj = 1./FLT_RADIX;
932 else
933 aadj *= 0.5;
934 dval(&aadj1) = -aadj;
935 }
936 }
937 else {
938 aadj *= 0.5;
939 dval(&aadj1) = dsign ? aadj : -aadj;
940 #ifdef Check_FLT_ROUNDS
941 /* CONSTCOND */
942 switch(Rounding) {
943 case 2: /* towards +infinity */
944 dval(&aadj1) -= 0.5;
945 break;
946 case 0: /* towards 0 */
947 case 3: /* towards -infinity */
948 dval(&aadj1) += 0.5;
949 }
950 #else
951 /* CONSTCOND */
952 if (Flt_Rounds == 0)
953 dval(&aadj1) += 0.5;
954 #endif /*Check_FLT_ROUNDS*/
955 }
956 y = word0(&rv) & Exp_mask;
957
958 /* Check for overflow */
959
960 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
961 dval(&rv0) = dval(&rv);
962 word0(&rv) -= P*Exp_msk1;
963 dval(&adj) = dval(&aadj1) * ulp(&rv);
964 dval(&rv) += dval(&adj);
965 if ((word0(&rv) & Exp_mask) >=
966 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
967 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
968 goto ovfl;
969 word0(&rv) = Big0;
970 word1(&rv) = Big1;
971 goto cont;
972 }
973 else
974 word0(&rv) += P*Exp_msk1;
975 }
976 else {
977 #ifdef Avoid_Underflow
978 if (scale && y <= 2*P*Exp_msk1) {
979 if (aadj <= 0x7fffffff) {
980 if ((z = aadj) == 0)
981 z = 1;
982 aadj = z;
983 dval(&aadj1) = dsign ? aadj : -aadj;
984 }
985 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
986 }
987 dval(&adj) = dval(&aadj1) * ulp(&rv);
988 dval(&rv) += dval(&adj);
989 #else
990 #ifdef Sudden_Underflow
991 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
992 dval(&rv0) = dval(&rv);
993 word0(&rv) += P*Exp_msk1;
994 dval(&adj) = dval(&aadj1) * ulp(&rv);
995 dval(&rv) += dval(&adj);
996 #ifdef IBM
997 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
998 #else
999 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
1000 #endif
1001 {
1002 if (word0(&rv0) == Tiny0
1003 && word1(&rv0) == Tiny1)
1004 goto undfl;
1005 word0(&rv) = Tiny0;
1006 word1(&rv) = Tiny1;
1007 goto cont;
1008 }
1009 else
1010 word0(&rv) -= P*Exp_msk1;
1011 }
1012 else {
1013 dval(&adj) = dval(&aadj1) * ulp(&rv);
1014 dval(&rv) += dval(&adj);
1015 }
1016 #else /*Sudden_Underflow*/
1017 /* Compute dval(&adj) so that the IEEE rounding rules will
1018 * correctly round rv + dval(&adj) in some half-way cases.
1019 * If rv * ulp(&rv) is denormalized (i.e.,
1020 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1021 * trouble from bits lost to denormalization;
1022 * example: 1.2e-307 .
1023 */
1024 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1025 dval(&aadj1) = (double)(int)(aadj + 0.5);
1026 if (!dsign)
1027 dval(&aadj1) = -dval(&aadj1);
1028 }
1029 dval(&adj) = dval(&aadj1) * ulp(&rv);
1030 dval(&rv) += adj;
1031 #endif /*Sudden_Underflow*/
1032 #endif /*Avoid_Underflow*/
1033 }
1034 z = word0(&rv) & Exp_mask;
1035 #ifndef SET_INEXACT
1036 #ifdef Avoid_Underflow
1037 if (!scale)
1038 #endif
1039 if (y == z) {
1040 /* Can we stop now? */
1041 L = (Long)aadj;
1042 aadj -= L;
1043 /* The tolerances below are conservative. */
1044 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1045 if (aadj < .4999999 || aadj > .5000001)
1046 break;
1047 }
1048 else if (aadj < .4999999/FLT_RADIX)
1049 break;
1050 }
1051 #endif
1052 cont:
1053 Bfree(bb);
1054 Bfree(bd);
1055 Bfree(bs);
1056 Bfree(delta);
1057 }
1058 Bfree(bb);
1059 Bfree(bd);
1060 Bfree(bs);
1061 Bfree(bd0);
1062 Bfree(delta);
1063 #ifdef SET_INEXACT
1064 if (inexact) {
1065 if (!oldinexact) {
1066 word0(&rv0) = Exp_1 + (70 << Exp_shift);
1067 word1(&rv0) = 0;
1068 dval(&rv0) += 1.;
1069 }
1070 }
1071 else if (!oldinexact)
1072 clear_inexact();
1073 #endif
1074 #ifdef Avoid_Underflow
1075 if (scale) {
1076 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1077 word1(&rv0) = 0;
1078 dval(&rv) *= dval(&rv0);
1079 #ifndef NO_ERRNO
1080 /* try to avoid the bug of testing an 8087 register value */
1081 #ifdef IEEE_Arith
1082 if (!(word0(&rv) & Exp_mask))
1083 #else
1084 if (word0(&rv) == 0 && word1(&rv) == 0)
1085 #endif
1086 errno = ERANGE;
1087 #endif
1088 }
1089 #endif /* Avoid_Underflow */
1090 #ifdef SET_INEXACT
1091 if (inexact && !(word0(&rv) & Exp_mask)) {
1092 /* set underflow bit */
1093 dval(&rv0) = 1e-300;
1094 dval(&rv0) *= dval(&rv0);
1095 }
1096 #endif
1097 ret:
1098 if (se)
1099 *se = __UNCONST(s);
1100 return sign ? -dval(&rv) : dval(&rv);
1101 }
1102
1103 double
1104 strtod(CONST char *s, char **sp)
1105 {
1106 return _int_strtod_l(s, sp, _current_locale());
1107 }
1108
1109 #ifdef __weak_alias
1110 __weak_alias(strtod_l, _strtod_l)
1111 #endif
1112
1113 double
1114 strtod_l(CONST char *s, char **sp, locale_t loc)
1115 {
1116 return _int_strtod_l(s, sp, loc);
1117 }
1118