floatconv.c

早期freebsd实现

#ifndef RND_PRODQUOT
                && FLT_ROUNDS == 1
#endif
                        ) {
                if (!e)
                        goto ret;
                if (e > 0) {
                        if (e <= Ten_pmax) {
#ifdef VAX
                                goto vax_ovfl_check;
#else
                                /* rv = */ rounded_product(rv, tens[e]);
                                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:
                                word0(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;
                                word0(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;
#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. */
                                word0(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) */
                                        word0(rv) = Big0;
                                        word1(rv) = Big1;
                                        }
                                else
                                        word0(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;
                                        }
                                word0(rv) = Tiny0;
                                word1(rv) = 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(s0, nd0, nd, y);

        for(;;) {
                bd = Balloc(bd0->k);
                Bcopy(bd, bd0);
                bb = d2b(rv, &bbe, &bbbits);    /* rv = bb * 2^bbe */
                bs = i2b(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) {
                        bs = pow5mult(bs, bb5);
                        bb1 = mult(bs, bb);
                        Bfree(bb);
                        bb = bb1;
                        }
                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(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*/
                                        word0(rv) = (word0(rv) & Exp_mask)
                                                + Exp_msk1
#ifdef IBM
                                                | Exp_msk1 >> 4
#endif
                                                ;
                                        word1(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
                                word0(rv) = L | Bndry_mask1;
                                word1(rv) = 0xffffffff;
#ifdef IBM
                                goto cont;
#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;
                        word0(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;
                                word0(rv) = Big0;
                                word1(rv) = Big1;
                                goto cont;
                                }
                        else
                                word0(rv) += P*Exp_msk1;
                        }
                else {
#ifdef Sudden_Underflow
                        if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
                                rv0 = rv;
                                word0(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;
                                        word0(rv) = Tiny0;
                                        word1(rv) = Tiny1;
                                        goto cont;
                                        }
                                else
                                        word0(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 = (long)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;
                        }
 cont:
                Bfree(bb);
                Bfree(bd);
                Bfree(bs);
                Bfree(delta);
                }
        Bfree(bb);
        Bfree(bd);
        Bfree(bs);
        Bfree(bd0);
        Bfree(delta);
 ret:
        if (se)
                *se = (char *)s;
        return sign ? -rv : rv;
        }

 static int
quorem
#ifdef KR_headers
        (b, S) Bigint *b, *S;
#else
        (Bigint *b, Bigint *S)
#endif
{
        int n;
        long borrow, y;
        unsigned long carry, q, ys;
        unsigned long *bx, *bxe, *sx, *sxe;
#ifdef Pack_32
        long z;
        unsigned long si, zs;
#endif

        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 {
#ifdef Pack_32
                        si = *sx++;