?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix mdoc(7)/man(7) mix up.
[netbsd-mini2440.git] / lib / libc / gdtoa / strtod.c
blobaf1c1c7854bb9b2cafa40382dd7d95571bd4d527
1 /* $NetBSD: strtod.c,v 1.4 2006/06/02 19:46:56 mrg Exp $ */
2
3 /****************************************************************
4
5 The author of this software is David M. Gay.
6
7 Copyright (C) 1998-2001 by Lucent Technologies
8 All Rights Reserved
9
10 Permission to use, copy, modify, and distribute this software and
11 its documentation for any purpose and without fee is hereby
12 granted, provided that the above copyright notice appear in all
13 copies and that both that the copyright notice and this
14 permission notice and warranty disclaimer appear in supporting
15 documentation, and that the name of Lucent or any of its entities
16 not be used in advertising or publicity pertaining to
17 distribution of the software without specific, written prior
18 permission.
19
20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
27 THIS SOFTWARE.
28
29 ****************************************************************/
30
31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
32 * with " at " changed at "@" and " dot " changed to "."). */
33
34 #include "gdtoaimp.h"
35 #ifndef NO_FENV_H
36 #include <fenv.h>
37 #endif
38
39 #ifdef USE_LOCALE
40 #include "locale.h"
41 #endif
42
43 #ifdef IEEE_Arith
44 #ifndef NO_IEEE_Scale
45 #define Avoid_Underflow
46 #undef tinytens
47 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
48 /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
49 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
50 9007199254740992.e-256
51 };
52 #endif
53 #endif
54
55 #ifdef Honor_FLT_ROUNDS
56 #define Rounding rounding
57 #undef Check_FLT_ROUNDS
58 #define Check_FLT_ROUNDS
59 #else
60 #define Rounding Flt_Rounds
61 #endif
62
63 #ifndef __HAVE_LONG_DOUBLE
64 __strong_alias(_strtold, strtod)
65 __weak_alias(strtold, _strtold)
66 #endif
67
68 double
69 strtod
70 #ifdef KR_headers
71 (s00, se) CONST char *s00; char **se;
72 #else
73 (CONST char *s00, char **se)
74 #endif
75 {
76 #ifdef Avoid_Underflow
77 int scale;
78 #endif
79 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
80 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
81 CONST char *s, *s0, *s1;
82 double aadj, aadj1, adj, rv, rv0;
83 Long L;
84 ULong y, z;
85 Bigint *bb = NULL, *bb1, *bd0;
86 Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
87 #ifdef SET_INEXACT
88 int inexact, oldinexact;
89 #endif
90 #ifdef Honor_FLT_ROUNDS
91 int rounding;
92 #endif
93
94 sign = nz0 = nz = decpt = 0;
95 dval(rv) = 0.;
96 for(s = s00;;s++) switch(*s) {
97 case '-':
98 sign = 1;
99 /* FALLTHROUGH */
100 case '+':
101 if (*++s)
102 goto break2;
103 /* FALLTHROUGH */
104 case 0:
105 goto ret0;
106 case '\t':
107 case '\n':
108 case '\v':
109 case '\f':
110 case '\r':
111 case ' ':
112 continue;
113 default:
114 goto break2;
115 }
116 break2:
117 if (*s == '0') {
118 #ifndef NO_HEX_FP
119 {
120 static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
121 Long expt;
122 ULong bits[2];
123 switch(s[1]) {
124 case 'x':
125 case 'X':
126 {
127 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
128 FPI fpi1 = fpi;
129 switch(fegetround()) {
130 case FE_TOWARDZERO: fpi1.rounding = 0; break;
131 case FE_UPWARD: fpi1.rounding = 2; break;
132 case FE_DOWNWARD: fpi1.rounding = 3;
133 }
134 #else
135 #define fpi1 fpi
136 #endif
137 switch((i = gethex(&s, &fpi1, &expt, &bb, sign)) & STRTOG_Retmask) {
138 case STRTOG_NoNumber:
139 s = s00;
140 sign = 0;
141 /* FALLTHROUGH */
142 case STRTOG_Zero:
143 break;
144 default:
145 if (bb) {
146 copybits(bits, fpi.nbits, bb);
147 Bfree(bb);
148 }
149 ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
150 }}
151 goto ret;
152 }
153 }
154 #endif
155 nz0 = 1;
156 while(*++s == '0') ;
157 if (!*s)
158 goto ret;
159 }
160 s0 = s;
161 y = z = 0;
162 for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
163 if (nd < 9)
164 y = 10*y + c - '0';
165 else if (nd < 16)
166 z = 10*z + c - '0';
167 nd0 = nd;
168 #ifdef USE_LOCALE
169 if (c == *localeconv()->decimal_point)
170 #else
171 if (c == '.')
172 #endif
173 {
174 decpt = 1;
175 c = *++s;
176 if (!nd) {
177 for(; c == '0'; c = *++s)
178 nz++;
179 if (c > '0' && c <= '9') {
180 s0 = s;
181 nf += nz;
182 nz = 0;
183 goto have_dig;
184 }
185 goto dig_done;
186 }
187 for(; c >= '0' && c <= '9'; c = *++s) {
188 have_dig:
189 nz++;
190 if (c -= '0') {
191 nf += nz;
192 for(i = 1; i < nz; i++)
193 if (nd++ < 9)
194 y *= 10;
195 else if (nd <= DBL_DIG + 1)
196 z *= 10;
197 if (nd++ < 9)
198 y = 10*y + c;
199 else if (nd <= DBL_DIG + 1)
200 z = 10*z + c;
201 nz = 0;
202 }
203 }
204 }
205 dig_done:
206 e = 0;
207 if (c == 'e' || c == 'E') {
208 if (!nd && !nz && !nz0) {
209 goto ret0;
210 }
211 s00 = s;
212 esign = 0;
213 switch(c = *++s) {
214 case '-':
215 esign = 1;
216 /* FALLTHROUGH */
217 case '+':
218 c = *++s;
219 }
220 if (c >= '0' && c <= '9') {
221 while(c == '0')
222 c = *++s;
223 if (c > '0' && c <= '9') {
224 L = c - '0';
225 s1 = s;
226 while((c = *++s) >= '0' && c <= '9')
227 L = 10*L + c - '0';
228 if (s - s1 > 8 || L > 19999)
229 /* Avoid confusion from exponents
230 * so large that e might overflow.
231 */
232 e = 19999; /* safe for 16 bit ints */
233 else
234 e = (int)L;
235 if (esign)
236 e = -e;
237 }
238 else
239 e = 0;
240 }
241 else
242 s = s00;
243 }
244 if (!nd) {
245 if (!nz && !nz0) {
246 #ifdef INFNAN_CHECK
247 /* Check for Nan and Infinity */
248 ULong bits[2];
249 static FPI fpinan = /* only 52 explicit bits */
250 { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
251 if (!decpt)
252 switch(c) {
253 case 'i':
254 case 'I':
255 if (match(&s,"nf")) {
256 --s;
257 if (!match(&s,"inity"))
258 ++s;
259 word0(rv) = 0x7ff00000;
260 word1(rv) = 0;
261 goto ret;
262 }
263 break;
264 case 'n':
265 case 'N':
266 if (match(&s, "an")) {
267 #ifndef No_Hex_NaN
268 if (*s == '(' /*)*/
269 && hexnan(&s, &fpinan, bits)
270 == STRTOG_NaNbits) {
271 word0(rv) = 0x7ff00000 | bits[1];
272 word1(rv) = bits[0];
273 }
274 else {
275 #endif
276 word0(rv) = NAN_WORD0;
277 word1(rv) = NAN_WORD1;
278 #ifndef No_Hex_NaN
279 }
280 #endif
281 goto ret;
282 }
283 }
284 #endif /* INFNAN_CHECK */
285 ret0:
286 s = s00;
287 sign = 0;
288 }
289 goto ret;
290 }
291 e1 = e -= nf;
292
293 /* Now we have nd0 digits, starting at s0, followed by a
294 * decimal point, followed by nd-nd0 digits. The number we're
295 * after is the integer represented by those digits times
296 * 10**e */
297
298 if (!nd0)
299 nd0 = nd;
300 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
301 dval(rv) = y;
302 if (k > 9) {
303 #ifdef SET_INEXACT
304 if (k > DBL_DIG)
305 oldinexact = get_inexact();
306 #endif
307 dval(rv) = tens[k - 9] * dval(rv) + z;
308 }
309 bd0 = 0;
310 if (nd <= DBL_DIG
311 #ifndef RND_PRODQUOT
312 #ifndef Honor_FLT_ROUNDS
313 && Flt_Rounds == 1
314 #endif
315 #endif
316 ) {
317 if (!e)
318 goto ret;
319 if (e > 0) {
320 if (e <= Ten_pmax) {
321 #ifdef VAX
322 goto vax_ovfl_check;
323 #else
324 #ifdef Honor_FLT_ROUNDS
325 /* round correctly FLT_ROUNDS = 2 or 3 */
326 if (sign) {
327 rv = -rv;
328 sign = 0;
329 }
330 #endif
331 /* rv = */ rounded_product(dval(rv), tens[e]);
332 goto ret;
333 #endif
334 }
335 i = DBL_DIG - nd;
336 if (e <= Ten_pmax + i) {
337 /* A fancier test would sometimes let us do
338 * this for larger i values.
339 */
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 e -= i;
348 dval(rv) *= tens[i];
349 #ifdef VAX
350 /* VAX exponent range is so narrow we must
351 * worry about overflow here...
352 */
353 vax_ovfl_check:
354 word0(rv) -= P*Exp_msk1;
355 /* rv = */ rounded_product(dval(rv), tens[e]);
356 if ((word0(rv) & Exp_mask)
357 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
358 goto ovfl;
359 word0(rv) += P*Exp_msk1;
360 #else
361 /* rv = */ rounded_product(dval(rv), tens[e]);
362 #endif
363 goto ret;
364 }
365 }
366 #ifndef Inaccurate_Divide
367 else if (e >= -Ten_pmax) {
368 #ifdef Honor_FLT_ROUNDS
369 /* round correctly FLT_ROUNDS = 2 or 3 */
370 if (sign) {
371 rv = -rv;
372 sign = 0;
373 }
374 #endif
375 /* rv = */ rounded_quotient(dval(rv), tens[-e]);
376 goto ret;
377 }
378 #endif
379 }
380 e1 += nd - k;
381
382 #ifdef IEEE_Arith
383 #ifdef SET_INEXACT
384 inexact = 1;
385 if (k <= DBL_DIG)
386 oldinexact = get_inexact();
387 #endif
388 #ifdef Avoid_Underflow
389 scale = 0;
390 #endif
391 #ifdef Honor_FLT_ROUNDS
392 if ((rounding = Flt_Rounds) >= 2) {
393 if (sign)
394 rounding = rounding == 2 ? 0 : 2;
395 else
396 if (rounding != 2)
397 rounding = 0;
398 }
399 #endif
400 #endif /*IEEE_Arith*/
401
402 /* Get starting approximation = rv * 10**e1 */
403
404 if (e1 > 0) {
405 if ( (i = e1 & 15) !=0)
406 dval(rv) *= tens[i];
407 if (e1 &= ~15) {
408 if (e1 > DBL_MAX_10_EXP) {
409 ovfl:
410 #ifndef NO_ERRNO
411 errno = ERANGE;
412 #endif
413 /* Can't trust HUGE_VAL */
414 #ifdef IEEE_Arith
415 #ifdef Honor_FLT_ROUNDS
416 switch(rounding) {
417 case 0: /* toward 0 */
418 case 3: /* toward -infinity */
419 word0(rv) = Big0;
420 word1(rv) = Big1;
421 break;
422 default:
423 word0(rv) = Exp_mask;
424 word1(rv) = 0;
425 }
426 #else /*Honor_FLT_ROUNDS*/
427 word0(rv) = Exp_mask;
428 word1(rv) = 0;
429 #endif /*Honor_FLT_ROUNDS*/
430 #ifdef SET_INEXACT
431 /* set overflow bit */
432 dval(rv0) = 1e300;
433 dval(rv0) *= dval(rv0);
434 #endif
435 #else /*IEEE_Arith*/
436 word0(rv) = Big0;
437 word1(rv) = Big1;
438 #endif /*IEEE_Arith*/
439 if (bd0)
440 goto retfree;
441 goto ret;
442 }
443 e1 = (unsigned int)e1 >> 4;
444 for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
445 if (e1 & 1)
446 dval(rv) *= bigtens[j];
447 /* The last multiplication could overflow. */
448 word0(rv) -= P*Exp_msk1;
449 dval(rv) *= bigtens[j];
450 if ((z = word0(rv) & Exp_mask)
451 > Exp_msk1*(DBL_MAX_EXP+Bias-P))
452 goto ovfl;
453 if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
454 /* set to largest number */
455 /* (Can't trust DBL_MAX) */
456 word0(rv) = Big0;
457 word1(rv) = Big1;
458 }
459 else
460 word0(rv) += P*Exp_msk1;
461 }
462 }
463 else if (e1 < 0) {
464 e1 = -e1;
465 if ( (i = e1 & 15) !=0)
466 dval(rv) /= tens[i];
467 if (e1 >>= 4) {
468 if (e1 >= 1 << n_bigtens)
469 goto undfl;
470 #ifdef Avoid_Underflow
471 if (e1 & Scale_Bit)
472 scale = 2*P;
473 for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
474 if (e1 & 1)
475 dval(rv) *= tinytens[j];
476 if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
477 >> Exp_shift)) > 0) {
478 /* scaled rv is denormal; zap j low bits */
479 if (j >= 32) {
480 word1(rv) = 0;
481 if (j >= 53)
482 word0(rv) = (P+2)*Exp_msk1;
483 else
484 word0(rv) &= 0xffffffff << (j-32);
485 }
486 else
487 word1(rv) &= 0xffffffff << j;
488 }
489 #else
490 for(j = 0; e1 > 1; j++, e1 >>= 1)
491 if (e1 & 1)
492 dval(rv) *= tinytens[j];
493 /* The last multiplication could underflow. */
494 dval(rv0) = dval(rv);
495 dval(rv) *= tinytens[j];
496 if (!dval(rv)) {
497 dval(rv) = 2.*dval(rv0);
498 dval(rv) *= tinytens[j];
499 #endif
500 if (!dval(rv)) {
501 undfl:
502 dval(rv) = 0.;
503 #ifndef NO_ERRNO
504 errno = ERANGE;
505 #endif
506 if (bd0)
507 goto retfree;
508 goto ret;
509 }
510 #ifndef Avoid_Underflow
511 word0(rv) = Tiny0;
512 word1(rv) = Tiny1;
513 /* The refinement below will clean
514 * this approximation up.
515 */
516 }
517 #endif
518 }
519 }
520
521 /* Now the hard part -- adjusting rv to the correct value.*/
522
523 /* Put digits into bd: true value = bd * 10^e */
524
525 bd0 = s2b(s0, nd0, nd, y);
526 if (bd0 == NULL)
527 goto ovfl;
528
529 for(;;) {
530 bd = Balloc(bd0->k);
531 if (bd == NULL)
532 goto ovfl;
533 Bcopy(bd, bd0);
534 bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
535 if (bb == NULL)
536 goto ovfl;
537 bs = i2b(1);
538 if (bs == NULL)
539 goto ovfl;
540
541 if (e >= 0) {
542 bb2 = bb5 = 0;
543 bd2 = bd5 = e;
544 }
545 else {
546 bb2 = bb5 = -e;
547 bd2 = bd5 = 0;
548 }
549 if (bbe >= 0)
550 bb2 += bbe;
551 else
552 bd2 -= bbe;
553 bs2 = bb2;
554 #ifdef Honor_FLT_ROUNDS
555 if (rounding != 1)
556 bs2++;
557 #endif
558 #ifdef Avoid_Underflow
559 j = bbe - scale;
560 i = j + bbbits - 1; /* logb(rv) */
561 if (i < Emin) /* denormal */
562 j += P - Emin;
563 else
564 j = P + 1 - bbbits;
565 #else /*Avoid_Underflow*/
566 #ifdef Sudden_Underflow
567 #ifdef IBM
568 j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
569 #else
570 j = P + 1 - bbbits;
571 #endif
572 #else /*Sudden_Underflow*/
573 j = bbe;
574 i = j + bbbits - 1; /* logb(rv) */
575 if (i < Emin) /* denormal */
576 j += P - Emin;
577 else
578 j = P + 1 - bbbits;
579 #endif /*Sudden_Underflow*/
580 #endif /*Avoid_Underflow*/
581 bb2 += j;
582 bd2 += j;
583 #ifdef Avoid_Underflow
584 bd2 += scale;
585 #endif
586 i = bb2 < bd2 ? bb2 : bd2;
587 if (i > bs2)
588 i = bs2;
589 if (i > 0) {
590 bb2 -= i;
591 bd2 -= i;
592 bs2 -= i;
593 }
594 if (bb5 > 0) {
595 bs = pow5mult(bs, bb5);
596 if (bs == NULL)
597 goto ovfl;
598 bb1 = mult(bs, bb);
599 if (bb1 == NULL)
600 goto ovfl;
601 Bfree(bb);
602 bb = bb1;
603 }
604 if (bb2 > 0) {
605 bb = lshift(bb, bb2);
606 if (bb == NULL)
607 goto ovfl;
608 }
609 if (bd5 > 0) {
610 bd = pow5mult(bd, bd5);
611 if (bd == NULL)
612 goto ovfl;
613 }
614 if (bd2 > 0) {
615 bd = lshift(bd, bd2);
616 if (bd == NULL)
617 goto ovfl;
618 }
619 if (bs2 > 0) {
620 bs = lshift(bs, bs2);
621 if (bs == NULL)
622 goto ovfl;
623 }
624 delta = diff(bb, bd);
625 if (delta == NULL)
626 goto ovfl;
627 dsign = delta->sign;
628 delta->sign = 0;
629 i = cmp(delta, bs);
630 #ifdef Honor_FLT_ROUNDS
631 if (rounding != 1) {
632 if (i < 0) {
633 /* Error is less than an ulp */
634 if (!delta->x[0] && delta->wds <= 1) {
635 /* exact */
636 #ifdef SET_INEXACT
637 inexact = 0;
638 #endif
639 break;
640 }
641 if (rounding) {
642 if (dsign) {
643 adj = 1.;
644 goto apply_adj;
645 }
646 }
647 else if (!dsign) {
648 adj = -1.;
649 if (!word1(rv)
650 && !(word0(rv) & Frac_mask)) {
651 y = word0(rv) & Exp_mask;
652 #ifdef Avoid_Underflow
653 if (!scale || y > 2*P*Exp_msk1)
654 #else
655 if (y)
656 #endif
657 {
658 delta = lshift(delta,Log2P);
659 if (cmp(delta, bs) <= 0)
660 adj = -0.5;
661 }
662 }
663 apply_adj:
664 #ifdef Avoid_Underflow
665 if (scale && (y = word0(rv) & Exp_mask)
666 <= 2*P*Exp_msk1)
667 word0(adj) += (2*P+1)*Exp_msk1 - y;
668 #else
669 #ifdef Sudden_Underflow
670 if ((word0(rv) & Exp_mask) <=
671 P*Exp_msk1) {
672 word0(rv) += P*Exp_msk1;
673 dval(rv) += adj*ulp(dval(rv));
674 word0(rv) -= P*Exp_msk1;
675 }
676 else
677 #endif /*Sudden_Underflow*/
678 #endif /*Avoid_Underflow*/
679 dval(rv) += adj*ulp(dval(rv));
680 }
681 break;
682 }
683 adj = ratio(delta, bs);
684 if (adj < 1.)
685 adj = 1.;
686 if (adj <= 0x7ffffffe) {
687 /* adj = rounding ? ceil(adj) : floor(adj); */
688 y = adj;
689 if (y != adj) {
690 if (!((rounding>>1) ^ dsign))
691 y++;
692 adj = y;
693 }
694 }
695 #ifdef Avoid_Underflow
696 if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
697 word0(adj) += (2*P+1)*Exp_msk1 - y;
698 #else
699 #ifdef Sudden_Underflow
700 if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
701 word0(rv) += P*Exp_msk1;
702 adj *= ulp(dval(rv));
703 if (dsign)
704 dval(rv) += adj;
705 else
706 dval(rv) -= adj;
707 word0(rv) -= P*Exp_msk1;
708 goto cont;
709 }
710 #endif /*Sudden_Underflow*/
711 #endif /*Avoid_Underflow*/
712 adj *= ulp(dval(rv));
713 if (dsign)
714 dval(rv) += adj;
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 if (word0(rv) == 0 && word1(rv) == 0)
996 errno = ERANGE;
997 #endif
998 }
999 #endif /* Avoid_Underflow */
1000 #ifdef SET_INEXACT
1001 if (inexact && !(word0(rv) & Exp_mask)) {
1002 /* set underflow bit */
1003 dval(rv0) = 1e-300;
1004 dval(rv0) *= dval(rv0);
1005 }
1006 #endif
1007 retfree:
1008 Bfree(bb);
1009 Bfree(bd);
1010 Bfree(bs);
1011 Bfree(bd0);
1012 Bfree(delta);
1013 ret:
1014 if (se)
1015 *se = __UNCONST(s);
1016 return sign ? -dval(rv) : dval(rv);
1017 }
1018
NetBSD on the FriendlyARM MINI2440