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