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


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