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 if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
118 strlcpy(decimalpoint_cache, s0, strlen(s0) + 1);
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_NoMemory:
184 goto ovfl;
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) = 0x7ff00000 | 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 if (nd <= DBL_DIG
359 #ifndef RND_PRODQUOT
360 #ifndef Honor_FLT_ROUNDS
361 && Flt_Rounds == 1
362 #endif
363 #endif
364 ) {
365 if (!e)
366 goto ret;
367 #ifndef ROUND_BIASED_without_Round_Up
368 if (e > 0) {
369 if (e <= Ten_pmax) {
370 #ifdef VAX
371 goto vax_ovfl_check;
372 #else
373 #ifdef Honor_FLT_ROUNDS
374 /* round correctly FLT_ROUNDS = 2 or 3 */
375 if (sign) {
376 rv.d = -rv.d;
377 sign = 0;
378 }
379 #endif
380 /* rv = */ rounded_product(dval(&rv), tens[e]);
381 goto ret;
382 #endif
383 }
384 i = DBL_DIG - nd;
385 if (e <= Ten_pmax + i) {
386 /* A fancier test would sometimes let us do
387 * this for larger i values.
388 */
389 #ifdef Honor_FLT_ROUNDS
390 /* round correctly FLT_ROUNDS = 2 or 3 */
391 if (sign) {
392 rv.d = -rv.d;
393 sign = 0;
394 }
395 #endif
396 e -= i;
397 dval(&rv) *= tens[i];
398 #ifdef VAX
399 /* VAX exponent range is so narrow we must
400 * worry about overflow here...
401 */
402 vax_ovfl_check:
403 word0(&rv) -= P*Exp_msk1;
404 /* rv = */ rounded_product(dval(&rv), tens[e]);
405 if ((word0(&rv) & Exp_mask)
406 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
407 goto ovfl;
408 word0(&rv) += P*Exp_msk1;
409 #else
410 /* rv = */ rounded_product(dval(&rv), tens[e]);
411 #endif
412 goto ret;
413 }
414 }
415 #ifndef Inaccurate_Divide
416 else if (e >= -Ten_pmax) {
417 #ifdef Honor_FLT_ROUNDS
418 /* round correctly FLT_ROUNDS = 2 or 3 */
419 if (sign) {
420 rv.d = -rv.d;
421 sign = 0;
422 }
423 #endif
424 /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
425 goto ret;
426 }
427 #endif
428 #endif /* ROUND_BIASED_without_Round_Up */
429 }
430 e1 += nd - k;
431
432 #ifdef IEEE_Arith
433 #ifdef SET_INEXACT
434 inexact = 1;
435 if (k <= DBL_DIG)
436 oldinexact = get_inexact();
437 #endif
438 #ifdef Avoid_Underflow
439 scale = 0;
440 #endif
441 #ifdef Honor_FLT_ROUNDS
442 if (Rounding >= 2) {
443 if (sign)
444 Rounding = Rounding == 2 ? 0 : 2;
445 else
446 if (Rounding != 2)
447 Rounding = 0;
448 }
449 #endif
450 #endif /*IEEE_Arith*/
451
452 /* Get starting approximation = rv * 10**e1 */
453
454 if (e1 > 0) {
455 if ( (i = e1 & 15) !=0)
456 dval(&rv) *= tens[i];
457 if (e1 &= ~15) {
458 if (e1 > DBL_MAX_10_EXP) {
459 ovfl:
460 /* Can't trust HUGE_VAL */
461 #ifdef IEEE_Arith
462 #ifdef Honor_FLT_ROUNDS
463 switch(Rounding) {
464 case 0: /* toward 0 */
465 case 3: /* toward -infinity */
466 word0(&rv) = Big0;
467 word1(&rv) = Big1;
468 break;
469 default:
470 word0(&rv) = Exp_mask;
471 word1(&rv) = 0;
472 }
473 #else /*Honor_FLT_ROUNDS*/
474 word0(&rv) = Exp_mask;
475 word1(&rv) = 0;
476 #endif /*Honor_FLT_ROUNDS*/
477 #ifdef SET_INEXACT
478 /* set overflow bit */
479 dval(&rv0) = 1e300;
480 dval(&rv0) *= dval(&rv0);
481 #endif
482 #else /*IEEE_Arith*/
483 word0(&rv) = Big0;
484 word1(&rv) = Big1;
485 #endif /*IEEE_Arith*/
486 range_err:
487 if (bd0) {
488 Bfree(bb);
489 Bfree(bd);
490 Bfree(bs);
491 Bfree(bd0);
492 Bfree(delta);
493 }
494 #ifndef NO_ERRNO
495 errno = ERANGE;
496 #endif
497 goto ret;
498 }
499 e1 >>= 4;
500 for(j = 0; e1 > 1; j++, e1 >>= 1)
501 if (e1 & 1)
502 dval(&rv) *= bigtens[j];
503 /* The last multiplication could overflow. */
504 word0(&rv) -= P*Exp_msk1;
505 dval(&rv) *= bigtens[j];
506 if ((z = word0(&rv) & Exp_mask)
507 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
508 goto ovfl;
509 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
510 /* set to largest number */
511 /* (Can't trust DBL_MAX) */
512 word0(&rv) = Big0;
513 word1(&rv) = Big1;
514 }
515 else
516 word0(&rv) += P*Exp_msk1;
517 }
518 }
519 else if (e1 < 0) {
520 e1 = -e1;
521 if ( (i = e1 & 15) !=0)
522 dval(&rv) /= tens[i];
523 if (e1 >>= 4) {
524 if (e1 >= 1 << n_bigtens)
525 goto undfl;
526 #ifdef Avoid_Underflow
527 if (e1 & Scale_Bit)
528 scale = 2*P;
529 for(j = 0; e1 > 0; j++, e1 >>= 1)
530 if (e1 & 1)
531 dval(&rv) *= tinytens[j];
532 if (scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
533 >> Exp_shift)) > 0) {
534 /* scaled rv is denormal; zap j low bits */
535 if (j >= 32) {
536 word1(&rv) = 0;
537 if (j >= 53)
538 word0(&rv) = (P+2)*Exp_msk1;
539 else
540 word0(&rv) &= 0xffffffff << (j-32);
541 }
542 else
543 word1(&rv) &= 0xffffffff << j;
544 }
545 #else
546 for(j = 0; e1 > 1; j++, e1 >>= 1)
547 if (e1 & 1)
548 dval(&rv) *= tinytens[j];
549 /* The last multiplication could underflow. */
550 dval(&rv0) = dval(&rv);
551 dval(&rv) *= tinytens[j];
552 if (!dval(&rv)) {
553 dval(&rv) = 2.*dval(&rv0);
554 dval(&rv) *= tinytens[j];
555 #endif
556 if (!dval(&rv)) {
557 undfl:
558 dval(&rv) = 0.;
559 goto range_err;
560 }
561 #ifndef Avoid_Underflow
562 word0(&rv) = Tiny0;
563 word1(&rv) = Tiny1;
564 /* The refinement below will clean
565 * this approximation up.
566 */
567 }
568 #endif
569 }
570 }
571
572 /* Now the hard part -- adjusting rv to the correct value.*/
573
574 /* Put digits into bd: true value = bd * 10^e */
575
576 bd0 = s2b(s0, nd0, nd, y, dplen);
577 if (bd0 == NULL)
578 goto ovfl;
579
580 for(;;) {
581 bd = Balloc(bd0->k);
582 if (bd == NULL)
583 goto ovfl;
584 Bcopy(bd, bd0);
585 bb = d2b(dval(&rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
586 if (bb == NULL)
587 goto ovfl;
588 bs = i2b(1);
589 if (bs == NULL)
590 goto ovfl;
591
592 if (e >= 0) {
593 bb2 = bb5 = 0;
594 bd2 = bd5 = e;
595 }
596 else {
597 bb2 = bb5 = -e;
598 bd2 = bd5 = 0;
599 }
600 if (bbe >= 0)
601 bb2 += bbe;
602 else
603 bd2 -= bbe;
604 bs2 = bb2;
605 #ifdef Honor_FLT_ROUNDS
606 if (Rounding != 1)
607 bs2++;
608 #endif
609 #ifdef Avoid_Underflow
610 Lsb = LSB;
611 Lsb1 = 0;
612 j = bbe - scale;
613 i = j + bbbits - 1; /* logb(rv) */
614 j = P + 1 - bbbits;
615 if (i < Emin) { /* denormal */
616 i = Emin - i;
617 j -= i;
618 if (i < 32)
619 Lsb <<= i;
620 else
621 Lsb1 = Lsb << (i-32);
622 }
623 #else /*Avoid_Underflow*/
624 #ifdef Sudden_Underflow
625 #ifdef IBM
626 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
627 #else
628 j = P + 1 - bbbits;
629 #endif
630 #else /*Sudden_Underflow*/
631 j = bbe;
632 i = j + bbbits - 1; /* logb(&rv) */
633 if (i < Emin) /* denormal */
634 j += P - Emin;
635 else
636 j = P + 1 - bbbits;
637 #endif /*Sudden_Underflow*/
638 #endif /*Avoid_Underflow*/
639 bb2 += j;
640 bd2 += j;
641 #ifdef Avoid_Underflow
642 bd2 += scale;
643 #endif
644 i = bb2 < bd2 ? bb2 : bd2;
645 if (i > bs2)
646 i = bs2;
647 if (i > 0) {
648 bb2 -= i;
649 bd2 -= i;
650 bs2 -= i;
651 }
652 if (bb5 > 0) {
653 bs = pow5mult(bs, bb5);
654 if (bs == NULL)
655 goto ovfl;
656 bb1 = mult(bs, bb);
657 if (bb1 == NULL)
658 goto ovfl;
659 Bfree(bb);
660 bb = bb1;
661 }
662 if (bb2 > 0) {
663 bb = lshift(bb, bb2);
664 if (bb == NULL)
665 goto ovfl;
666 }
667 if (bd5 > 0) {
668 bd = pow5mult(bd, bd5);
669 if (bd == NULL)
670 goto ovfl;
671 }
672 if (bd2 > 0) {
673 bd = lshift(bd, bd2);
674 if (bd == NULL)
675 goto ovfl;
676 }
677 if (bs2 > 0) {
678 bs = lshift(bs, bs2);
679 if (bs == NULL)
680 goto ovfl;
681 }
682 delta = diff(bb, bd);
683 if (delta == NULL)
684 goto ovfl;
685 dsign = delta->sign;
686 delta->sign = 0;
687 i = cmp(delta, bs);
688 #ifdef Honor_FLT_ROUNDS
689 if (Rounding != 1) {
690 if (i < 0) {
691 /* Error is less than an ulp */
692 if (!delta->x[0] && delta->wds <= 1) {
693 /* exact */
694 #ifdef SET_INEXACT
695 inexact = 0;
696 #endif
697 break;
698 }
699 if (Rounding) {
700 if (dsign) {
701 dval(&adj) = 1.;
702 goto apply_adj;
703 }
704 }
705 else if (!dsign) {
706 dval(&adj) = -1.;
707 if (!word1(&rv)
708 && !(word0(&rv) & Frac_mask)) {
709 y = word0(&rv) & Exp_mask;
710 #ifdef Avoid_Underflow
711 if (!scale || y > 2*P*Exp_msk1)
712 #else
713 if (y)
714 #endif
715 {
716 delta = lshift(delta,Log2P);
717 if (delta == NULL)
718 goto ovfl;
719 if (cmp(delta, bs) <= 0)
720 dval(&adj) = -0.5;
721 }
722 }
723 apply_adj:
724 #ifdef Avoid_Underflow
725 if (scale && (y = word0(&rv) & Exp_mask)
726 <= 2*P*Exp_msk1)
727 word0(&adj) += (2*P+1)*Exp_msk1 - y;
728 #else
729 #ifdef Sudden_Underflow
730 if ((word0(&rv) & Exp_mask) <=
731 P*Exp_msk1) {
732 word0(&rv) += P*Exp_msk1;
733 dval(&rv) += adj*ulp(&rv);
734 word0(&rv) -= P*Exp_msk1;
735 }
736 else
737 #endif /*Sudden_Underflow*/
738 #endif /*Avoid_Underflow*/
739 dval(&rv) += adj.d*ulp(&rv);
740 }
741 break;
742 }
743 dval(&adj) = ratio(delta, bs);
744 if (adj.d < 1.)
745 dval(&adj) = 1.;
746 if (adj.d <= 0x7ffffffe) {
747 /* dval(&adj) = Rounding ? ceil(&adj) : floor(&adj); */
748 y = adj.d;
749 if (y != adj.d) {
750 if (!((Rounding>>1) ^ dsign))
751 y++;
752 dval(&adj) = y;
753 }
754 }
755 #ifdef Avoid_Underflow
756 if (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
757 word0(&adj) += (2*P+1)*Exp_msk1 - y;
758 #else
759 #ifdef Sudden_Underflow
760 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
761 word0(&rv) += P*Exp_msk1;
762 dval(&adj) *= ulp(&rv);
763 if (dsign)
764 dval(&rv) += adj;
765 else
766 dval(&rv) -= adj;
767 word0(&rv) -= P*Exp_msk1;
768 goto cont;
769 }
770 #endif /*Sudden_Underflow*/
771 #endif /*Avoid_Underflow*/
772 dval(&adj) *= ulp(&rv);
773 if (dsign) {
774 if (word0(&rv) == Big0 && word1(&rv) == Big1)
775 goto ovfl;
776 dval(&rv) += adj.d;
777 }
778 else
779 dval(&rv) -= adj.d;
780 goto cont;
781 }
782 #endif /*Honor_FLT_ROUNDS*/
783
784 if (i < 0) {
785 /* Error is less than half an ulp -- check for
786 * special case of mantissa a power of two.
787 */
788 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask
789 #ifdef IEEE_Arith
790 #ifdef Avoid_Underflow
791 || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
792 #else
793 || (word0(&rv) & Exp_mask) <= Exp_msk1
794 #endif
795 #endif
796 ) {
797 #ifdef SET_INEXACT
798 if (!delta->x[0] && delta->wds <= 1)
799 inexact = 0;
800 #endif
801 break;
802 }
803 if (!delta->x[0] && delta->wds <= 1) {
804 /* exact result */
805 #ifdef SET_INEXACT
806 inexact = 0;
807 #endif
808 break;
809 }
810 delta = lshift(delta,Log2P);
811 if (delta == NULL)
812 goto ovfl;
813 if (cmp(delta, bs) > 0)
814 goto drop_down;
815 break;
816 }
817 if (i == 0) {
818 /* exactly half-way between */
819 if (dsign) {
820 if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
821 && word1(&rv) == (
822 #ifdef Avoid_Underflow
823 (scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
824 ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
825 #endif
826 0xffffffff)) {
827 /*boundary case -- increment exponent*/
828 if (word0(&rv) == Big0 && word1(&rv) == Big1)
829 goto ovfl;
830 word0(&rv) = (word0(&rv) & Exp_mask)
831 + Exp_msk1
832 #ifdef IBM
833 | Exp_msk1 >> 4
834 #endif
835 ;
836 word1(&rv) = 0;
837 #ifdef Avoid_Underflow
838 dsign = 0;
839 #endif
840 break;
841 }
842 }
843 else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
844 drop_down:
845 /* boundary case -- decrement exponent */
846 #ifdef Sudden_Underflow /*{{*/
847 L = word0(&rv) & Exp_mask;
848 #ifdef IBM
849 if (L < Exp_msk1)
850 #else
851 #ifdef Avoid_Underflow
852 if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
853 #else
854 if (L <= Exp_msk1)
855 #endif /*Avoid_Underflow*/
856 #endif /*IBM*/
857 goto undfl;
858 L -= Exp_msk1;
859 #else /*Sudden_Underflow}{*/
860 #ifdef Avoid_Underflow
861 if (scale) {
862 L = word0(&rv) & Exp_mask;
863 if (L <= (2*P+1)*Exp_msk1) {
864 if (L > (P+2)*Exp_msk1)
865 /* round even ==> */
866 /* accept rv */
867 break;
868 /* rv = smallest denormal */
869 goto undfl;
870 }
871 }
872 #endif /*Avoid_Underflow*/
873 L = (word0(&rv) & Exp_mask) - Exp_msk1;
874 #endif /*Sudden_Underflow}}*/
875 word0(&rv) = L | Bndry_mask1;
876 word1(&rv) = 0xffffffff;
877 #ifdef IBM
878 goto cont;
879 #else
880 break;
881 #endif
882 }
883 #ifndef ROUND_BIASED
884 #ifdef Avoid_Underflow
885 if (Lsb1) {
886 if (!(word0(&rv) & Lsb1))
887 break;
888 }
889 else if (!(word1(&rv) & Lsb))
890 break;
891 #else
892 if (!(word1(&rv) & LSB))
893 break;
894 #endif
895 #endif
896 if (dsign)
897 #ifdef Avoid_Underflow
898 dval(&rv) += sulp(&rv, scale);
899 #else
900 dval(&rv) += ulp(&rv);
901 #endif
902 #ifndef ROUND_BIASED
903 else {
904 #ifdef Avoid_Underflow
905 dval(&rv) -= sulp(&rv, scale);
906 #else
907 dval(&rv) -= ulp(&rv);
908 #endif
909 #ifndef Sudden_Underflow
910 if (!dval(&rv))
911 goto undfl;
912 #endif
913 }
914 #ifdef Avoid_Underflow
915 dsign = 1 - dsign;
916 #endif
917 #endif
918 break;
919 }
920 if ((aadj = ratio(delta, bs)) <= 2.) {
921 if (dsign)
922 aadj = dval(&aadj1) = 1.;
923 else if (word1(&rv) || word0(&rv) & Bndry_mask) {
924 #ifndef Sudden_Underflow
925 if (word1(&rv) == Tiny1 && !word0(&rv))
926 goto undfl;
927 #endif
928 aadj = 1.;
929 dval(&aadj1) = -1.;
930 }
931 else {
932 /* special case -- power of FLT_RADIX to be */
933 /* rounded down... */
934
935 if (aadj < 2./FLT_RADIX)
936 aadj = 1./FLT_RADIX;
937 else
938 aadj *= 0.5;
939 dval(&aadj1) = -aadj;
940 }
941 }
942 else {
943 aadj *= 0.5;
944 dval(&aadj1) = dsign ? aadj : -aadj;
945 #ifdef Check_FLT_ROUNDS
946 switch(Rounding) {
947 case 2: /* towards +infinity */
948 dval(&aadj1) -= 0.5;
949 break;
950 case 0: /* towards 0 */
951 case 3: /* towards -infinity */
952 dval(&aadj1) += 0.5;
953 }
954 #else
955 if (Flt_Rounds == 0)
956 dval(&aadj1) += 0.5;
957 #endif /*Check_FLT_ROUNDS*/
958 }
959 y = word0(&rv) & Exp_mask;
960
961 /* Check for overflow */
962
963 if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
964 dval(&rv0) = dval(&rv);
965 word0(&rv) -= P*Exp_msk1;
966 dval(&adj) = dval(&aadj1) * ulp(&rv);
967 dval(&rv) += dval(&adj);
968 if ((word0(&rv) & Exp_mask) >=
969 Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
970 if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
971 goto ovfl;
972 word0(&rv) = Big0;
973 word1(&rv) = Big1;
974 goto cont;
975 }
976 else
977 word0(&rv) += P*Exp_msk1;
978 }
979 else {
980 #ifdef Avoid_Underflow
981 if (scale && y <= 2*P*Exp_msk1) {
982 if (aadj <= 0x7fffffff) {
983 if ((z = aadj) <= 0)
984 z = 1;
985 aadj = z;
986 dval(&aadj1) = dsign ? aadj : -aadj;
987 }
988 word0(&aadj1) += (2*P+1)*Exp_msk1 - y;
989 }
990 dval(&adj) = dval(&aadj1) * ulp(&rv);
991 dval(&rv) += dval(&adj);
992 #else
993 #ifdef Sudden_Underflow
994 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
995 dval(&rv0) = dval(&rv);
996 word0(&rv) += P*Exp_msk1;
997 dval(&adj) = dval(&aadj1) * ulp(&rv);
998 dval(&rv) += dval(&adj);
999 #ifdef IBM
1000 if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
1001 #else
1002 if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
1003 #endif
1004 {
1005 if (word0(&rv0) == Tiny0
1006 && word1(&rv0) == Tiny1)
1007 goto undfl;
1008 word0(&rv) = Tiny0;
1009 word1(&rv) = Tiny1;
1010 goto cont;
1011 }
1012 else
1013 word0(&rv) -= P*Exp_msk1;
1014 }
1015 else {
1016 dval(&adj) = dval(&aadj1) * ulp(&rv);
1017 dval(&rv) += dval(&adj);
1018 }
1019 #else /*Sudden_Underflow*/
1020 /* Compute dval(&adj) so that the IEEE rounding rules will
1021 * correctly round rv + dval(&adj) in some half-way cases.
1022 * If rv * ulp(&rv) is denormalized (i.e.,
1023 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1024 * trouble from bits lost to denormalization;
1025 * example: 1.2e-307 .
1026 */
1027 if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1028 dval(&aadj1) = (double)(int)(aadj + 0.5);
1029 if (!dsign)
1030 dval(&aadj1) = -dval(&aadj1);
1031 }
1032 dval(&adj) = dval(&aadj1) * ulp(&rv);
1033 dval(&rv) += adj;
1034 #endif /*Sudden_Underflow*/
1035 #endif /*Avoid_Underflow*/
1036 }
1037 z = word0(&rv) & Exp_mask;
1038 #ifndef SET_INEXACT
1039 #ifdef Avoid_Underflow
1040 if (!scale)
1041 #endif
1042 if (y == z) {
1043 /* Can we stop now? */
1044 L = (Long)aadj;
1045 aadj -= L;
1046 /* The tolerances below are conservative. */
1047 if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
1048 if (aadj < .4999999 || aadj > .5000001)
1049 break;
1050 }
1051 else if (aadj < .4999999/FLT_RADIX)
1052 break;
1053 }
1054 #endif
1055 cont:
1056 Bfree(bb);
1057 Bfree(bd);
1058 Bfree(bs);
1059 Bfree(delta);
1060 }
1061 Bfree(bb);
1062 Bfree(bd);
1063 Bfree(bs);
1064 Bfree(bd0);
1065 Bfree(delta);
1066 #ifdef SET_INEXACT
1067 if (inexact) {
1068 if (!oldinexact) {
1069 word0(&rv0) = Exp_1 + (70 << Exp_shift);
1070 word1(&rv0) = 0;
1071 dval(&rv0) += 1.;
1072 }
1073 }
1074 else if (!oldinexact)
1075 clear_inexact();
1076 #endif
1077 #ifdef Avoid_Underflow
1078 if (scale) {
1079 word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
1080 word1(&rv0) = 0;
1081 dval(&rv) *= dval(&rv0);
1082 #ifndef NO_ERRNO
1083 /* try to avoid the bug of testing an 8087 register value */
1084 #ifdef IEEE_Arith
1085 if (!(word0(&rv) & Exp_mask))
1086 #else
1087 if (word0(&rv) == 0 && word1(&rv) == 0)
1088 #endif
1089 errno = ERANGE;
1090 #endif
1091 }
1092 #endif /* Avoid_Underflow */
1093 #ifdef SET_INEXACT
1094 if (inexact && !(word0(&rv) & Exp_mask)) {
1095 /* set underflow bit */
1096 dval(&rv0) = 1e-300;
1097 dval(&rv0) *= dval(&rv0);
1098 }
1099 #endif
1100 ret:
1101 if (se)
1102 *se = (char *)s;
1103 return sign ? -dval(&rv) : dval(&rv);
1104 }
1105
1106