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 #if !defined(NO_FENV_H) && !defined(_MSC_VER)
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, *bb1, *bd, *bd0, *bs, *delta;
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 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
118 strcpy(decimalpoint_cache, s0);
119 s0 = decimalpoint_cache;
120 }
121 dplen = strlen(s0);
122 }
123 decimalpoint = (char*)s0;
124 #endif /*NO_LOCALE_CACHE*/
125 #else /*USE_LOCALE}{*/
126 #define dplen 1
127 #endif /*USE_LOCALE}}*/
128
129 #ifdef Honor_FLT_ROUNDS /*{*/
130 int Rounding;
131 #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
132 Rounding = Flt_Rounds;
133 #else /*}{*/
134 Rounding = 1;
135 switch(fegetround()) {
136 case FE_TOWARDZERO: Rounding = 0; break;
137 case FE_UPWARD: Rounding = 2; break;
138 case FE_DOWNWARD: Rounding = 3;
139 }
140 #endif /*}}*/
141 #endif /*}*/
142
143 sign = nz0 = nz = decpt = 0;
144 dval(&rv) = 0.;
145 for(s = s00;;s++) switch(*s) {
146 case '-':
147 sign = 1;
148 /* no break */
149 case '+':
150 if (*++s)
151 goto break2;
152 /* no break */
153 case 0:
154 goto ret0;
155 case '\t':
156 case '\n':
157 case '\v':
158 case '\f':
159 case '\r':
160 case ' ':
161 continue;
162 default:
163 goto break2;
164 }
165 break2:
166 if (*s == '0') {
167 #ifndef NO_HEX_FP /*{*/
168 {
169 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
170 Long exp;
171 ULong bits[2];
172 switch(s[1]) {
173 case 'x':
174 case 'X':
175 {
176 #ifdef Honor_FLT_ROUNDS
177 FPI fpi1 = fpi;
178 fpi1.rounding = Rounding;
179 #else
180 #define fpi1 fpi
181 #endif
182 switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
183 case STRTOG_NoNumber:
184 s = s00;
185 sign = 0;
186 case STRTOG_Zero:
187 break;
188 default:
189 if (bb) {
190 copybits(bits, fpi.nbits, bb);
191 Bfree(bb);
192 }
193 ULtod(((U*)&rv)->L, bits, exp, i);
194 }}
195 goto ret;
196 }
197 }
198 #endif /*}*/
199 nz0 = 1;
200 while(*++s == '0') ;
201 if (!*s)
202 goto ret;
203 }
204 s0 = s;
205 y = z = 0;
206 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
207 if (nd < 9)
208 y = 10*y + c - '0';
209 else if (nd < 16)
210 z = 10*z + c - '0';
211 nd0 = nd;
212 #ifdef USE_LOCALE
213 if (c == *decimalpoint) {
214 for(i = 1; decimalpoint[i]; ++i)
215 if (s[i] != decimalpoint[i])
216 goto dig_done;
217 s += i;
218 c = *s;
219 #else
220 if (c == '.') {
221 c = *++s;
222 #endif
223 decpt = 1;
224 if (!nd) {
225 for(; c == '0'; c = *++s)
226 nz++;
227 if (c > '0' && c <= '9') {
228 s0 = s;
229 nf += nz;
230 nz = 0;
231 goto have_dig;
232 }
233 goto dig_done;
234 }
235 for(; c >= '0' && c <= '9'; c = *++s) {
236 have_dig:
237 nz++;
238 if (c -= '0') {
239 nf += nz;
240 for(i = 1; i < nz; i++)
241 if (nd++ < 9)
242 y *= 10;
243 else if (nd <= DBL_DIG + 1)
244 z *= 10;
245 if (nd++ < 9)
246 y = 10*y + c;
247 else if (nd <= DBL_DIG + 1)
248 z = 10*z + c;
249 nz = 0;
250 }
251 }
252 }/*}*/
253 dig_done:
254 e = 0;
255 if (c == 'e' || c == 'E') {
256 if (!nd && !nz && !nz0) {
257 goto ret0;
258 }
259 s00 = s;
260 esign = 0;
261 switch(c = *++s) {
262 case '-':
263 esign = 1;
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 >>= 4;
499 for(j = 0; e1 > 1; j++, 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 >>= 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) &= 0xffffffff << (j-32);
540 }
541 else
542 word1(&rv) &= 0xffffffff << j;
543 }
544 #else
545 for(j = 0; e1 > 1; j++, 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
577 for(;;) {
578 bd = Balloc(bd0->k);
579 Bcopy(bd, bd0);
580 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
581 bs = i2b(1);
582
583 if (e >= 0) {
584 bb2 = bb5 = 0;
585 bd2 = bd5 = e;
586 }
587 else {
588 bb2 = bb5 = -e;
589 bd2 = bd5 = 0;
590 }
591 if (bbe >= 0)
592 bb2 += bbe;
593 else
594 bd2 -= bbe;
595 bs2 = bb2;
596 #ifdef Honor_FLT_ROUNDS
597 if (Rounding != 1)
598 bs2++;
599 #endif
600 #ifdef Avoid_Underflow
601 Lsb = LSB;
602 Lsb1 = 0;
603 j = bbe - scale;
604 i = j + bbbits - 1; /* logb(rv) */
605 j = P + 1 - bbbits;
606 if (i < Emin) { /* denormal */
607 i = Emin - i;
608 j -= i;
609 if (i < 32)
610 Lsb <<= i;
611 else
612 Lsb1 = Lsb << (i-32);
613 }
614 #else /*Avoid_Underflow*/
615 #ifdef Sudden_Underflow
616 #ifdef IBM
617 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
618 #else
619 j = P + 1 - bbbits;
620 #endif
621 #else /*Sudden_Underflow*/
622 j = bbe;
623 i = j + bbbits - 1; /* logb(&rv) */
624 if (i < Emin) /* denormal */
625 j += P - Emin;
626 else
627 j = P + 1 - bbbits;
628 #endif /*Sudden_Underflow*/
629 #endif /*Avoid_Underflow*/
630 bb2 += j;
631 bd2 += j;
632 #ifdef Avoid_Underflow
633 bd2 += scale;
634 #endif
635 i = bb2 < bd2 ? bb2 : bd2;
636 if (i > bs2)
637 i = bs2;
638 if (i > 0) {
639 bb2 -= i;
640 bd2 -= i;
641 bs2 -= i;
642 }
643 if (bb5 > 0) {
644 bs = pow5mult(bs, bb5);
645 bb1 = mult(bs, bb);
646 Bfree(bb);
647 bb = bb1;
648 }
649 if (bb2 > 0)
650 bb = lshift(bb, bb2);
651 if (bd5 > 0)
652 bd = pow5mult(bd, bd5);
653 if (bd2 > 0)
654 bd = lshift(bd, bd2);
655 if (bs2 > 0)
656 bs = lshift(bs, bs2);
657 delta = diff(bb, bd);
658 dsign = delta->sign;
659 delta->sign = 0;
660 i = cmp(delta, bs);
661 #ifdef Honor_FLT_ROUNDS
662 if (Rounding != 1) {
663 if (i < 0) {
664 /* Error is less than an ulp */
665 if (!delta->x[0] && delta->wds <= 1) {
666 /* exact */
667 #ifdef SET_INEXACT
668 inexact = 0;
669 #endif
670 break;
671 }
672 if (Rounding) {
673 if (dsign) {
674 dval(&adj) = 1.;
675 goto apply_adj;
676 }
677 }
678 else if (!dsign) {
679 dval(&adj) = -1.;
680 if (!word1(&rv)
681 && !(word0(&rv) & Frac_mask)) {
682 y = word0(&rv) & Exp_mask;
683 #ifdef Avoid_Underflow
684 if (!scale || y > 2*P*Exp_msk1)
685 #else
686 if (y)
687 #endif
688 {
689 delta = lshift(delta,Log2P);
690 if (cmp(delta, bs) <= 0)
691 dval(&adj) = -0.5;
692 }
693 }
694 apply_adj:
695 #ifdef Avoid_Underflow
696 if (scale && (y = word0(&rv) & Exp_mask)
697 <= 2*P*Exp_msk1)
698 word0(&adj) += (2*P+1)*Exp_msk1 - y;
699 #else
700 #ifdef Sudden_Underflow
701 if ((word0(&rv) & Exp_mask) <=
702 P*Exp_msk1) {
703 word0(&rv) += P*Exp_msk1;
704 dval(&rv) += adj*ulp(&rv);
705 word0(&rv) -= P*Exp_msk1;
706 }
707 else
708 #endif /*Sudden_Underflow*/
709 #endif /*Avoid_Underflow*/
710 dval(&rv) += adj.d*ulp(&rv);
711 }
712 break;
713 }
714 dval(&adj) = ratio(delta, bs);
715 if (adj.d < 1.)
716 dval(&adj) = 1.;
717 if (adj.d <= 0x7ffffffe) {
718 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
719 y = adj.d;
720 if (y != adj.d) {
721 if (!((Rounding>>1) ^ dsign))
722 y++;
723 dval(&adj) = y;
724 }
725 }
726 #ifdef Avoid_Underflow
727 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
728 word0(&adj) += (2*P+1)*Exp_msk1 - y;
729 #else
730 #ifdef Sudden_Underflow
731 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
732 word0(&rv) += P*Exp_msk1;
733 dval(&adj) *= ulp(&rv);
734 if (dsign)
735 dval(&rv) += adj;
736 else
737 dval(&rv) -= adj;
738 word0(&rv) -= P*Exp_msk1;
739 goto cont;
740 }
741 #endif /*Sudden_Underflow*/
742 #endif /*Avoid_Underflow*/
743 dval(&adj) *= ulp(&rv);
744 if (dsign) {
745 if (word0(&rv) == Big0 && word1(&rv) == Big1)
746 goto ovfl;
747 dval(&rv) += adj.d;
748 }
749 else
750 dval(&rv) -= adj.d;
751 goto cont;
752 }
753 #endif /*Honor_FLT_ROUNDS*/
754
755 if (i < 0) {
756 /* Error is less than half an ulp -- check for
757 * special case of mantissa a power of two.
758 */
759 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
760 #ifdef IEEE_Arith
761 #ifdef Avoid_Underflow
762 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
763 #else
764 || (word0(&rv) & Exp_mask) <= Exp_msk1
765 #endif
766 #endif
767 ) {
768 #ifdef SET_INEXACT
769 if (!delta->x[0] && delta->wds <= 1)
770 inexact = 0;
771 #endif
772 break;
773 }
774 if (!delta->x[0] && delta->wds <= 1) {
775 /* exact result */
776 #ifdef SET_INEXACT
777 inexact = 0;
778 #endif
779 break;
780 }
781 delta = lshift(delta,Log2P);
782 if (cmp(delta, bs) > 0)
783 goto drop_down;
784 break;
785 }
786 if (i == 0) {
787 /* exactly half-way between */
788 if (dsign) {
789 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
790 && word1(&rv) == (
791 #ifdef Avoid_Underflow
792 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
793 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
794 #endif
795 0xffffffff)) {
796 /*boundary case -- increment exponent*/
797 if (word0(&rv) == Big0 && word1(&rv) == Big1)
798 goto ovfl;
799 word0(&rv) = (word0(&rv) & Exp_mask)
800 + Exp_msk1
801 #ifdef IBM
802 | Exp_msk1 >> 4
803 #endif
804 ;
805 word1(&rv) = 0;
806 #ifdef Avoid_Underflow
807 dsign = 0;
808 #endif
809 break;
810 }
811 }
812 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
813 drop_down:
814 /* boundary case -- decrement exponent */
815 #ifdef Sudden_Underflow /*{{*/
816 L = word0(&rv) & Exp_mask;
817 #ifdef IBM
818 if (L < Exp_msk1)
819 #else
820 #ifdef Avoid_Underflow
821 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
822 #else
823 if (L <= Exp_msk1)
824 #endif /*Avoid_Underflow*/
825 #endif /*IBM*/
826 goto undfl;
827 L -= Exp_msk1;
828 #else /*Sudden_Underflow}{*/
829 #ifdef Avoid_Underflow
830 if (scale) {
831 L = word0(&rv) & Exp_mask;
832 if (L <= (2*P+1)*Exp_msk1) {
833 if (L > (P+2)*Exp_msk1)
834 /* round even ==> */
835 /* accept rv */
836 break;
837 /* rv = smallest denormal */
838 goto undfl;
839 }
840 }
841 #endif /*Avoid_Underflow*/
842 L = (word0(&rv) & Exp_mask) - Exp_msk1;
843 #endif /*Sudden_Underflow}}*/
844 word0(&rv) = L | Bndry_mask1;
845 word1(&rv) = 0xffffffff;
846 #ifdef IBM
847 goto cont;
848 #else
849 break;
850 #endif
851 }
852 #ifndef ROUND_BIASED
853 #ifdef Avoid_Underflow
854 if (Lsb1) {
855 if (!(word0(&rv) & Lsb1))
856 break;
857 }
858 else if (!(word1(&rv) & Lsb))
859 break;
860 #else
861 if (!(word1(&rv) & LSB))
862 break;
863 #endif
864 #endif
865 if (dsign)
866 #ifdef Avoid_Underflow
867 dval(&rv) += sulp(&rv, scale);
868 #else
869 dval(&rv) += ulp(&rv);
870 #endif
871 #ifndef ROUND_BIASED
872 else {
873 #ifdef Avoid_Underflow
874 dval(&rv) -= sulp(&rv, scale);
875 #else
876 dval(&rv) -= ulp(&rv);
877 #endif
878 #ifndef Sudden_Underflow
879 if (!dval(&rv))
880 goto undfl;
881 #endif
882 }
883 #ifdef Avoid_Underflow
884 dsign = 1 - dsign;
885 #endif
886 #endif
887 break;
888 }
889 if ((aadj = ratio(delta, bs)) <= 2.) {
890 if (dsign)
891 aadj = dval(&aadj1) = 1.;
892 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
893 #ifndef Sudden_Underflow
894 if (word1(&rv) == Tiny1 && !word0(&rv))
895 goto undfl;
896 #endif
897 aadj = 1.;
898 dval(&aadj1) = -1.;
899 }
900 else {
901 /* special case -- power of FLT_RADIX to be */
902 /* rounded down... */
903
904 if (aadj < 2./FLT_RADIX)
905 aadj = 1./FLT_RADIX;
906 else
907 aadj *= 0.5;
908 dval(&aadj1) = -aadj;
909 }
910 }
911 else {
912 aadj *= 0.5;
913 dval(&aadj1) = dsign ? aadj : -aadj;
914 #ifdef Check_FLT_ROUNDS
915 switch(Rounding) {
916 case 2: /* towards +infinity */
917 dval(&aadj1) -= 0.5;
918 break;
919 case 0: /* towards 0 */
920 case 3: /* towards -infinity */
921 dval(&aadj1) += 0.5;
922 }
923 #else
924 if (Flt_Rounds == 0)
925 dval(&aadj1) += 0.5;
926 #endif /*Check_FLT_ROUNDS*/
927 }
928 y = word0(&rv) & Exp_mask;
929
930 /* Check for overflow */
931
932 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
933 dval(&rv0) = dval(&rv);
934 word0(&rv) -= P*Exp_msk1;
935 dval(&adj) = dval(&aadj1) * ulp(&rv);
936 dval(&rv) += dval(&adj);
937 if ((word0(&rv) & Exp_mask) >=
938 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
939 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
940 goto ovfl;
941 word0(&rv) = Big0;
942 word1(&rv) = Big1;
943 goto cont;
944 }
945 else
946 word0(&rv) += P*Exp_msk1;
947 }
948 else {
949 #ifdef Avoid_Underflow
950 if (scale && y <= 2*P*Exp_msk1) {
951 if (aadj <= 0x7fffffff) {
952 if ((z = (ULong)aadj) <= 0)
953 z = 1;
954 aadj = z;
955 dval(&aadj1) = dsign ? aadj : -aadj;
956 }
957 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
958 }
959 dval(&adj) = dval(&aadj1) * ulp(&rv);
960 dval(&rv) += dval(&adj);
961 #else
962 #ifdef Sudden_Underflow
963 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
964 dval(&rv0) = dval(&rv);
965 word0(&rv) += P*Exp_msk1;
966 dval(&adj) = dval(&aadj1) * ulp(&rv);
967 dval(&rv) += adj;
968 #ifdef IBM
969 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
970 #else
971 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
972 #endif
973 {
974 if (word0(&rv0) == Tiny0
975 && word1(&rv0) == Tiny1)
976 goto undfl;
977 word0(&rv) = Tiny0;
978 word1(&rv) = Tiny1;
979 goto cont;
980 }
981 else
982 word0(&rv) -= P*Exp_msk1;
983 }
984 else {
985 dval(&adj) = dval(&aadj1) * ulp(&rv);
986 dval(&rv) += adj;
987 }
988 #else /*Sudden_Underflow*/
989 /* Compute dval(&adj) so that the IEEE rounding rules will
990 * correctly round rv + dval(&adj) in some half-way cases.
991 * If rv * ulp(&rv) is denormalized (i.e.,
992 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
993 * trouble from bits lost to denormalization;
994 * example: 1.2e-307 .
995 */
996 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
997 dval(&aadj1) = (double)(int)(aadj + 0.5);
998 if (!dsign)
999 dval(&aadj1) = -dval(&aadj1);
1000 }
1001 dval(&adj) = dval(&aadj1) * ulp(&rv);
1002 dval(&rv) += adj;
1003 #endif /*Sudden_Underflow*/
1004 #endif /*Avoid_Underflow*/
1005 }
1006 z = word0(&rv) & Exp_mask;
1007 #ifndef SET_INEXACT
1008 #ifdef Avoid_Underflow
1009 if (!scale)
1010 #endif
1011 if (y == z) {
1012 /* Can we stop now? */
1013 L = (Long)aadj;
1014 aadj -= L;
1015 /* The tolerances below are conservative. */
1016 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1017 if (aadj < .4999999 || aadj > .5000001)
1018 break;
1019 }
1020 else if (aadj < .4999999/FLT_RADIX)
1021 break;
1022 }
1023 #endif
1024 cont:
1025 Bfree(bb);
1026 Bfree(bd);
1027 Bfree(bs);
1028 Bfree(delta);
1029 }
1030 Bfree(bb);
1031 Bfree(bd);
1032 Bfree(bs);
1033 Bfree(bd0);
1034 Bfree(delta);
1035 #ifdef SET_INEXACT
1036 if (inexact) {
1037 if (!oldinexact) {
1038 word0(&rv0) = Exp_1 + (70 << Exp_shift);
1039 word1(&rv0) = 0;
1040 dval(&rv0) += 1.;
1041 }
1042 }
1043 else if (!oldinexact)
1044 clear_inexact();
1045 #endif
1046 #ifdef Avoid_Underflow
1047 if (scale) {
1048 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1049 word1(&rv0) = 0;
1050 dval(&rv) *= dval(&rv0);
1051 #ifndef NO_ERRNO
1052 /* try to avoid the bug of testing an 8087 register value */
1053 #ifdef IEEE_Arith
1054 if (!(word0(&rv) & Exp_mask))
1055 #else
1056 if (word0(&rv) == 0 && word1(&rv) == 0)
1057 #endif
1058 errno = ERANGE;
1059 #endif
1060 }
1061 #endif /* Avoid_Underflow */
1062 #ifdef SET_INEXACT
1063 if (inexact && !(word0(&rv) & Exp_mask)) {
1064 /* set underflow bit */
1065 dval(&rv0) = 1e-300;
1066 dval(&rv0) *= dval(&rv0);
1067 }
1068 #endif
1069 ret:
1070 if (se)
1071 *se = (char *)s;
1072 return sign ? -dval(&rv) : dval(&rv);
1073 }
1074
1075