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