![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to _dtoa.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 _dtoa.c CVS log | Up to [Research Unix] / researchv10no / libc / stdio |
1.1 ! root 1: #include "fconv.h" ! 2: ! 3: static int quorem(Bigint *, Bigint *); ! 4: ! 5: /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. ! 6: * ! 7: * Inspired by "How to Print Floating-Point Numbers Accurately" by ! 8: * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101]. ! 9: * ! 10: * Modifications: ! 11: * 1. Rather than iterating, we use a simple numeric overestimate ! 12: * to determine k = floor(log10(d)). We scale relevant ! 13: * quantities using O(log2(k)) rather than O(k) multiplications. ! 14: * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't ! 15: * try to generate digits strictly left to right. Instead, we ! 16: * compute with fewer bits and propagate the carry if necessary ! 17: * when rounding the final digit up. This is often faster. ! 18: * 3. Under the assumption that input will be rounded nearest, ! 19: * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. ! 20: * That is, we allow equality in stopping tests when the ! 21: * round-nearest rule will give the same floating-point value ! 22: * as would satisfaction of the stopping test with strict ! 23: * inequality. ! 24: * 4. We remove common factors of powers of 2 from relevant ! 25: * quantities. ! 26: * 5. When converting floating-point integers less than 1e16, ! 27: * we use floating-point arithmetic rather than resorting ! 28: * to multiple-precision integers. ! 29: * 6. When asked to produce fewer than 15 digits, we first try ! 30: * to get by with floating-point arithmetic; we resort to ! 31: * multiple-precision integer arithmetic only if we cannot ! 32: * guarantee that the floating-point calculation has given ! 33: * the correctly rounded result. For k requested digits and ! 34: * "uniformly" distributed input, the probability is ! 35: * something like 10^(k-15) that we must resort to the long ! 36: * calculation. ! 37: */ ! 38: ! 39: char * ! 40: _dtoa(double darg, int mode, int ndigits, int *decpt, int *sign, char **rve) ! 41: { ! 42: /* Arguments ndigits, decpt, sign are similar to those ! 43: of ecvt and fcvt; trailing zeros are suppressed from ! 44: the returned string. If not null, *rve is set to point ! 45: to the end of the return value. If d is +-Infinity or NaN, ! 46: then *decpt is set to 9999. ! 47: ! 48: mode: ! 49: 0 ==> shortest string that yields d when read in ! 50: and rounded to nearest. ! 51: 1 ==> like 0, but with Steele & White stopping rule; ! 52: e.g. with IEEE P754 arithmetic , mode 0 gives ! 53: 1e23 whereas mode 1 gives 9.999999999999999e22. ! 54: 2 ==> max(1,ndigits) significant digits. This gives a ! 55: return value similar to that of ecvt, except ! 56: that trailing zeros are suppressed. ! 57: 3 ==> through ndigits past the decimal point. This ! 58: gives a return value similar to that from fcvt, ! 59: except that trailing zeros are suppressed, and ! 60: ndigits can be negative. ! 61: 4-9 should give the same return values as 2-3, i.e., ! 62: 4 <= mode <= 9 ==> same return as mode ! 63: 2 + (mode & 1). These modes are mainly for ! 64: debugging; often they run slower but sometimes ! 65: faster than modes 2-3. ! 66: 4,5,8,9 ==> left-to-right digit generation. ! 67: 6-9 ==> don't try fast floating-point estimate ! 68: (if applicable). ! 69: ! 70: Values of mode other than 0-9 are treated as mode 0. ! 71: ! 72: Sufficient space is allocated to the return value ! 73: to hold the suppressed trailing zeros. ! 74: */ ! 75: ! 76: int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, ! 77: j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, ! 78: spec_case, try_quick; ! 79: long L; ! 80: #ifndef Sudden_Underflow ! 81: int denorm; ! 82: unsigned long x; ! 83: #endif ! 84: Bigint *b, *b1, *delta, *mlo, *mhi, *S; ! 85: double ds; ! 86: Dul d2, eps; ! 87: char *s, *s0; ! 88: static Bigint *result; ! 89: static int result_k; ! 90: Dul d; ! 91: ! 92: d.d = darg; ! 93: if (result) { ! 94: result->k = result_k; ! 95: result->maxwds = 1 << result_k; ! 96: Bfree(result); ! 97: result = 0; ! 98: } ! 99: ! 100: if (word0(d) & Sign_bit) { ! 101: /* set sign for everything, including 0's and NaNs */ ! 102: *sign = 1; ! 103: word0(d) &= ~Sign_bit; /* clear sign bit */ ! 104: } ! 105: else ! 106: *sign = 0; ! 107: ! 108: #if defined(IEEE_Arith) + defined(VAX) ! 109: #ifdef IEEE_Arith ! 110: if ((word0(d) & Exp_mask) == Exp_mask) ! 111: #else ! 112: if (word0(d) == 0x8000) ! 113: #endif ! 114: { ! 115: /* Infinity or NaN */ ! 116: *decpt = 9999; ! 117: s = ! 118: #ifdef IEEE_Arith ! 119: !word1(d) && !(word0(d) & 0xfffff) ? "Infinity" : ! 120: #endif ! 121: "NaN"; ! 122: if (rve) ! 123: *rve = ! 124: #ifdef IEEE_Arith ! 125: s[3] ? s + 8 : ! 126: #endif ! 127: s + 3; ! 128: return s; ! 129: } ! 130: #endif ! 131: #ifdef IBM ! 132: d.d += 0; /* normalize */ ! 133: #endif ! 134: if (!d.d) { ! 135: *decpt = 1; ! 136: s = "0"; ! 137: if (rve) ! 138: *rve = s + 1; ! 139: return s; ! 140: } ! 141: ! 142: b = d2b(d.d, &be, &bbits); ! 143: #ifdef Sudden_Underflow ! 144: i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1)); ! 145: #else ! 146: if (i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1))) { ! 147: #endif ! 148: d2.d = d.d; ! 149: word0(d2) &= Frac_mask1; ! 150: word0(d2) |= Exp_11; ! 151: #ifdef IBM ! 152: if (j = 11 - hi0bits(word0(d2) & Frac_mask)) ! 153: d2.d /= 1 << j; ! 154: #endif ! 155: ! 156: /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 ! 157: * log10(x) = log(x) / log(10) ! 158: * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) ! 159: * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) ! 160: * ! 161: * This suggests computing an approximation k to log10(d) by ! 162: * ! 163: * k = (i - Bias)*0.301029995663981 ! 164: * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); ! 165: * ! 166: * We want k to be too large rather than too small. ! 167: * The error in the first-order Taylor series approximation ! 168: * is in our favor, so we just round up the constant enough ! 169: * to compensate for any error in the multiplication of ! 170: * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, ! 171: * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, ! 172: * adding 1e-13 to the constant term more than suffices. ! 173: * Hence we adjust the constant term to 0.1760912590558. ! 174: * (We could get a more accurate k by invoking log10, ! 175: * but this is probably not worthwhile.) ! 176: */ ! 177: ! 178: i -= Bias; ! 179: #ifdef IBM ! 180: i <<= 2; ! 181: i += j; ! 182: #endif ! 183: #ifndef Sudden_Underflow ! 184: denorm = 0; ! 185: } ! 186: else { ! 187: /* d is denormalized */ ! 188: ! 189: i = bbits + be + (Bias + (P-1) - 1); ! 190: x = i > 32 ? word0(d) << 64 - i | word1(d) >> i - 32 ! 191: : word1(d) << 32 - i; ! 192: d2.d = x; ! 193: word0(d2) -= 31*Exp_msk1; /* adjust exponent */ ! 194: i -= (Bias + (P-1) - 1) + 1; ! 195: denorm = 1; ! 196: } ! 197: #endif ! 198: ds = (d2.d-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981; ! 199: k = (int)ds; ! 200: if (ds < 0. && ds != k) ! 201: k--; /* want k = floor(ds) */ ! 202: k_check = 1; ! 203: if (k >= 0 && k <= Ten_pmax) { ! 204: if (d.d < tens[k]) ! 205: k--; ! 206: k_check = 0; ! 207: } ! 208: j = bbits - i - 1; ! 209: if (j >= 0) { ! 210: b2 = 0; ! 211: s2 = j; ! 212: } ! 213: else { ! 214: b2 = -j; ! 215: s2 = 0; ! 216: } ! 217: if (k >= 0) { ! 218: b5 = 0; ! 219: s5 = k; ! 220: s2 += k; ! 221: } ! 222: else { ! 223: b2 -= k; ! 224: b5 = -k; ! 225: s5 = 0; ! 226: } ! 227: if (mode < 0 || mode > 9) ! 228: mode = 0; ! 229: try_quick = 1; ! 230: if (mode > 5) { ! 231: mode -= 4; ! 232: try_quick = 0; ! 233: } ! 234: leftright = 1; ! 235: switch(mode) { ! 236: case 0: ! 237: case 1: ! 238: ilim = ilim1 = -1; ! 239: i = 18; ! 240: ndigits = 0; ! 241: break; ! 242: case 2: ! 243: leftright = 0; ! 244: /* no break */ ! 245: case 4: ! 246: if (ndigits <= 0) ! 247: ndigits = 1; ! 248: ilim = ilim1 = i = ndigits; ! 249: break; ! 250: case 3: ! 251: leftright = 0; ! 252: /* no break */ ! 253: case 5: ! 254: i = ndigits + k + 1; ! 255: ilim = i; ! 256: ilim1 = i - 1; ! 257: if (i <= 0) ! 258: i = 1; ! 259: } ! 260: j = sizeof(unsigned long); ! 261: for(result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i; ! 262: j <<= 1) result_k++; ! 263: result = Balloc(result_k); ! 264: s = s0 = (char *)result; ! 265: ! 266: if (ilim >= 0 && ilim <= Quick_max && try_quick) { ! 267: ! 268: /* Try to get by with floating-point arithmetic. */ ! 269: ! 270: i = 0; ! 271: d2.d = d.d; ! 272: k0 = k; ! 273: ilim0 = ilim; ! 274: ieps = 2; /* conservative */ ! 275: if (k > 0) { ! 276: ds = tens[k&0xf]; ! 277: j = k >> 4; ! 278: if (j & Bletch) { ! 279: /* prevent overflows */ ! 280: j &= Bletch - 1; ! 281: d.d /= bigtens[n_bigtens-1]; ! 282: ieps++; ! 283: } ! 284: for(; j; j >>= 1, i++) ! 285: if (j & 1) { ! 286: ieps++; ! 287: ds *= bigtens[i]; ! 288: } ! 289: d.d /= ds; ! 290: } ! 291: else if (j1 = -k) { ! 292: d.d *= tens[j1 & 0xf]; ! 293: for(j = j1 >> 4; j; j >>= 1, i++) ! 294: if (j & 1) { ! 295: ieps++; ! 296: d.d *= bigtens[i]; ! 297: } ! 298: } ! 299: if (k_check && d.d < 1. && ilim > 0) { ! 300: if (ilim1 <= 0) ! 301: goto fast_failed; ! 302: ilim = ilim1; ! 303: k--; ! 304: d.d *= 10.; ! 305: ieps++; ! 306: } ! 307: eps.d = ieps*d.d + 7.; ! 308: word0(eps) -= (P-1)*Exp_msk1; ! 309: if (ilim == 0) { ! 310: S = mhi = 0; ! 311: d.d -= 5.; ! 312: if (d.d > eps.d) ! 313: goto one_digit; ! 314: if (d.d < -eps.d) ! 315: goto no_digits; ! 316: goto fast_failed; ! 317: } ! 318: #ifndef No_leftright ! 319: if (leftright) { ! 320: /* Use Steele & White method of only ! 321: * generating digits needed. ! 322: */ ! 323: eps.d = 0.5/tens[ilim-1] - eps.d; ! 324: for(i = 0;;) { ! 325: L = d.d; ! 326: d.d -= L; ! 327: *s++ = '0' + (int)L; ! 328: if (d.d < eps.d) ! 329: goto ret1; ! 330: if (1. - d.d < eps.d) ! 331: goto bump_up; ! 332: if (++i >= ilim) ! 333: break; ! 334: eps.d *= 10.; ! 335: d.d *= 10.; ! 336: } ! 337: } ! 338: else { ! 339: #endif ! 340: /* Generate ilim digits, then fix them up. */ ! 341: eps.d *= tens[ilim-1]; ! 342: for(i = 1;; i++, d.d *= 10.) { ! 343: L = d.d; ! 344: d.d -= L; ! 345: *s++ = '0' + (int)L; ! 346: if (i == ilim) { ! 347: if (d.d > 0.5 + eps.d) ! 348: goto bump_up; ! 349: else if (d.d < 0.5 - eps.d) { ! 350: while(*--s == '0'); ! 351: s++; ! 352: goto ret1; ! 353: } ! 354: break; ! 355: } ! 356: } ! 357: #ifndef No_leftright ! 358: } ! 359: #endif ! 360: fast_failed: ! 361: s = s0; ! 362: d.d = d2.d; ! 363: k = k0; ! 364: ilim = ilim0; ! 365: } ! 366: ! 367: /* Do we have a "small" integer? */ ! 368: ! 369: if (be >= 0 && k <= Int_max) { ! 370: /* Yes. */ ! 371: ds = tens[k]; ! 372: if (ndigits < 0 && ilim <= 0) { ! 373: S = mhi = 0; ! 374: if (ilim < 0 || d.d <= 5*ds) ! 375: goto no_digits; ! 376: goto one_digit; ! 377: } ! 378: for(i = 1;; i++) { ! 379: L = d.d / ds; ! 380: d.d -= L*ds; ! 381: #ifdef Check_FLT_ROUNDS ! 382: /* If FLT_ROUNDS == 2, L will usually be high by 1 */ ! 383: if (d.d < 0) { ! 384: L--; ! 385: d.d += ds; ! 386: } ! 387: #endif ! 388: *s++ = '0' + (int)L; ! 389: if (i == ilim) { ! 390: d.d += d.d; ! 391: if (d.d > ds || d.d == ds && L & 1) { ! 392: bump_up: ! 393: while(*--s == '9') ! 394: if (s == s0) { ! 395: k++; ! 396: *s = '0'; ! 397: break; ! 398: } ! 399: ++*s++; ! 400: } ! 401: break; ! 402: } ! 403: if (!(d.d *= 10.)) ! 404: break; ! 405: } ! 406: goto ret1; ! 407: } ! 408: ! 409: m2 = b2; ! 410: m5 = b5; ! 411: mhi = mlo = 0; ! 412: if (leftright) { ! 413: if (mode < 2) { ! 414: i = ! 415: #ifndef Sudden_Underflow ! 416: denorm ? be + (Bias + (P-1) - 1 + 1) : ! 417: #endif ! 418: #ifdef IBM ! 419: 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3); ! 420: #else ! 421: 1 + P - bbits; ! 422: #endif ! 423: } ! 424: else { ! 425: j = ilim - 1; ! 426: if (m5 >= j) ! 427: m5 -= j; ! 428: else { ! 429: s5 += j -= m5; ! 430: b5 += j; ! 431: m5 = 0; ! 432: } ! 433: if ((i = ilim) < 0) { ! 434: m2 -= i; ! 435: i = 0; ! 436: } ! 437: } ! 438: b2 += i; ! 439: s2 += i; ! 440: mhi = i2b(1); ! 441: } ! 442: if (m2 > 0 && s2 > 0) { ! 443: i = m2 < s2 ? m2 : s2; ! 444: b2 -= i; ! 445: m2 -= i; ! 446: s2 -= i; ! 447: } ! 448: if (b5 > 0) { ! 449: if (leftright) { ! 450: if (m5 > 0) { ! 451: mhi = pow5mult(mhi, m5); ! 452: b1 = mult(mhi, b); ! 453: Bfree(b); ! 454: b = b1; ! 455: } ! 456: if (j = b5 - m5) ! 457: b = pow5mult(b, j); ! 458: } ! 459: else ! 460: b = pow5mult(b, b5); ! 461: } ! 462: S = i2b(1); ! 463: if (s5 > 0) ! 464: S = pow5mult(S, s5); ! 465: ! 466: /* Check for special case that d is a normalized power of 2. */ ! 467: ! 468: if (mode < 2) { ! 469: if (!word1(d) && !(word0(d) & Bndry_mask) ! 470: #ifndef Sudden_Underflow ! 471: && word0(d) & Exp_mask ! 472: #endif ! 473: ) { ! 474: /* The special case */ ! 475: b2 += Log2P; ! 476: s2 += Log2P; ! 477: spec_case = 1; ! 478: } ! 479: else ! 480: spec_case = 0; ! 481: } ! 482: ! 483: /* Arrange for convenient computation of quotients: ! 484: * shift left if necessary so divisor has 4 leading 0 bits. ! 485: * ! 486: * Perhaps we should just compute leading 28 bits of S once ! 487: * and for all and pass them and a shift to quorem, so it ! 488: * can do shifts and ors to compute the numerator for q. ! 489: */ ! 490: #ifdef Pack_32 ! 491: if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) ! 492: i = 32 - i; ! 493: #else ! 494: if (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) ! 495: i = 16 - i; ! 496: #endif ! 497: if (i > 4) { ! 498: i -= 4; ! 499: b2 += i; ! 500: m2 += i; ! 501: s2 += i; ! 502: } ! 503: else if (i < 4) { ! 504: i += 28; ! 505: b2 += i; ! 506: m2 += i; ! 507: s2 += i; ! 508: } ! 509: if (b2 > 0) ! 510: b = lshift(b, b2); ! 511: if (s2 > 0) ! 512: S = lshift(S, s2); ! 513: if (k_check) { ! 514: if (cmp(b,S) < 0) { ! 515: k--; ! 516: b = multadd(b, 10, 0); /* we botched the k estimate */ ! 517: if (leftright) ! 518: mhi = multadd(mhi, 10, 0); ! 519: ilim = ilim1; ! 520: } ! 521: } ! 522: if (ilim <= 0 && mode > 2) { ! 523: if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) { ! 524: /* no digits, fcvt style */ ! 525: no_digits: ! 526: k = -1 - ndigits; ! 527: goto ret; ! 528: } ! 529: one_digit: ! 530: *s++ = '1'; ! 531: k++; ! 532: goto ret; ! 533: } ! 534: if (leftright) { ! 535: if (m2 > 0) ! 536: mhi = lshift(mhi, m2); ! 537: ! 538: /* Compute mlo -- check for special case ! 539: * that d is a normalized power of 2. ! 540: */ ! 541: ! 542: mlo = mhi; ! 543: if (spec_case) { ! 544: mhi = Balloc(mhi->k); ! 545: Bcopy(mhi, mlo); ! 546: mhi = lshift(mhi, Log2P); ! 547: } ! 548: ! 549: for(i = 1;;i++) { ! 550: dig = quorem(b,S) + '0'; ! 551: /* Do we yet have the shortest decimal string ! 552: * that will round to d? ! 553: */ ! 554: j = cmp(b, mlo); ! 555: delta = diff(S, mhi); ! 556: j1 = delta->sign ? 1 : cmp(b, delta); ! 557: Bfree(delta); ! 558: #ifndef ROUND_BIASED ! 559: if (j1 == 0 && !mode && !(word1(d) & 1)) { ! 560: if (dig == '9') ! 561: goto round_9_up; ! 562: if (j > 0) ! 563: dig++; ! 564: *s++ = dig; ! 565: goto ret; ! 566: } ! 567: #endif ! 568: if (j < 0 || j == 0 && !mode ! 569: #ifndef ROUND_BIASED ! 570: && !(word1(d) & 1) ! 571: #endif ! 572: ) { ! 573: if (j1 > 0) { ! 574: b = lshift(b, 1); ! 575: j1 = cmp(b, S); ! 576: if ((j1 > 0 || j1 == 0 && dig & 1) ! 577: && dig++ == '9') ! 578: goto round_9_up; ! 579: } ! 580: *s++ = dig; ! 581: goto ret; ! 582: } ! 583: if (j1 > 0) { ! 584: if (dig == '9') { /* possible if i == 1 */ ! 585: round_9_up: ! 586: *s++ = '9'; ! 587: goto roundoff; ! 588: } ! 589: *s++ = dig + 1; ! 590: goto ret; ! 591: } ! 592: *s++ = dig; ! 593: if (i == ilim) ! 594: break; ! 595: b = multadd(b, 10, 0); ! 596: if (mlo == mhi) ! 597: mlo = mhi = multadd(mhi, 10, 0); ! 598: else { ! 599: mlo = multadd(mlo, 10, 0); ! 600: mhi = multadd(mhi, 10, 0); ! 601: } ! 602: } ! 603: } ! 604: else ! 605: for(i = 1;; i++) { ! 606: *s++ = dig = quorem(b,S) + '0'; ! 607: if (i >= ilim) ! 608: break; ! 609: b = multadd(b, 10, 0); ! 610: } ! 611: ! 612: /* Round off last digit */ ! 613: ! 614: b = lshift(b, 1); ! 615: j = cmp(b, S); ! 616: if (j > 0 || j == 0 && dig & 1) { ! 617: roundoff: ! 618: while(*--s == '9') ! 619: if (s == s0) { ! 620: k++; ! 621: *s++ = '1'; ! 622: goto ret; ! 623: } ! 624: ++*s++; ! 625: } ! 626: else { ! 627: while(*--s == '0'); ! 628: s++; ! 629: } ! 630: ret: ! 631: Bfree(S); ! 632: if (mhi) { ! 633: if (mlo && mlo != mhi) ! 634: Bfree(mlo); ! 635: Bfree(mhi); ! 636: } ! 637: ret1: ! 638: Bfree(b); ! 639: *s = 0; ! 640: *decpt = k + 1; ! 641: if (rve) ! 642: *rve = s; ! 643: return s0; ! 644: } ! 645: ! 646: static int ! 647: quorem(Bigint *b, Bigint *S) ! 648: { ! 649: int n; ! 650: long borrow, y; ! 651: unsigned long carry, q, ys; ! 652: unsigned long *bx, *bxe, *sx, *sxe; ! 653: #ifdef Pack_32 ! 654: long z; ! 655: unsigned long si, zs; ! 656: #endif ! 657: ! 658: n = S->wds; ! 659: #ifdef DEBUG ! 660: /*debug*/ if (b->wds > n) ! 661: /*debug*/ Bug("oversize b in quorem"); ! 662: #endif ! 663: if (b->wds < n) ! 664: return 0; ! 665: sx = S->x; ! 666: sxe = sx + --n; ! 667: bx = b->x; ! 668: bxe = bx + n; ! 669: q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ ! 670: #ifdef DEBUG ! 671: /*debug*/ if (q > 9) ! 672: /*debug*/ Bug("oversized quotient in quorem"); ! 673: #endif ! 674: if (q) { ! 675: borrow = 0; ! 676: carry = 0; ! 677: do { ! 678: #ifdef Pack_32 ! 679: si = *sx++; ! 680: ys = (si & 0xffff) * q + carry; ! 681: zs = (si >> 16) * q + (ys >> 16); ! 682: carry = zs >> 16; ! 683: y = (*bx & 0xffff) - (ys & 0xffff) + borrow; ! 684: borrow = y >> 16; ! 685: Sign_Extend(borrow, y); ! 686: z = (*bx >> 16) - (zs & 0xffff) + borrow; ! 687: borrow = z >> 16; ! 688: Sign_Extend(borrow, z); ! 689: Storeinc(bx, z, y); ! 690: #else ! 691: ys = *sx++ * q + carry; ! 692: carry = ys >> 16; ! 693: y = *bx - (ys & 0xffff) + borrow; ! 694: borrow = y >> 16; ! 695: Sign_Extend(borrow, y); ! 696: *bx++ = y & 0xffff; ! 697: #endif ! 698: } ! 699: while(sx <= sxe); ! 700: if (!*bxe) { ! 701: bx = b->x; ! 702: while(--bxe > bx && !*bxe) ! 703: --n; ! 704: b->wds = n; ! 705: } ! 706: } ! 707: if (cmp(b, S) >= 0) { ! 708: q++; ! 709: borrow = 0; ! 710: carry = 0; ! 711: bx = b->x; ! 712: sx = S->x; ! 713: do { ! 714: #ifdef Pack_32 ! 715: si = *sx++; ! 716: ys = (si & 0xffff) + carry; ! 717: zs = (si >> 16) + (ys >> 16); ! 718: carry = zs >> 16; ! 719: y = (*bx & 0xffff) - (ys & 0xffff) + borrow; ! 720: borrow = y >> 16; ! 721: Sign_Extend(borrow, y); ! 722: z = (*bx >> 16) - (zs & 0xffff) + borrow; ! 723: borrow = z >> 16; ! 724: Sign_Extend(borrow, z); ! 725: Storeinc(bx, z, y); ! 726: #else ! 727: ys = *sx++ + carry; ! 728: carry = ys >> 16; ! 729: y = *bx - (ys & 0xffff) + borrow; ! 730: borrow = y >> 16; ! 731: Sign_Extend(borrow, y); ! 732: *bx++ = y & 0xffff; ! 733: #endif ! 734: } ! 735: while(sx <= sxe); ! 736: bx = b->x; ! 737: bxe = bx + n; ! 738: if (!*bxe) { ! 739: while(--bxe > bx && !*bxe) ! 740: --n; ! 741: b->wds = n; ! 742: } ! 743: } ! 744: return q; ! 745: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.