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 #include "gdtoaimp.h"
33 #ifndef NO_FENV_H
34 #include <fenv.h>
35 #endif
36
37 #ifdef USE_LOCALE
38 #include "locale.h"
39 #endif
40
41 #ifdef IEEE_Arith
42 #ifndef NO_IEEE_Scale
43 #define Avoid_Underflow
44 #undef tinytens
45 /* The factor of 2^106 in tinytens[4] helps us avoid setting the underflow */
46 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
47 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
48 9007199254740992.*9007199254740992.e-256
49 };
50 #endif
51 #endif
52
53 #ifdef Honor_FLT_ROUNDS
54 #undef Check_FLT_ROUNDS
55 #define Check_FLT_ROUNDS
56 #else
57 #define Rounding Flt_Rounds
58 #endif
59
60 #ifdef Avoid_Underflow /*{*/
61 static double
sulp(x,scale)62 sulp
63 #ifdef KR_headers
64 (x, scale) U *x; int scale;
65 #else
66 (U *x, int scale)
67 #endif
68 {
69 U u;
70 double rv;
71 int i;
72
73 rv = ulp(x);
74 if (!scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
75 return rv; /* Is there an example where i <= 0 ? */
76 word0(&u) = Exp_1 + (i << Exp_shift);
77 word1(&u) = 0;
78 return rv * u.d;
79 }
80 #endif /*}*/
81
82 double
strtod(s00,se)83 strtod
84 #ifdef KR_headers
85 (s00, se) CONST char *s00; char **se;
86 #else
87 (CONST char *s00, char **se)
88 #endif
89 {
90 #ifdef Avoid_Underflow
91 int scale;
92 #endif
93 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
94 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
95 CONST char *s, *s0, *s1;
96 double aadj;
97 Long L;
98 U adj, aadj1, rv, rv0;
99 ULong y, z;
100 Bigint *bb = NULL, *bb1, *bd = NULL, *bd0 = NULL, *bs = NULL, *delta = NULL;
101 #ifdef Avoid_Underflow
102 ULong Lsb, Lsb1;
103 #endif
104 #ifdef SET_INEXACT
105 int inexact, oldinexact;
106 #endif
107 #ifdef USE_LOCALE /*{{*/
108 #ifdef NO_LOCALE_CACHE
109 char *decimalpoint = localeconv()->decimal_point;
110 int dplen = strlen(decimalpoint);
111 #else
112 char *decimalpoint;
113 static char *decimalpoint_cache;
114 static int dplen;
115 if (!(s0 = decimalpoint_cache)) {
116 s0 = localeconv()->decimal_point;
117 decimalpoint_cache = strdup(s0);
118 dplen = strlen(s0);
119 }
120 decimalpoint = (char*)s0;
121 #endif /*NO_LOCALE_CACHE*/
122 #else /*USE_LOCALE}{*/
123 #define dplen 1
124 #endif /*USE_LOCALE}}*/
125
126 #ifdef Honor_FLT_ROUNDS /*{*/
127 int Rounding;
128 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
129 Rounding = Flt_Rounds;
130 #else /*}{*/
131 Rounding = 1;
132 switch(fegetround()) {
133 case FE_TOWARDZERO: Rounding = 0; break;
134 case FE_UPWARD: Rounding = 2; break;
135 case FE_DOWNWARD: Rounding = 3;
136 }
137 #endif /*}}*/
138 #endif /*}*/
139
140 sign = nz0 = nz = decpt = 0;
141 dval(&rv) = 0.;
142 for(s = s00;;s++) switch(*s) {
143 case '-':
144 sign = 1;
145 /* no break */
146 case '+':
147 if (*++s)
148 goto break2;
149 /* no break */
150 case 0:
151 goto ret0;
152 case '\t':
153 case '\n':
154 case '\v':
155 case '\f':
156 case '\r':
157 case ' ':
158 continue;
159 default:
160 goto break2;
161 }
162 break2:
163 if (*s == '0') {
164 #ifndef NO_HEX_FP /*{*/
165 {
166 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
167 Long exp;
168 ULong bits[2];
169 switch(s[1]) {
170 case 'x':
171 case 'X':
172 {
173 #ifdef Honor_FLT_ROUNDS
174 FPI fpi1 = fpi;
175 fpi1.rounding = Rounding;
176 #else
177 #define fpi1 fpi
178 #endif
179 switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
180 case STRTOG_NoMemory:
181 goto ovfl;
182 case STRTOG_NoNumber:
183 s = s00;
184 sign = 0;
185 case STRTOG_Zero:
186 break;
187 default:
188 if (bb) {
189 copybits(bits, fpi.nbits, bb);
190 Bfree(bb);
191 }
192 ULtod(((U*)&rv)->L, bits, exp, i);
193 }}
194 goto ret;
195 }
196 }
197 #endif /*}*/
198 nz0 = 1;
199 while(*++s == '0') ;
200 if (!*s)
201 goto ret;
202 }
203 s0 = s;
204 y = z = 0;
205 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
206 if (nd < 9)
207 y = 10*y + c - '0';
208 else if (nd < 16)
209 z = 10*z + c - '0';
210 nd0 = nd;
211 #ifdef USE_LOCALE
212 if (c == *decimalpoint) {
213 for(i = 1; decimalpoint[i]; ++i)
214 if (s[i] != decimalpoint[i])
215 goto dig_done;
216 s += i;
217 c = *s;
218 #else
219 if (c == '.') {
220 c = *++s;
221 #endif
222 decpt = 1;
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 case '+':
264 c = *++s;
265 }
266 if (c >= '0' && c <= '9') {
267 while(c == '0')
268 c = *++s;
269 if (c > '0' && c <= '9') {
270 L = c - '0';
271 s1 = s;
272 while((c = *++s) >= '0' && c <= '9')
273 L = 10*L + c - '0';
274 if (s - s1 > 8 || L > 19999)
275 /* Avoid confusion from exponents
276 * so large that e might overflow.
277 */
278 e = 19999; /* safe for 16 bit ints */
279 else
280 e = (int)L;
281 if (esign)
282 e = -e;
283 }
284 else
285 e = 0;
286 }
287 else
288 s = s00;
289 }
290 if (!nd) {
291 if (!nz && !nz0) {
292 #ifdef INFNAN_CHECK
293 /* Check for Nan and Infinity */
294 ULong bits[2];
295 static FPI fpinan = /* only 52 explicit bits */
296 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
297 if (!decpt)
298 switch(c) {
299 case 'i':
300 case 'I':
301 if (match(&s,"nf")) {
302 --s;
303 if (!match(&s,"inity"))
304 ++s;
305 word0(&rv) = 0x7ff00000;
306 word1(&rv) = 0;
307 goto ret;
308 }
309 break;
310 case 'n':
311 case 'N':
312 if (match(&s, "an")) {
313 #ifndef No_Hex_NaN
314 if (*s == '(' /*)*/
315 && hexnan(&s, &fpinan, bits)
316 == STRTOG_NaNbits) {
317 word0(&rv) = 0x7ff00000 | bits[1];
318 word1(&rv) = bits[0];
319 }
320 else {
321 #endif
322 word0(&rv) = NAN_WORD0;
323 word1(&rv) = NAN_WORD1;
324 #ifndef No_Hex_NaN
325 }
326 #endif
327 goto ret;
328 }
329 }
330 #endif /* INFNAN_CHECK */
331 ret0:
332 s = s00;
333 sign = 0;
334 }
335 goto ret;
336 }
337 e1 = e -= nf;
338
339 /* Now we have nd0 digits, starting at s0, followed by a
340 * decimal point, followed by nd-nd0 digits. The number we're
341 * after is the integer represented by those digits times
342 * 10**e */
343
344 if (!nd0)
345 nd0 = nd;
346 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
347 dval(&rv) = y;
348 if (k > 9) {
349 #ifdef SET_INEXACT
350 if (k > DBL_DIG)
351 oldinexact = get_inexact();
352 #endif
353 dval(&rv) = tens[k - 9] * dval(&rv) + z;
354 }
355 if (nd <= DBL_DIG
356 #ifndef RND_PRODQUOT
357 #ifndef Honor_FLT_ROUNDS
358 && Flt_Rounds == 1
359 #endif
360 #endif
361 ) {
362 if (!e)
363 goto ret;
364 #ifndef ROUND_BIASED_without_Round_Up
365 if (e > 0) {
366 if (e <= Ten_pmax) {
367 #ifdef VAX
368 goto vax_ovfl_check;
369 #else
370 #ifdef Honor_FLT_ROUNDS
371 /* round correctly FLT_ROUNDS = 2 or 3 */
372 if (sign) {
373 rv.d = -rv.d;
374 sign = 0;
375 }
376 #endif
377 /* rv = */ rounded_product(dval(&rv), tens[e]);
378 goto ret;
379 #endif
380 }
381 i = DBL_DIG - nd;
382 if (e <= Ten_pmax + i) {
383 /* A fancier test would sometimes let us do
384 * this for larger i values.
385 */
386 #ifdef Honor_FLT_ROUNDS
387 /* round correctly FLT_ROUNDS = 2 or 3 */
388 if (sign) {
389 rv.d = -rv.d;
390 sign = 0;
391 }
392 #endif
393 e -= i;
394 dval(&rv) *= tens[i];
395 #ifdef VAX
396 /* VAX exponent range is so narrow we must
397 * worry about overflow here...
398 */
399 vax_ovfl_check:
400 word0(&rv) -= P*Exp_msk1;
401 /* rv = */ rounded_product(dval(&rv), tens[e]);
402 if ((word0(&rv) & Exp_mask)
403 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
404 goto ovfl;
405 word0(&rv) += P*Exp_msk1;
406 #else
407 /* rv = */ rounded_product(dval(&rv), tens[e]);
408 #endif
409 goto ret;
410 }
411 }
412 #ifndef Inaccurate_Divide
413 else if (e >= -Ten_pmax) {
414 #ifdef Honor_FLT_ROUNDS
415 /* round correctly FLT_ROUNDS = 2 or 3 */
416 if (sign) {
417 rv.d = -rv.d;
418 sign = 0;
419 }
420 #endif
421 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
422 goto ret;
423 }
424 #endif
425 #endif /* ROUND_BIASED_without_Round_Up */
426 }
427 e1 += nd - k;
428
429 #ifdef IEEE_Arith
430 #ifdef SET_INEXACT
431 inexact = 1;
432 if (k <= DBL_DIG)
433 oldinexact = get_inexact();
434 #endif
435 #ifdef Avoid_Underflow
436 scale = 0;
437 #endif
438 #ifdef Honor_FLT_ROUNDS
439 if (Rounding >= 2) {
440 if (sign)
441 Rounding = Rounding == 2 ? 0 : 2;
442 else
443 if (Rounding != 2)
444 Rounding = 0;
445 }
446 #endif
447 #endif /*IEEE_Arith*/
448
449 /* Get starting approximation = rv * 10**e1 */
450
451 if (e1 > 0) {
452 if ( (i = e1 & 15) !=0)
453 dval(&rv) *= tens[i];
454 if (e1 &= ~15) {
455 if (e1 > DBL_MAX_10_EXP) {
456 ovfl:
457 /* Can't trust HUGE_VAL */
458 #ifdef IEEE_Arith
459 #ifdef Honor_FLT_ROUNDS
460 switch(Rounding) {
461 case 0: /* toward 0 */
462 case 3: /* toward -infinity */
463 word0(&rv) = Big0;
464 word1(&rv) = Big1;
465 break;
466 default:
467 word0(&rv) = Exp_mask;
468 word1(&rv) = 0;
469 }
470 #else /*Honor_FLT_ROUNDS*/
471 word0(&rv) = Exp_mask;
472 word1(&rv) = 0;
473 #endif /*Honor_FLT_ROUNDS*/
474 #ifdef SET_INEXACT
475 /* set overflow bit */
476 dval(&rv0) = 1e300;
477 dval(&rv0) *= dval(&rv0);
478 #endif
479 #else /*IEEE_Arith*/
480 word0(&rv) = Big0;
481 word1(&rv) = Big1;
482 #endif /*IEEE_Arith*/
483 range_err:
484 if (bd0) {
485 Bfree(bb);
486 Bfree(bd);
487 Bfree(bs);
488 Bfree(bd0);
489 Bfree(delta);
490 }
491 #ifndef NO_ERRNO
492 errno = ERANGE;
493 #endif
494 goto ret;
495 }
496 e1 >>= 4;
497 for(j = 0; e1 > 1; j++, e1 >>= 1)
498 if (e1 & 1)
499 dval(&rv) *= bigtens[j];
500 /* The last multiplication could overflow. */
501 word0(&rv) -= P*Exp_msk1;
502 dval(&rv) *= bigtens[j];
503 if ((z = word0(&rv) & Exp_mask)
504 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
505 goto ovfl;
506 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
507 /* set to largest number */
508 /* (Can't trust DBL_MAX) */
509 word0(&rv) = Big0;
510 word1(&rv) = Big1;
511 }
512 else
513 word0(&rv) += P*Exp_msk1;
514 }
515 }
516 else if (e1 < 0) {
517 e1 = -e1;
518 if ( (i = e1 & 15) !=0)
519 dval(&rv) /= tens[i];
520 if (e1 >>= 4) {
521 if (e1 >= 1 << n_bigtens)
522 goto undfl;
523 #ifdef Avoid_Underflow
524 if (e1 & Scale_Bit)
525 scale = 2*P;
526 for(j = 0; e1 > 0; j++, e1 >>= 1)
527 if (e1 & 1)
528 dval(&rv) *= tinytens[j];
529 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
530 >> Exp_shift)) > 0) {
531 /* scaled rv is denormal; zap j low bits */
532 if (j >= 32) {
533 word1(&rv) = 0;
534 if (j >= 53)
535 word0(&rv) = (P+2)*Exp_msk1;
536 else
537 word0(&rv) &= 0xffffffff << (j-32);
538 }
539 else
540 word1(&rv) &= 0xffffffff << j;
541 }
542 #else
543 for(j = 0; e1 > 1; j++, e1 >>= 1)
544 if (e1 & 1)
545 dval(&rv) *= tinytens[j];
546 /* The last multiplication could underflow. */
547 dval(&rv0) = dval(&rv);
548 dval(&rv) *= tinytens[j];
549 if (!dval(&rv)) {
550 dval(&rv) = 2.*dval(&rv0);
551 dval(&rv) *= tinytens[j];
552 #endif
553 if (!dval(&rv)) {
554 undfl:
555 dval(&rv) = 0.;
556 goto range_err;
557 }
558 #ifndef Avoid_Underflow
559 word0(&rv) = Tiny0;
560 word1(&rv) = Tiny1;
561 /* The refinement below will clean
562 * this approximation up.
563 */
564 }
565 #endif
566 }
567 }
568
569 /* Now the hard part -- adjusting rv to the correct value.*/
570
571 /* Put digits into bd: true value = bd * 10^e */
572
573 bd0 = s2b(s0, nd0, nd, y, dplen);
574 if (bd0 == NULL)
575 goto ovfl;
576
577 for(;;) {
578 bd = Balloc(bd0->k);
579 if (bd == NULL)
580 goto ovfl;
581 Bcopy(bd, bd0);
582 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
583 if (bb == NULL)
584 goto ovfl;
585 bs = i2b(1);
586 if (bs == NULL)
587 goto ovfl;
588
589 if (e >= 0) {
590 bb2 = bb5 = 0;
591 bd2 = bd5 = e;
592 }
593 else {
594 bb2 = bb5 = -e;
595 bd2 = bd5 = 0;
596 }
597 if (bbe >= 0)
598 bb2 += bbe;
599 else
600 bd2 -= bbe;
601 bs2 = bb2;
602 #ifdef Honor_FLT_ROUNDS
603 if (Rounding != 1)
604 bs2++;
605 #endif
606 #ifdef Avoid_Underflow
607 Lsb = LSB;
608 Lsb1 = 0;
609 j = bbe - scale;
610 i = j + bbbits - 1; /* logb(rv) */
611 j = P + 1 - bbbits;
612 if (i < Emin) { /* denormal */
613 i = Emin - i;
614 j -= i;
615 if (i < 32)
616 Lsb <<= i;
617 else
618 Lsb1 = Lsb << (i-32);
619 }
620 #else /*Avoid_Underflow*/
621 #ifdef Sudden_Underflow
622 #ifdef IBM
623 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
624 #else
625 j = P + 1 - bbbits;
626 #endif
627 #else /*Sudden_Underflow*/
628 j = bbe;
629 i = j + bbbits - 1; /* logb(&rv) */
630 if (i < Emin) /* denormal */
631 j += P - Emin;
632 else
633 j = P + 1 - bbbits;
634 #endif /*Sudden_Underflow*/
635 #endif /*Avoid_Underflow*/
636 bb2 += j;
637 bd2 += j;
638 #ifdef Avoid_Underflow
639 bd2 += scale;
640 #endif
641 i = bb2 < bd2 ? bb2 : bd2;
642 if (i > bs2)
643 i = bs2;
644 if (i > 0) {
645 bb2 -= i;
646 bd2 -= i;
647 bs2 -= i;
648 }
649 if (bb5 > 0) {
650 bs = pow5mult(bs, bb5);
651 if (bs == NULL)
652 goto ovfl;
653 bb1 = mult(bs, bb);
654 if (bb1 == NULL)
655 goto ovfl;
656 Bfree(bb);
657 bb = bb1;
658 }
659 if (bb2 > 0) {
660 bb = lshift(bb, bb2);
661 if (bb == NULL)
662 goto ovfl;
663 }
664 if (bd5 > 0) {
665 bd = pow5mult(bd, bd5);
666 if (bd == NULL)
667 goto ovfl;
668 }
669 if (bd2 > 0) {
670 bd = lshift(bd, bd2);
671 if (bd == NULL)
672 goto ovfl;
673 }
674 if (bs2 > 0) {
675 bs = lshift(bs, bs2);
676 if (bs == NULL)
677 goto ovfl;
678 }
679 delta = diff(bb, bd);
680 if (delta == NULL)
681 goto ovfl;
682 dsign = delta->sign;
683 delta->sign = 0;
684 i = cmp(delta, bs);
685 #ifdef Honor_FLT_ROUNDS
686 if (Rounding != 1) {
687 if (i < 0) {
688 /* Error is less than an ulp */
689 if (!delta->x[0] && delta->wds <= 1) {
690 /* exact */
691 #ifdef SET_INEXACT
692 inexact = 0;
693 #endif
694 break;
695 }
696 if (Rounding) {
697 if (dsign) {
698 dval(&adj) = 1.;
699 goto apply_adj;
700 }
701 }
702 else if (!dsign) {
703 dval(&adj) = -1.;
704 if (!word1(&rv)
705 && !(word0(&rv) & Frac_mask)) {
706 y = word0(&rv) & Exp_mask;
707 #ifdef Avoid_Underflow
708 if (!scale || y > 2*P*Exp_msk1)
709 #else
710 if (y)
711 #endif
712 {
713 delta = lshift(delta,Log2P);
714 if (delta == NULL)
715 goto ovfl;
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 (delta == NULL)
809 goto ovfl;
810 if (cmp(delta, bs) > 0)
811 goto drop_down;
812 break;
813 }
814 if (i == 0) {
815 /* exactly half-way between */
816 if (dsign) {
817 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
818 && word1(&rv) == (
819 #ifdef Avoid_Underflow
820 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
821 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
822 #endif
823 0xffffffff)) {
824 /*boundary case -- increment exponent*/
825 if (word0(&rv) == Big0 && word1(&rv) == Big1)
826 goto ovfl;
827 word0(&rv) = (word0(&rv) & Exp_mask)
828 + Exp_msk1
829 #ifdef IBM
830 | Exp_msk1 >> 4
831 #endif
832 ;
833 word1(&rv) = 0;
834 #ifdef Avoid_Underflow
835 dsign = 0;
836 #endif
837 break;
838 }
839 }
840 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
841 drop_down:
842 /* boundary case -- decrement exponent */
843 #ifdef Sudden_Underflow /*{{*/
844 L = word0(&rv) & Exp_mask;
845 #ifdef IBM
846 if (L < Exp_msk1)
847 #else
848 #ifdef Avoid_Underflow
849 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
850 #else
851 if (L <= Exp_msk1)
852 #endif /*Avoid_Underflow*/
853 #endif /*IBM*/
854 goto undfl;
855 L -= Exp_msk1;
856 #else /*Sudden_Underflow}{*/
857 #ifdef Avoid_Underflow
858 if (scale) {
859 L = word0(&rv) & Exp_mask;
860 if (L <= (2*P+1)*Exp_msk1) {
861 if (L > (P+2)*Exp_msk1)
862 /* round even ==> */
863 /* accept rv */
864 break;
865 /* rv = smallest denormal */
866 goto undfl;
867 }
868 }
869 #endif /*Avoid_Underflow*/
870 L = (word0(&rv) & Exp_mask) - Exp_msk1;
871 #endif /*Sudden_Underflow}}*/
872 word0(&rv) = L | Bndry_mask1;
873 word1(&rv) = 0xffffffff;
874 #ifdef IBM
875 goto cont;
876 #else
877 break;
878 #endif
879 }
880 #ifndef ROUND_BIASED
881 #ifdef Avoid_Underflow
882 if (Lsb1) {
883 if (!(word0(&rv) & Lsb1))
884 break;
885 }
886 else if (!(word1(&rv) & Lsb))
887 break;
888 #else
889 if (!(word1(&rv) & LSB))
890 break;
891 #endif
892 #endif
893 if (dsign)
894 #ifdef Avoid_Underflow
895 dval(&rv) += sulp(&rv, scale);
896 #else
897 dval(&rv) += ulp(&rv);
898 #endif
899 #ifndef ROUND_BIASED
900 else {
901 #ifdef Avoid_Underflow
902 dval(&rv) -= sulp(&rv, scale);
903 #else
904 dval(&rv) -= ulp(&rv);
905 #endif
906 #ifndef Sudden_Underflow
907 if (!dval(&rv))
908 goto undfl;
909 #endif
910 }
911 #ifdef Avoid_Underflow
912 dsign = 1 - dsign;
913 #endif
914 #endif
915 break;
916 }
917 if ((aadj = ratio(delta, bs)) <= 2.) {
918 if (dsign)
919 aadj = dval(&aadj1) = 1.;
920 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
921 #ifndef Sudden_Underflow
922 if (word1(&rv) == Tiny1 && !word0(&rv))
923 goto undfl;
924 #endif
925 aadj = 1.;
926 dval(&aadj1) = -1.;
927 }
928 else {
929 /* special case -- power of FLT_RADIX to be */
930 /* rounded down... */
931
932 if (aadj < 2./FLT_RADIX)
933 aadj = 1./FLT_RADIX;
934 else
935 aadj *= 0.5;
936 dval(&aadj1) = -aadj;
937 }
938 }
939 else {
940 aadj *= 0.5;
941 dval(&aadj1) = dsign ? aadj : -aadj;
942 #ifdef Check_FLT_ROUNDS
943 switch(Rounding) {
944 case 2: /* towards +infinity */
945 dval(&aadj1) -= 0.5;
946 break;
947 case 0: /* towards 0 */
948 case 3: /* towards -infinity */
949 dval(&aadj1) += 0.5;
950 }
951 #else
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 = (char *)s;
1100 return sign ? -dval(&rv) : dval(&rv);
1101 }
1102 DEF_STRONG(strtod);
1103