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