![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to strtod.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / libc / stdio|
Return to strtod.c CVS log | Up to [Research Unix] / researchv10no / libc / stdio |
1.1 ! root 1: #include "fconv.h" ! 2: ! 3: /* strtod for IEEE-, VAX-, and IBM-arithmetic machines (dmg). ! 4: * ! 5: * This strtod returns a nearest machine number to the input decimal ! 6: * string (or sets errno to ERANGE). With IEEE arithmetic, ties are ! 7: * broken by the IEEE round-even rule. Otherwise ties are broken by ! 8: * biased rounding (add half and chop). ! 9: * ! 10: * Inspired loosely by William D. Clinger's paper "How to Read Floating ! 11: * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. ! 12: * ! 13: * Modifications: ! 14: * ! 15: * 1. We only require IEEE, IBM, or VAX double-precision ! 16: * arithmetic (not IEEE double-extended). ! 17: * 2. We get by with floating-point arithmetic in a case that ! 18: * Clinger missed -- when we're computing d * 10^n ! 19: * for a small integer d and the integer n is not too ! 20: * much larger than 22 (the maximum integer k for which ! 21: * we can represent 10^k exactly), we may be able to ! 22: * compute (d*10^k) * 10^(e-k) with just one roundoff. ! 23: * 3. Rather than a bit-at-a-time adjustment of the binary ! 24: * result in the hard case, we use floating-point ! 25: * arithmetic to determine the adjustment to within ! 26: * one bit; only in really hard cases do we need to ! 27: * compute a second residual. ! 28: * 4. Because of 3., we don't need a large table of powers of 10 ! 29: * for ten-to-e (just some small tables, e.g. of 10^k ! 30: * for 0 <= k <= 22). ! 31: */ ! 32: ! 33: #ifdef RND_PRODQUOT ! 34: #define rounded_product(a,b) a = rnd_prod(a, b) ! 35: #define rounded_quotient(a,b) a = rnd_quot(a, b) ! 36: extern double rnd_prod(double, double), rnd_quot(double, double); ! 37: #else ! 38: #define rounded_product(a,b) a *= b ! 39: #define rounded_quotient(a,b) a /= b ! 40: #endif ! 41: ! 42: static double ! 43: ulp(double xarg) ! 44: { ! 45: register long L; ! 46: Dul a; ! 47: Dul x; ! 48: ! 49: x.d = xarg; ! 50: L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; ! 51: #ifndef Sudden_Underflow ! 52: if (L > 0) { ! 53: #endif ! 54: #ifdef IBM ! 55: L |= Exp_msk1 >> 4; ! 56: #endif ! 57: word0(a) = L; ! 58: word1(a) = 0; ! 59: #ifndef Sudden_Underflow ! 60: } ! 61: else { ! 62: L = -L >> Exp_shift; ! 63: if (L < Exp_shift) { ! 64: word0(a) = 0x80000 >> L; ! 65: word1(a) = 0; ! 66: } ! 67: else { ! 68: word0(a) = 0; ! 69: L -= Exp_shift; ! 70: word1(a) = L >= 31 ? 1 : 1 << 31 - L; ! 71: } ! 72: } ! 73: #endif ! 74: return a.d; ! 75: } ! 76: ! 77: static Bigint * ! 78: s2b(CONST char *s, int nd0, int nd, unsigned long y9) ! 79: { ! 80: Bigint *b; ! 81: int i, k; ! 82: long x, y; ! 83: ! 84: x = (nd + 8) / 9; ! 85: for(k = 0, y = 1; x > y; y <<= 1, k++) ; ! 86: #ifdef Pack_32 ! 87: b = Balloc(k); ! 88: b->x[0] = y9; ! 89: b->wds = 1; ! 90: #else ! 91: b = Balloc(k+1); ! 92: b->x[0] = y9 & 0xffff; ! 93: b->wds = (b->x[1] = y9 >> 16) ? 2 : 1; ! 94: #endif ! 95: ! 96: i = 9; ! 97: if (9 < nd0) { ! 98: s += 9; ! 99: do b = multadd(b, 10, *s++ - '0'); ! 100: while(++i < nd0); ! 101: s++; ! 102: } ! 103: else ! 104: s += 10; ! 105: for(; i < nd; i++) ! 106: b = multadd(b, 10, *s++ - '0'); ! 107: return b; ! 108: } ! 109: ! 110: static double ! 111: b2d(Bigint *a, int *e) ! 112: { ! 113: unsigned long *xa, *xa0, w, y, z; ! 114: int k; ! 115: Dul d; ! 116: #ifdef VAX ! 117: unsigned long d0, d1; ! 118: #else ! 119: #define d0 word0(d) ! 120: #define d1 word1(d) ! 121: #endif ! 122: ! 123: xa0 = a->x; ! 124: xa = xa0 + a->wds; ! 125: y = *--xa; ! 126: #ifdef DEBUG ! 127: if (!y) Bug("zero y in b2d"); ! 128: #endif ! 129: k = hi0bits(y); ! 130: *e = 32 - k; ! 131: #ifdef Pack_32 ! 132: if (k < Ebits) { ! 133: d0 = Exp_1 | y >> Ebits - k; ! 134: w = xa > xa0 ? *--xa : 0; ! 135: d1 = y << (32-Ebits) + k | w >> Ebits - k; ! 136: goto ret_d; ! 137: } ! 138: z = xa > xa0 ? *--xa : 0; ! 139: if (k -= Ebits) { ! 140: d0 = Exp_1 | y << k | z >> 32 - k; ! 141: y = xa > xa0 ? *--xa : 0; ! 142: d1 = z << k | y >> 32 - k; ! 143: } ! 144: else { ! 145: d0 = Exp_1 | y; ! 146: d1 = z; ! 147: } ! 148: #else ! 149: if (k < Ebits + 16) { ! 150: z = xa > xa0 ? *--xa : 0; ! 151: d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; ! 152: w = xa > xa0 ? *--xa : 0; ! 153: y = xa > xa0 ? *--xa : 0; ! 154: d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; ! 155: goto ret_d; ! 156: } ! 157: z = xa > xa0 ? *--xa : 0; ! 158: w = xa > xa0 ? *--xa : 0; ! 159: k -= Ebits + 16; ! 160: d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; ! 161: y = xa > xa0 ? *--xa : 0; ! 162: d1 = w << k + 16 | y << k; ! 163: #endif ! 164: ret_d: ! 165: #ifdef VAX ! 166: word0(d) = d0 >> 16 | d0 << 16; ! 167: word1(d) = d1 >> 16 | d1 << 16; ! 168: #else ! 169: #undef d0 ! 170: #undef d1 ! 171: #endif ! 172: return d.d; ! 173: } ! 174: ! 175: static double ! 176: ratio(Bigint *a, Bigint *b) ! 177: { ! 178: Dul da, db; ! 179: int k, ka, kb; ! 180: ! 181: da.d = b2d(a, &ka); ! 182: db.d = b2d(b, &kb); ! 183: #ifdef Pack_32 ! 184: k = ka - kb + 32*(a->wds - b->wds); ! 185: #else ! 186: k = ka - kb + 16*(a->wds - b->wds); ! 187: #endif ! 188: #ifdef IBM ! 189: if (k > 0) { ! 190: word0(da) += (k >> 2)*Exp_msk1; ! 191: if (k &= 3) ! 192: da *= 1 << k; ! 193: } ! 194: else { ! 195: k = -k; ! 196: word0(db) += (k >> 2)*Exp_msk1; ! 197: if (k &= 3) ! 198: db *= 1 << k; ! 199: } ! 200: #else ! 201: if (k > 0) ! 202: word0(da) += k*Exp_msk1; ! 203: else { ! 204: k = -k; ! 205: word0(db) += k*Exp_msk1; ! 206: } ! 207: #endif ! 208: return da.d / db.d; ! 209: } ! 210: ! 211: double ! 212: strtod(CONST char *s00, char **se) ! 213: { ! 214: int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, ! 215: e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; ! 216: CONST char *s, *s0, *s1; ! 217: double aadj, aadj1, adj; ! 218: Dul rv, rv0; ! 219: long L; ! 220: unsigned long y, z; ! 221: Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; ! 222: sign = nz0 = nz = 0; ! 223: rv.d = 0.; ! 224: for(s = s00;;s++) switch(*s) { ! 225: case '-': ! 226: sign = 1; ! 227: /* no break */ ! 228: case '+': ! 229: if (*++s) ! 230: goto break2; ! 231: /* no break */ ! 232: case 0: ! 233: s = s00; ! 234: goto ret; ! 235: case '\t': ! 236: case '\n': ! 237: case '\v': ! 238: case '\f': ! 239: case '\r': ! 240: case ' ': ! 241: continue; ! 242: default: ! 243: goto break2; ! 244: } ! 245: break2: ! 246: if (*s == '0') { ! 247: nz0 = 1; ! 248: while(*++s == '0') ; ! 249: if (!*s) ! 250: goto ret; ! 251: } ! 252: s0 = s; ! 253: y = z = 0; ! 254: for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) ! 255: if (nd < 9) ! 256: y = 10*y + c - '0'; ! 257: else if (nd < 16) ! 258: z = 10*z + c - '0'; ! 259: nd0 = nd; ! 260: if (c == '.') { ! 261: c = *++s; ! 262: if (!nd) { ! 263: for(; c == '0'; c = *++s) ! 264: nz++; ! 265: if (c > '0' && c <= '9') { ! 266: s0 = s; ! 267: nf += nz; ! 268: nz = 0; ! 269: goto have_dig; ! 270: } ! 271: goto dig_done; ! 272: } ! 273: for(; c >= '0' && c <= '9'; c = *++s) { ! 274: have_dig: ! 275: nz++; ! 276: if (c -= '0') { ! 277: nf += nz; ! 278: for(i = 1; i < nz; i++) ! 279: if (nd++ < 9) ! 280: y *= 10; ! 281: else if (nd <= DBL_DIG + 1) ! 282: z *= 10; ! 283: if (nd++ < 9) ! 284: y = 10*y + c; ! 285: else if (nd <= DBL_DIG + 1) ! 286: z = 10*z + c; ! 287: nz = 0; ! 288: } ! 289: } ! 290: } ! 291: dig_done: ! 292: e = 0; ! 293: if (c == 'e' || c == 'E') { ! 294: if (!nd && !nz && !nz0) { ! 295: s = s00; ! 296: goto ret; ! 297: } ! 298: s00 = s; ! 299: esign = 0; ! 300: switch(c = *++s) { ! 301: case '-': ! 302: esign = 1; ! 303: case '+': ! 304: c = *++s; ! 305: } ! 306: if (c >= '0' && c <= '9') { ! 307: while(c == '0') ! 308: c = *++s; ! 309: if (c > '0' && c <= '9') { ! 310: e = c - '0'; ! 311: s1 = s; ! 312: while((c = *++s) >= '0' && c <= '9') ! 313: e = 10*e + c - '0'; ! 314: if (s - s1 > 8) ! 315: /* Avoid confusion from exponents ! 316: * so large that e might overflow. ! 317: */ ! 318: e = 9999999; ! 319: if (esign) ! 320: e = -e; ! 321: } ! 322: else ! 323: e = 0; ! 324: } ! 325: else ! 326: s = s00; ! 327: } ! 328: if (!nd) { ! 329: if (!nz && !nz0) ! 330: s = s00; ! 331: goto ret; ! 332: } ! 333: e1 = e -= nf; ! 334: ! 335: /* Now we have nd0 digits, starting at s0, followed by a ! 336: * decimal point, followed by nd-nd0 digits. The number we're ! 337: * after is the integer represented by those digits times ! 338: * 10**e */ ! 339: ! 340: if (!nd0) ! 341: nd0 = nd; ! 342: k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; ! 343: rv.d = y; ! 344: if (k > 9) ! 345: rv.d = tens[k - 9] * rv.d + z; ! 346: if (nd <= DBL_DIG ! 347: #ifndef RND_PRODQUOT ! 348: && FLT_ROUNDS == 1 ! 349: #endif ! 350: ) { ! 351: if (!e) ! 352: goto ret; ! 353: if (e > 0) { ! 354: if (e <= Ten_pmax) { ! 355: #ifdef VAX ! 356: goto vax_ovfl_check; ! 357: #else ! 358: /* rv = */ rounded_product(rv.d, tens[e]); ! 359: goto ret; ! 360: #endif ! 361: } ! 362: i = DBL_DIG - nd; ! 363: if (e <= Ten_pmax + i) { ! 364: /* A fancier test would sometimes let us do ! 365: * this for larger i values. ! 366: */ ! 367: e -= i; ! 368: rv.d *= tens[i]; ! 369: #ifdef VAX ! 370: /* VAX exponent range is so narrow we must ! 371: * worry about overflow here... ! 372: */ ! 373: vax_ovfl_check: ! 374: word0(rv) -= P*Exp_msk1; ! 375: /* rv = */ rounded_product(rv.d, tens[e]); ! 376: if ((word0(rv) & Exp_mask) ! 377: > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) ! 378: goto ovfl; ! 379: word0(rv) += P*Exp_msk1; ! 380: #else ! 381: /* rv = */ rounded_product(rv.d, tens[e]); ! 382: #endif ! 383: goto ret; ! 384: } ! 385: } ! 386: else if (e >= -Ten_pmax) { ! 387: /* rv = */ rounded_quotient(rv.d, tens[-e]); ! 388: goto ret; ! 389: } ! 390: } ! 391: e1 += nd - k; ! 392: ! 393: /* Get starting approximation = rv * 10**e1 */ ! 394: ! 395: if (e1 > 0) { ! 396: if (i = e1 & 15) ! 397: rv.d *= tens[i]; ! 398: if (e1 &= ~15) { ! 399: if (e1 > DBL_MAX_10_EXP) { ! 400: ovfl: ! 401: errno = ERANGE; ! 402: rv.d = HUGE_VAL; ! 403: goto ret; ! 404: } ! 405: if (e1 >>= 4) { ! 406: for(j = 0; e1 > 1; j++, e1 >>= 1) ! 407: if (e1 & 1) ! 408: rv.d *= bigtens[j]; ! 409: /* The last multiplication could overflow. */ ! 410: word0(rv) -= P*Exp_msk1; ! 411: rv.d *= bigtens[j]; ! 412: if ((z = word0(rv) & Exp_mask) ! 413: > Exp_msk1*(DBL_MAX_EXP+Bias-P)) ! 414: goto ovfl; ! 415: if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { ! 416: /* set to largest number */ ! 417: /* (Can't trust DBL_MAX) */ ! 418: word0(rv) = Big0; ! 419: word1(rv) = Big1; ! 420: } ! 421: else ! 422: word0(rv) += P*Exp_msk1; ! 423: } ! 424: ! 425: } ! 426: } ! 427: else if (e1 < 0) { ! 428: e1 = -e1; ! 429: if (i = e1 & 15) ! 430: rv.d /= tens[i]; ! 431: if (e1 &= ~15) { ! 432: e1 >>= 4; ! 433: for(j = 0; e1 > 1; j++, e1 >>= 1) ! 434: if (e1 & 1) ! 435: rv.d *= tinytens[j]; ! 436: /* The last multiplication could underflow. */ ! 437: rv0.d = rv.d; ! 438: rv.d *= tinytens[j]; ! 439: if (rv.d == 0) { ! 440: rv.d = 2.*rv0.d; ! 441: rv.d *= tinytens[j]; ! 442: if (rv.d == 0) { ! 443: undfl: ! 444: rv.d = 0.; ! 445: errno = ERANGE; ! 446: goto ret; ! 447: } ! 448: word0(rv) = Tiny0; ! 449: word1(rv) = Tiny1; ! 450: /* The refinement below will clean ! 451: * this approximation up. ! 452: */ ! 453: } ! 454: } ! 455: } ! 456: ! 457: /* Now the hard part -- adjusting rv to the correct value.*/ ! 458: ! 459: /* Put digits into bd: true value = bd * 10^e */ ! 460: ! 461: bd0 = s2b(s0, nd0, nd, y); ! 462: ! 463: for(;;) { ! 464: bd = Balloc(bd0->k); ! 465: Bcopy(bd, bd0); ! 466: bb = d2b(rv.d, &bbe, &bbbits); /* rv = bb * 2^bbe */ ! 467: bs = i2b(1); ! 468: ! 469: if (e >= 0) { ! 470: bb2 = bb5 = 0; ! 471: bd2 = bd5 = e; ! 472: } ! 473: else { ! 474: bb2 = bb5 = -e; ! 475: bd2 = bd5 = 0; ! 476: } ! 477: if (bbe >= 0) ! 478: bb2 += bbe; ! 479: else ! 480: bd2 -= bbe; ! 481: bs2 = bb2; ! 482: #ifdef Sudden_Underflow ! 483: #ifdef IBM ! 484: j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); ! 485: #else ! 486: j = P + 1 - bbbits; ! 487: #endif ! 488: #else ! 489: i = bbe + bbbits - 1; /* logb(rv) */ ! 490: if (i < Emin) /* denormal */ ! 491: j = bbe + (P-Emin); ! 492: else ! 493: j = P + 1 - bbbits; ! 494: #endif ! 495: bb2 += j; ! 496: bd2 += j; ! 497: i = bb2 < bd2 ? bb2 : bd2; ! 498: if (i > bs2) ! 499: i = bs2; ! 500: if (i > 0) { ! 501: bb2 -= i; ! 502: bd2 -= i; ! 503: bs2 -= i; ! 504: } ! 505: if (bb5 > 0) { ! 506: bs = pow5mult(bs, bb5); ! 507: bb1 = mult(bs, bb); ! 508: Bfree(bb); ! 509: bb = bb1; ! 510: } ! 511: if (bb2 > 0) ! 512: bb = lshift(bb, bb2); ! 513: if (bd5 > 0) ! 514: bd = pow5mult(bd, bd5); ! 515: if (bd2 > 0) ! 516: bd = lshift(bd, bd2); ! 517: if (bs2 > 0) ! 518: bs = lshift(bs, bs2); ! 519: delta = diff(bb, bd); ! 520: dsign = delta->sign; ! 521: delta->sign = 0; ! 522: i = cmp(delta, bs); ! 523: if (i < 0) { ! 524: /* Error is less than half an ulp -- check for ! 525: * special case of mantissa a power of two. ! 526: */ ! 527: if (dsign || word1(rv) || word0(rv) & Bndry_mask) ! 528: break; ! 529: delta = lshift(delta,Log2P); ! 530: if (cmp(delta, bs) > 0) ! 531: goto drop_down; ! 532: break; ! 533: } ! 534: if (i == 0) { ! 535: /* exactly half-way between */ ! 536: if (dsign) { ! 537: if ((word0(rv) & Bndry_mask1) == Bndry_mask1 ! 538: && word1(rv) == 0xffffffff) { ! 539: /*boundary case -- increment exponent*/ ! 540: word0(rv) = (word0(rv) & Exp_mask) ! 541: + Exp_msk1 ! 542: #ifdef IBM ! 543: | Exp_msk1 >> 4 ! 544: #endif ! 545: ; ! 546: word1(rv) = 0; ! 547: break; ! 548: } ! 549: } ! 550: else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { ! 551: drop_down: ! 552: /* boundary case -- decrement exponent */ ! 553: #ifdef Sudden_Underflow ! 554: L = word0(rv) & Exp_mask; ! 555: #ifdef IBM ! 556: if (L < Exp_msk1) ! 557: #else ! 558: if (L <= Exp_msk1) ! 559: #endif ! 560: goto undfl; ! 561: L -= Exp_msk1; ! 562: #else ! 563: L = (word0(rv) & Exp_mask) - Exp_msk1; ! 564: #endif ! 565: word0(rv) = L | Bndry_mask1; ! 566: word1(rv) = 0xffffffff; ! 567: #ifdef IBM ! 568: goto cont; ! 569: #else ! 570: break; ! 571: #endif ! 572: } ! 573: #ifndef ROUND_BIASED ! 574: if (!(word1(rv) & LSB)) ! 575: break; ! 576: #endif ! 577: if (dsign) ! 578: rv.d += ulp(rv.d); ! 579: #ifndef ROUND_BIASED ! 580: else { ! 581: rv.d -= ulp(rv.d); ! 582: #ifndef Sudden_Underflow ! 583: if (rv.d == 0) ! 584: goto undfl; ! 585: #endif ! 586: } ! 587: #endif ! 588: break; ! 589: } ! 590: if ((aadj = ratio(delta, bs)) <= 2.) { ! 591: if (dsign) ! 592: aadj = aadj1 = 1.; ! 593: else if (word1(rv) || word0(rv) & Bndry_mask) { ! 594: #ifndef Sudden_Underflow ! 595: if (word1(rv) == Tiny1 && !word0(rv)) ! 596: goto undfl; ! 597: #endif ! 598: aadj = 1.; ! 599: aadj1 = -1.; ! 600: } ! 601: else { ! 602: /* special case -- power of FLT_RADIX to be */ ! 603: /* rounded down... */ ! 604: ! 605: if (aadj < 2./FLT_RADIX) ! 606: aadj = 1./FLT_RADIX; ! 607: else ! 608: aadj *= 0.5; ! 609: aadj1 = -aadj; ! 610: } ! 611: } ! 612: else { ! 613: aadj *= 0.5; ! 614: aadj1 = dsign ? aadj : -aadj; ! 615: #ifdef Check_FLT_ROUNDS ! 616: switch(FLT_ROUNDS) { ! 617: case 2: /* towards +infinity */ ! 618: aadj1 -= 0.5; ! 619: break; ! 620: case 0: /* towards 0 */ ! 621: case 3: /* towards -infinity */ ! 622: aadj1 += 0.5; ! 623: } ! 624: #else ! 625: if (FLT_ROUNDS == 0) ! 626: aadj1 += 0.5; ! 627: #endif ! 628: } ! 629: y = word0(rv) & Exp_mask; ! 630: ! 631: /* Check for overflow */ ! 632: ! 633: if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { ! 634: rv0.d = rv.d; ! 635: word0(rv) -= P*Exp_msk1; ! 636: adj = aadj1 * ulp(rv.d); ! 637: rv.d += adj; ! 638: if ((word0(rv) & Exp_mask) >= ! 639: Exp_msk1*(DBL_MAX_EXP+Bias-P)) { ! 640: if (word0(rv0) == Big0 && word1(rv0) == Big1) ! 641: goto ovfl; ! 642: word0(rv) = Big0; ! 643: word1(rv) = Big1; ! 644: goto cont; ! 645: } ! 646: else ! 647: word0(rv) += P*Exp_msk1; ! 648: } ! 649: else { ! 650: #ifdef Sudden_Underflow ! 651: if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { ! 652: rv0.d = rv.d; ! 653: word0(rv) += P*Exp_msk1; ! 654: adj = aadj1 * ulp(rv.d); ! 655: rv.d += adj; ! 656: #ifdef IBM ! 657: if ((word0(rv) & Exp_mask) < P*Exp_msk1) ! 658: #else ! 659: if ((word0(rv) & Exp_mask) <= P*Exp_msk1) ! 660: #endif ! 661: { ! 662: if (word0(rv0) == Tiny0 ! 663: && word1(rv0) == Tiny1) ! 664: goto undfl; ! 665: word0(rv) = Tiny0; ! 666: word1(rv) = Tiny1; ! 667: goto cont; ! 668: } ! 669: else ! 670: word0(rv) -= P*Exp_msk1; ! 671: } ! 672: else { ! 673: adj = aadj1 * ulp(rv.d); ! 674: rv.d += adj; ! 675: } ! 676: #else ! 677: /* Compute adj so that the IEEE rounding rules will ! 678: * correctly round rv + adj in some half-way cases. ! 679: * If rv * ulp(rv) is denormalized (i.e., ! 680: * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid ! 681: * trouble from bits lost to denormalization; ! 682: * example: 1.2e-307 . ! 683: */ ! 684: if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { ! 685: aadj1 = (double)(int)(aadj + 0.5); ! 686: if (!dsign) ! 687: aadj1 = -aadj1; ! 688: } ! 689: adj = aadj1 * ulp(rv.d); ! 690: rv.d += adj; ! 691: #endif ! 692: } ! 693: z = word0(rv) & Exp_mask; ! 694: if (y == z) { ! 695: /* Can we stop now? */ ! 696: L = aadj; ! 697: aadj -= L; ! 698: /* The tolerances below are conservative. */ ! 699: if (dsign || word1(rv) || word0(rv) & Bndry_mask) { ! 700: if (aadj < .4999999 || aadj > .5000001) ! 701: break; ! 702: } ! 703: else if (aadj < .4999999/FLT_RADIX) ! 704: break; ! 705: } ! 706: cont: ! 707: Bfree(bb); ! 708: Bfree(bd); ! 709: Bfree(bs); ! 710: Bfree(delta); ! 711: } ! 712: Bfree(bb); ! 713: Bfree(bd); ! 714: Bfree(bs); ! 715: Bfree(bd0); ! 716: Bfree(delta); ! 717: ret: ! 718: if (se) ! 719: *se = (char *)s; ! 720: return sign ? -rv.d : rv.d; ! 721: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.