📄 floatconv.c
字号:
goto ret;#endif } i = DBL_DIG - nd; if (e <= Ten_pmax + i) { /* A fancier test would sometimes let us do * this for larger i values. */ e -= i; rv *= tens[i];#ifdef VAX /* VAX exponent range is so narrow we must * worry about overflow here... */ vax_ovfl_check: addword0(rv, - P*Exp_msk1); /* rv = */ rounded_product(rv, tens[e]); if ((word0(rv) & Exp_mask) > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) goto ovfl; addword0(rv, P*Exp_msk1);#else /* rv = */ rounded_product(rv, tens[e]);#endif goto ret; } }#ifndef Inaccurate_Divide else if (e >= -Ten_pmax) { /* rv = */ rounded_quotient(rv, tens[-e]); goto ret; }#endif } e1 += nd - k; /* Get starting approximation = rv * 10**e1 */ if (e1 > 0) { if ((i = e1 & 15)) rv *= tens[i]; if (e1 &= ~15) { if (e1 > DBL_MAX_10_EXP) { ovfl: errno = ERANGE;#if defined(sun) && !defined(__svr4__)/* SunOS defines HUGE_VAL as __infinity(), which is in libm. */#undef HUGE_VAL#endif#ifndef HUGE_VAL#define HUGE_VAL 1.7976931348623157E+308#endif rv = HUGE_VAL; goto ret; } if (e1 >>= 4) { for(j = 0; e1 > 1; j++, e1 >>= 1) if (e1 & 1) rv *= bigtens[j]; /* The last multiplication could overflow. */ addword0(rv, -P*Exp_msk1); rv *= bigtens[j]; if ((z = word0(rv) & Exp_mask) > Exp_msk1*(DBL_MAX_EXP+Bias-P)) goto ovfl; if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { /* set to largest number */ /* (Can't trust DBL_MAX) */ setwords(rv, Big0, Big1); } else addword0(rv, P*Exp_msk1); } } } else if (e1 < 0) { e1 = -e1; if ((i = e1 & 15)) rv /= tens[i]; if (e1 &= ~15) { e1 >>= 4; for(j = 0; e1 > 1; j++, e1 >>= 1) if (e1 & 1) rv *= tinytens[j]; /* The last multiplication could underflow. */ rv0 = rv; rv *= tinytens[j]; if (!rv) { rv = 2.*rv0; rv *= tinytens[j]; if (!rv) { undfl: rv = 0.; errno = ERANGE; goto ret; } setwords(rv, Tiny0, Tiny1); /* The refinement below will clean * this approximation up. */ } } } /* Now the hard part -- adjusting rv to the correct value.*/ /* Put digits into bd: true value = bd * 10^e */ bd0 = s2b(bd0, s0, nd0, nd, y); bd = Brealloc(bd, bd0->k); for(;;) { Bcopy(bd, bd0); bb = d2b(bb, rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ bs = i2b(bs, 1); if (e >= 0) { bb2 = bb5 = 0; bd2 = bd5 = e; } else { bb2 = bb5 = -e; bd2 = bd5 = 0; } if (bbe >= 0) bb2 += bbe; else bd2 -= bbe; bs2 = bb2;#ifdef Sudden_Underflow#ifdef IBM j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);#else j = P + 1 - bbbits;#endif#else i = bbe + bbbits - 1; /* logb(rv) */ if (i < Emin) /* denormal */ j = bbe + (P-Emin); else j = P + 1 - bbbits;#endif bb2 += j; bd2 += j; i = bb2 < bd2 ? bb2 : bd2; if (i > bs2) i = bs2; if (i > 0) { bb2 -= i; bd2 -= i; bs2 -= i; } if (bb5 > 0) { Bigint *b_tmp; bs = pow5mult(bs, bb5); b_tmp = mult(b_avail, bs, bb); b_avail = bb; bb = b_tmp; } if (bb2 > 0) bb = lshift(bb, bb2); if (bd5 > 0) bd = pow5mult(bd, bd5); if (bd2 > 0) bd = lshift(bd, bd2); if (bs2 > 0) bs = lshift(bs, bs2); delta = diff(delta, bb, bd); dsign = delta->sign; delta->sign = 0; i = cmp(delta, bs); if (i < 0) { /* Error is less than half an ulp -- check for * special case of mantissa a power of two. */ if (dsign || word1(rv) || word0(rv) & Bndry_mask) break; delta = lshift(delta,Log2P); if (cmp(delta, bs) > 0) goto drop_down; break; } if (i == 0) { /* exactly half-way between */ if (dsign) { if ((word0(rv) & Bndry_mask1) == Bndry_mask1 && word1(rv) == 0xffffffff) { /*boundary case -- increment exponent*/ setword0(rv, (word0(rv) & Exp_mask) + Exp_msk1);#ifdef IBM setword0 (rv, word0(rv) | (Exp_msk1 >> 4));#endif setword1(rv, 0); break; } } else if (!(word0(rv) & Bndry_mask) && !word1(rv)) { drop_down: /* boundary case -- decrement exponent */#ifdef Sudden_Underflow L = word0(rv) & Exp_mask;#ifdef IBM if (L < Exp_msk1)#else if (L <= Exp_msk1)#endif goto undfl; L -= Exp_msk1;#else L = (word0(rv) & Exp_mask) - Exp_msk1;#endif setwords(rv, L | Bndry_mask1, 0xffffffff);#ifdef IBM continue;#else break;#endif }#ifndef ROUND_BIASED if (!(word1(rv) & LSB)) break;#endif if (dsign) rv += ulp(rv);#ifndef ROUND_BIASED else { rv -= ulp(rv);#ifndef Sudden_Underflow if (!rv) goto undfl;#endif }#endif break; } if ((aadj = ratio(delta, bs)) <= 2.) { if (dsign) aadj = aadj1 = 1.; else if (word1(rv) || word0(rv) & Bndry_mask) {#ifndef Sudden_Underflow if (word1(rv) == Tiny1 && !word0(rv)) goto undfl;#endif aadj = 1.; aadj1 = -1.; } else { /* special case -- power of FLT_RADIX to be */ /* rounded down... */ if (aadj < 2./FLT_RADIX) aadj = 1./FLT_RADIX; else aadj *= 0.5; aadj1 = -aadj; } } else { aadj *= 0.5; aadj1 = dsign ? aadj : -aadj;#ifdef Check_FLT_ROUNDS switch(FLT_ROUNDS) { case 2: /* towards +infinity */ aadj1 -= 0.5; break; case 0: /* towards 0 */ case 3: /* towards -infinity */ aadj1 += 0.5; }#else if (FLT_ROUNDS == 0) aadj1 += 0.5;#endif } y = word0(rv) & Exp_mask; /* Check for overflow */ if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { rv0 = rv; addword0(rv, - P*Exp_msk1); adj = aadj1 * ulp(rv); rv += adj; if ((word0(rv) & Exp_mask) >= Exp_msk1*(DBL_MAX_EXP+Bias-P)) { if (word0(rv0) == Big0 && word1(rv0) == Big1) goto ovfl; setwords(rv, Big0, Big1); continue; } else addword0(rv, P*Exp_msk1); } else {#ifdef Sudden_Underflow if ((word0(rv) & Exp_mask) <= P*Exp_msk1) { rv0 = rv; addword0(rv, P*Exp_msk1); adj = aadj1 * ulp(rv); rv += adj;#ifdef IBM if ((word0(rv) & Exp_mask) < P*Exp_msk1)#else if ((word0(rv) & Exp_mask) <= P*Exp_msk1)#endif { if (word0(rv0) == Tiny0 && word1(rv0) == Tiny1) goto undfl; setwords(rv, Tiny0, Tiny1); continue; } else addword0(rv, -P*Exp_msk1); } else { adj = aadj1 * ulp(rv); rv += adj; }#else /* Compute adj so that the IEEE rounding rules will * correctly round rv + adj in some half-way cases. * If rv * ulp(rv) is denormalized (i.e., * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid * trouble from bits lost to denormalization; * example: 1.2e-307 . */ if (y <= (P-1)*Exp_msk1 && aadj >= 1.) { aadj1 = (double)(int)(aadj + 0.5); if (!dsign) aadj1 = -aadj1; } adj = aadj1 * ulp(rv); rv += adj;#endif } z = word0(rv) & Exp_mask; if (y == z) { /* Can we stop now? */ L = (_G_int32_t)aadj; aadj -= L; /* The tolerances below are conservative. */ if (dsign || word1(rv) || word0(rv) & Bndry_mask) { if (aadj < .4999999 || aadj > .5000001) break; } else if (aadj < .4999999/FLT_RADIX) break; } } Bfree(bb); Bfree(bd); Bfree(bs); Bfree(bd0); Bfree(delta); Bfree(b_avail); ret: if (se) *se = (char *)s; return sign ? -rv : rv; }static intquorem#ifdef KR_headers (b, S) Bigint *b, *S;#else (Bigint *b, Bigint *S)#endif{ int n; _G_int32_t borrow, y; unsigned32 carry, q, ys; unsigned32 *bx, *bxe, *sx, *sxe; _G_int32_t z; unsigned32 si, zs; n = S->wds;#ifdef DEBUG /*debug*/ if (b->wds > n) /*debug*/ Bug("oversize b in quorem");#endif if (b->wds < n) return 0; sx = S->x; sxe = sx + --n; bx = b->x; bxe = bx + n; q = *bxe / (*sxe + 1); /* ensure q <= true quotient */#ifdef DEBUG /*debug*/ if (q > 9) /*debug*/ Bug("oversized quotient in quorem");#endif if (q) { borrow = 0; carry = 0; do { si = *sx++; ys = (si & 0xffff) * q + carry; zs = (si >> 16) * q + (ys >> 16); carry = zs >> 16; y = (*bx & 0xffff) - (ys & 0xffff) + borrow; borrow = y >> 16; Sign_Extend(borrow, y); z = (*bx >> 16) - (zs & 0xffff) + borrow; borrow = z >> 16; Sign_Extend(borrow, z); Storeinc(bx, z, y); } while(sx <= sxe); if (!*bxe) { bx = b->x; while(--bxe > bx && !*bxe) --n; b->wds = n; } } if (cmp(b, S) >= 0) { q++; borrow = 0; carry = 0; bx = b->x; sx = S->x; do {
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -