![[BACK]](/icons/back.gif){kind=icon}
Return to strtod.c CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [cvs.NetBSD.org] / src / lib / libc / gdtoa|
Return to strtod.c CVS log | Up to [cvs.NetBSD.org] / src / lib / libc / gdtoa |
1.5 ! christos 1: /* $NetBSD: strtod.c,v 1.4 2006/06/02 19:46:56 mrg Exp $ */
1.1 kleink 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:
1.3 kleink 63: #ifndef __HAVE_LONG_DOUBLE
64: __strong_alias(_strtold, strtod)
65: __weak_alias(strtold, _strtold)
66: #endif
67:
1.1 kleink 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;
1.4 mrg 85: Bigint *bb = NULL, *bb1, *bd0;
1.2 kleink 86: Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
1.1 kleink 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;
1.2 kleink 99: /* FALLTHROUGH */
1.1 kleink 100: case '+':
101: if (*++s)
102: goto break2;
1.2 kleink 103: /* FALLTHROUGH */
1.1 kleink 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 };
1.2 kleink 121: Long expt;
1.1 kleink 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
1.2 kleink 137: switch((i = gethex(&s, &fpi1, &expt, &bb, sign)) & STRTOG_Retmask) {
1.1 kleink 138: case STRTOG_NoNumber:
139: s = s00;
140: sign = 0;
1.2 kleink 141: /* FALLTHROUGH */
1.1 kleink 142: case STRTOG_Zero:
143: break;
144: default:
145: if (bb) {
146: copybits(bits, fpi.nbits, bb);
147: Bfree(bb);
148: }
1.2 kleink 149: ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
1.1 kleink 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;
1.2 kleink 216: /* FALLTHROUGH */
1.1 kleink 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: }
1.2 kleink 443: e1 = (unsigned int)e1 >> 4;
444: for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
1.1 kleink 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;
1.2 kleink 473: for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
1.1 kleink 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
1.2 kleink 484: word0(rv) &= 0xffffffff << (j-32);
1.1 kleink 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);
1.5 ! christos 526: if (bd0 == NULL)
! 527: goto ovfl;
1.1 kleink 528:
529: for(;;) {
530: bd = Balloc(bd0->k);
1.5 ! christos 531: if (bd == NULL)
! 532: goto ovfl;
1.1 kleink 533: Bcopy(bd, bd0);
534: bb = d2b(dval(rv), &bbe, &bbbits); /* rv = bb * 2^bbe */
1.5 ! christos 535: if (bb == NULL)
! 536: goto ovfl;
1.1 kleink 537: bs = i2b(1);
1.5 ! christos 538: if (bs == NULL)
! 539: goto ovfl;
1.1 kleink 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);
1.5 ! christos 596: if (bs == NULL)
! 597: goto ovfl;
1.1 kleink 598: bb1 = mult(bs, bb);
1.5 ! christos 599: if (bb1 == NULL)
! 600: goto ovfl;
1.1 kleink 601: Bfree(bb);
602: bb = bb1;
603: }
1.5 ! christos 604: if (bb2 > 0) {
1.1 kleink 605: bb = lshift(bb, bb2);
1.5 ! christos 606: if (bb == NULL)
! 607: goto ovfl;
! 608: }
! 609: if (bd5 > 0) {
1.1 kleink 610: bd = pow5mult(bd, bd5);
1.5 ! christos 611: if (bd == NULL)
! 612: goto ovfl;
! 613: }
! 614: if (bd2 > 0) {
1.1 kleink 615: bd = lshift(bd, bd2);
1.5 ! christos 616: if (bd == NULL)
! 617: goto ovfl;
! 618: }
! 619: if (bs2 > 0) {
1.1 kleink 620: bs = lshift(bs, bs2);
1.5 ! christos 621: if (bs == NULL)
! 622: goto ovfl;
! 623: }
1.1 kleink 624: delta = diff(bb, bd);
1.5 ! christos 625: if (delta == NULL)
! 626: goto ovfl;
1.1 kleink 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) {
1.2 kleink 899: if ((z = aadj) == 0)
1.1 kleink 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)
1.2 kleink 1015: *se = __UNCONST(s);
1.1 kleink 1016: return sign ? -dval(rv) : dval(rv);
1017: }
1018: