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