]> git.saurik.com Git - apple/libc.git/blob - gdtoa/FreeBSD/gdtoa-strtod.c
projects / apple / libc.git / blob
? search:


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