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