double
ruby_strtod(const char *s00, char **se)
{
#ifdef Avoid_Underflow
    int scale;
#endif
    int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
         e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
    const char *s, *s0, *s1;
    double aadj, adj;
    double_u aadj1, rv, rv0;
    Long L;
    ULong y, z;
    Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
#ifdef SET_INEXACT
    int inexact, oldinexact;
#endif
#ifdef Honor_FLT_ROUNDS
    int rounding;
#endif
#ifdef USE_LOCALE
    const char *s2;
#endif

    errno = 0;
    sign = nz0 = nz = 0;
    dval(rv) = 0.;
    for (s = s00;;s++)
        switch (*s) {
          case '-':
            sign = 1;
            /* no break */
          case '+':
            if (*++s)
                goto break2;
            /* no break */
          case 0:
            goto ret0;
          case '\t':
          case '\n':
          case '\v':
          case '\f':
          case '\r':
          case ' ':
            continue;
          default:
            goto break2;
        }
break2:
    if (*s == '0') {
	if (s[1] == 'x' || s[1] == 'X') {
	    static const char hexdigit[] = "0123456789abcdef0123456789ABCDEF";
	    s0 = ++s;
	    adj = 0;
	    aadj = 1.0;
	    nd0 = -4;

	    if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0;
	    if (*s == '0') {
		while (*++s == '0');
		s1 = strchr(hexdigit, *s);
	    }
	    if (s1 != NULL) {
		do {
		    adj += aadj * ((s1 - hexdigit) & 15);
		    nd0 += 4;
		    aadj /= 16;
		} while (*++s && (s1 = strchr(hexdigit, *s)));
	    }

	    if (*s == '.') {
		dsign = 1;
		if (!*++s || !(s1 = strchr(hexdigit, *s))) goto ret0;
		if (nd0 < 0) {
		    while (*s == '0') {
			s++;
			nd0 -= 4;
		    }
		}
		for (; *s && (s1 = strchr(hexdigit, *s)); ++s) {
		    adj += aadj * ((s1 - hexdigit) & 15);
		    if ((aadj /= 16) == 0.0) {
			while (strchr(hexdigit, *++s));
			break;
		    }
		}
	    }
	    else {
		dsign = 0;
	    }

	    if (*s == 'P' || *s == 'p') {
		dsign = 0x2C - *++s; /* +: 2B, -: 2D */
		if (abs(dsign) == 1) s++;
		else dsign = 1;

		nd = 0;
		c = *s;
		if (c < '0' || '9' < c) goto ret0;
		do {
		    nd *= 10;
		    nd += c;
		    nd -= '0';
		    c = *++s;
		    /* Float("0x0."+("0"*267)+"1fp2095") */
		    if (nd + dsign * nd0 > 2095) {
			while ('0' <= c && c <= '9') c = *++s;
			break;
		    }
		} while ('0' <= c && c <= '9');
		nd0 += nd * dsign;
	    }
	    else {
		if (dsign) goto ret0;
	    }
	    dval(rv) = ldexp(adj, nd0);
	    goto ret;
	}
        nz0 = 1;
        while (*++s == '0') ;
        if (!*s)
            goto ret;
    }
    s0 = s;
    y = z = 0;
    for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
        if (nd < 9)
            y = 10*y + c - '0';
        else if (nd < 16)
            z = 10*z + c - '0';
    nd0 = nd;
#ifdef USE_LOCALE
    s1 = localeconv()->decimal_point;
    if (c == *s1) {
        c = '.';
        if (*++s1) {
            s2 = s;
            for (;;) {
                if (*++s2 != *s1) {
                    c = 0;
                    break;
                }
                if (!*++s1) {
                    s = s2;
                    break;
                }
            }
        }
    }
#endif
    if (c == '.') {
        if (!ISDIGIT(s[1]))
            goto dig_done;
        c = *++s;
        if (!nd) {
            for (; c == '0'; c = *++s)
                nz++;
            if (c > '0' && c <= '9') {
                s0 = s;
                nf += nz;
                nz = 0;
                goto have_dig;
            }
            goto dig_done;
        }
        for (; c >= '0' && c <= '9'; c = *++s) {
have_dig:
            nz++;
            if (c -= '0') {
                nf += nz;
                for (i = 1; i < nz; i++)
                    if (nd++ < 9)
                        y *= 10;
                    else if (nd <= DBL_DIG + 1)
                        z *= 10;
                if (nd++ < 9)
                    y = 10*y + c;
                else if (nd <= DBL_DIG + 1)
                    z = 10*z + c;
                nz = 0;
            }
        }
    }
dig_done:
    e = 0;
    if (c == 'e' || c == 'E') {
        if (!nd && !nz && !nz0) {
            goto ret0;
        }
        s00 = s;
        esign = 0;
        switch (c = *++s) {
          case '-':
            esign = 1;
          case '+':
            c = *++s;
        }
        if (c >= '0' && c <= '9') {
            while (c == '0')
                c = *++s;
            if (c > '0' && c <= '9') {
                L = c - '0';
                s1 = s;
                while ((c = *++s) >= '0' && c <= '9')
                    L = 10*L + c - '0';
                if (s - s1 > 8 || L > 19999)
                    /* Avoid confusion from exponents
                     * so large that e might overflow.
                     */
                    e = 19999; /* safe for 16 bit ints */
                else
                    e = (int)L;
                if (esign)
                    e = -e;
            }
            else
                e = 0;
        }
        else
            s = s00;
    }
    if (!nd) {
        if (!nz && !nz0) {
#ifdef INFNAN_CHECK
            /* Check for Nan and Infinity */
            switch (c) {
              case 'i':
              case 'I':
                if (match(&s,"nf")) {
                    --s;
                    if (!match(&s,"inity"))
                        ++s;
                    word0(rv) = 0x7ff00000;
                    word1(rv) = 0;
                    goto ret;
                }
                break;
              case 'n':
              case 'N':
                if (match(&s, "an")) {
                    word0(rv) = NAN_WORD0;
                    word1(rv) = NAN_WORD1;
#ifndef No_Hex_NaN
                    if (*s == '(') /*)*/
                        hexnan(&rv, &s);
#endif
                    goto ret;
                }
            }
#endif /* INFNAN_CHECK */
ret0:
            s = s00;
            sign = 0;
        }
        goto ret;
    }
    e1 = e -= nf;

    /* Now we have nd0 digits, starting at s0, followed by a
     * decimal point, followed by nd-nd0 digits.  The number we're
     * after is the integer represented by those digits times
     * 10**e */

    if (!nd0)
        nd0 = nd;
    k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
    dval(rv) = y;
    if (k > 9) {
#ifdef SET_INEXACT
        if (k > DBL_DIG)
            oldinexact = get_inexact();
#endif
        dval(rv) = tens[k - 9] * dval(rv) + z;
    }
    bd0 = bb = bd = bs = delta = 0;
    if (nd <= DBL_DIG
#ifndef RND_PRODQUOT
#ifndef Honor_FLT_ROUNDS
        && Flt_Rounds == 1
#endif
#endif
    ) {
        if (!e)
            goto ret;
        if (e > 0) {
            if (e <= Ten_pmax) {
#ifdef VAX
                goto vax_ovfl_check;
#else
#ifdef Honor_FLT_ROUNDS
                /* round correctly FLT_ROUNDS = 2 or 3 */
                if (sign) {
                    dval(rv) = -dval(rv);
                    sign = 0;
                }
#endif
                /* rv = */ rounded_product(dval(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.
                 */
#ifdef Honor_FLT_ROUNDS
                /* round correctly FLT_ROUNDS = 2 or 3 */
                if (sign) {
                    dval(rv) = -dval(rv);
                    sign = 0;
                }
#endif
                e -= i;
                dval(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(dval(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(dval(rv), tens[e]);
#endif
                goto ret;
            }
        }
#ifndef Inaccurate_Divide
        else if (e >= -Ten_pmax) {
#ifdef Honor_FLT_ROUNDS
            /* round correctly FLT_ROUNDS = 2 or 3 */
            if (sign) {
                dval(rv) = -dval(rv);
                sign = 0;
            }
#endif
            /* rv = */ rounded_quotient(dval(rv), tens[-e]);
            goto ret;
        }
#endif
    }
    e1 += nd - k;

#ifdef IEEE_Arith
#ifdef SET_INEXACT
    inexact = 1;
    if (k <= DBL_DIG)
        oldinexact = get_inexact();
#endif
#ifdef Avoid_Underflow
    scale = 0;
#endif
#ifdef Honor_FLT_ROUNDS
    if ((rounding = Flt_Rounds) >= 2) {
        if (sign)
            rounding = rounding == 2 ? 0 : 2;
        else
            if (rounding != 2)
                rounding = 0;
    }
#endif
#endif /*IEEE_Arith*/

    /* Get starting approximation = rv * 10**e1 */

    if (e1 > 0) {
        if ((i = e1 & 15) != 0)
            dval(rv) *= tens[i];
        if (e1 &= ~15) {
            if (e1 > DBL_MAX_10_EXP) {
ovfl:
#ifndef NO_ERRNO
                errno = ERANGE;
#endif
                /* Can't trust HUGE_VAL */
#ifdef IEEE_Arith
#ifdef Honor_FLT_ROUNDS
                switch (rounding) {
                  case 0: /* toward 0 */
                  case 3: /* toward -infinity */
                    word0(rv) = Big0;
                    word1(rv) = Big1;
                    break;
                  default:
                    word0(rv) = Exp_mask;
                    word1(rv) = 0;
                }
#else /*Honor_FLT_ROUNDS*/
                word0(rv) = Exp_mask;
                word1(rv) = 0;
#endif /*Honor_FLT_ROUNDS*/
#ifdef SET_INEXACT
                /* set overflow bit */
                dval(rv0) = 1e300;
                dval(rv0) *= dval(rv0);
#endif
#else /*IEEE_Arith*/
                word0(rv) = Big0;
                word1(rv) = Big1;
#endif /*IEEE_Arith*/
                if (bd0)
                    goto retfree;
                goto ret;
            }
            e1 >>= 4;
            for (j = 0; e1 > 1; j++, e1 >>= 1)
                if (e1 & 1)
                    dval(rv) *= bigtens[j];
            /* The last multiplication could overflow. */
            word0(rv) -= P*Exp_msk1;
            dval(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) != 0)
            dval(rv) /= tens[i];
        if (e1 >>= 4) {
            if (e1 >= 1 << n_bigtens)
                goto undfl;
#ifdef Avoid_Underflow
            if (e1 & Scale_Bit)
                scale = 2*P;
            for (j = 0; e1 > 0; j++, e1 >>= 1)
                if (e1 & 1)
                    dval(rv) *= tinytens[j];
            if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
                    >> Exp_shift)) > 0) {
                /* scaled rv is denormal; zap j low bits */
                if (j >= 32) {
                    word1(rv) = 0;
                    if (j >= 53)
                        word0(rv) = (P+2)*Exp_msk1;
                    else
                        word0(rv) &= 0xffffffff << (j-32);
                }
                else
                    word1(rv) &= 0xffffffff << j;
            }
#else
            for (j = 0; e1 > 1; j++, e1 >>= 1)
                if (e1 & 1)
                    dval(rv) *= tinytens[j];
            /* The last multiplication could underflow. */
            dval(rv0) = dval(rv);
            dval(rv) *= tinytens[j];
            if (!dval(rv)) {
                dval(rv) = 2.*dval(rv0);
                dval(rv) *= tinytens[j];
#endif
                if (!dval(rv)) {
undfl:
                    dval(rv) = 0.;
#ifndef NO_ERRNO
                    errno = ERANGE;
#endif
                    if (bd0)
                        goto retfree;
                    goto ret;
                }
#ifndef Avoid_Underflow
                word0(rv) = Tiny0;
                word1(rv) = Tiny1;
                /* The refinement below will clean
                 * this approximation up.
                 */
            }
#endif
        }
    }

    /* 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(dval(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 Honor_FLT_ROUNDS
        if (rounding != 1)
            bs2++;
#endif
#ifdef Avoid_Underflow
        j = bbe - scale;
        i = j + bbbits - 1; /* logb(rv) */
        if (i < Emin)   /* denormal */
            j += P - Emin;
        else
            j = P + 1 - bbbits;
#else /*Avoid_Underflow*/
#ifdef Sudden_Underflow
#ifdef IBM
        j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
#else
        j = P + 1 - bbbits;
#endif
#else /*Sudden_Underflow*/
        j = bbe;
        i = j + bbbits - 1; /* logb(rv) */
        if (i < Emin)   /* denormal */
            j += P - Emin;
        else
            j = P + 1 - bbbits;
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
        bb2 += j;
        bd2 += j;
#ifdef Avoid_Underflow
        bd2 += scale;
#endif
        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);
#ifdef Honor_FLT_ROUNDS
        if (rounding != 1) {
            if (i < 0) {
                /* Error is less than an ulp */
                if (!delta->x[0] && delta->wds <= 1) {
                    /* exact */
#ifdef SET_INEXACT
                    inexact = 0;
#endif
                    break;
                }
                if (rounding) {
                    if (dsign) {
                        adj = 1.;
                        goto apply_adj;
                    }
                }
                else if (!dsign) {
                    adj = -1.;
                    if (!word1(rv)
                     && !(word0(rv) & Frac_mask)) {
                        y = word0(rv) & Exp_mask;
#ifdef Avoid_Underflow
                        if (!scale || y > 2*P*Exp_msk1)
#else
                        if (y)
#endif
                        {
                            delta = lshift(delta,Log2P);
                            if (cmp(delta, bs) <= 0)
                                adj = -0.5;
                        }
                    }
apply_adj:
#ifdef Avoid_Underflow
                    if (scale && (y = word0(rv) & Exp_mask)
                            <= 2*P*Exp_msk1)
                        word0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
                    if ((word0(rv) & Exp_mask) <=
                            P*Exp_msk1) {
                        word0(rv) += P*Exp_msk1;
                        dval(rv) += adj*ulp(dval(rv));
                        word0(rv) -= P*Exp_msk1;
                    }
                    else
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
                    dval(rv) += adj*ulp(dval(rv));
                }
                break;
            }
            adj = ratio(delta, bs);
            if (adj < 1.)
                adj = 1.;
            if (adj <= 0x7ffffffe) {
                /* adj = rounding ? ceil(adj) : floor(adj); */
                y = adj;
                if (y != adj) {
                    if (!((rounding>>1) ^ dsign))
                        y++;
                    adj = y;
                }
            }
#ifdef Avoid_Underflow
            if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
                word0(adj) += (2*P+1)*Exp_msk1 - y;
#else
#ifdef Sudden_Underflow
            if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
                word0(rv) += P*Exp_msk1;
                adj *= ulp(dval(rv));
                if (dsign)
                    dval(rv) += adj;
                else
                    dval(rv) -= adj;
                word0(rv) -= P*Exp_msk1;
                goto cont;
            }
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
            adj *= ulp(dval(rv));
            if (dsign)
                dval(rv) += adj;
            else
                dval(rv) -= adj;
            goto cont;
        }
#endif /*Honor_FLT_ROUNDS*/

        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
#ifdef IEEE_Arith
#ifdef Avoid_Underflow
                || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
#else
                || (word0(rv) & Exp_mask) <= Exp_msk1
#endif
#endif
            ) {
#ifdef SET_INEXACT
                if (!delta->x[0] && delta->wds <= 1)
                    inexact = 0;
#endif
                break;
            }
            if (!delta->x[0] && delta->wds <= 1) {
                /* exact result */
#ifdef SET_INEXACT
                inexact = 0;
#endif
                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) == (
#ifdef Avoid_Underflow
                        (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
                        ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
#endif
                        0xffffffff)) {
                    /*boundary case -- increment exponent*/
                    word0(rv) = (word0(rv) & Exp_mask)
                                + Exp_msk1
#ifdef IBM
                                | Exp_msk1 >> 4
#endif
                    ;
                    word1(rv) = 0;
#ifdef Avoid_Underflow
                    dsign = 0;
#endif
                    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
#ifdef Avoid_Underflow
                if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
#else
                if (L <= Exp_msk1)
#endif /*Avoid_Underflow*/
#endif /*IBM*/
                    goto undfl;
                L -= Exp_msk1;
#else /*Sudden_Underflow}{*/
#ifdef Avoid_Underflow
                if (scale) {
                    L = word0(rv) & Exp_mask;
                    if (L <= (2*P+1)*Exp_msk1) {
                        if (L > (P+2)*Exp_msk1)
                            /* round even ==> */
                            /* accept rv */
                            break;
                        /* rv = smallest denormal */
                        goto undfl;
                    }
                }
#endif /*Avoid_Underflow*/
                L = (word0(rv) & Exp_mask) - Exp_msk1;
#endif /*Sudden_Underflow}}*/
                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)
                dval(rv) += ulp(dval(rv));
#ifndef ROUND_BIASED
            else {
                dval(rv) -= ulp(dval(rv));
#ifndef Sudden_Underflow
                if (!dval(rv))
                    goto undfl;
#endif
            }
#ifdef Avoid_Underflow
            dsign = 1 - dsign;
#endif
#endif
            break;
        }
        if ((aadj = ratio(delta, bs)) <= 2.) {
            if (dsign)
                aadj = dval(aadj1) = 1.;
            else if (word1(rv) || word0(rv) & Bndry_mask) {
#ifndef Sudden_Underflow
                if (word1(rv) == Tiny1 && !word0(rv))
                    goto undfl;
#endif
                aadj = 1.;
                dval(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;
                dval(aadj1) = -aadj;
            }
        }
        else {
            aadj *= 0.5;
            dval(aadj1) = dsign ? aadj : -aadj;
#ifdef Check_FLT_ROUNDS
            switch (Rounding) {
              case 2: /* towards +infinity */
                dval(aadj1) -= 0.5;
                break;
              case 0: /* towards 0 */
              case 3: /* towards -infinity */
                dval(aadj1) += 0.5;
            }
#else
            if (Flt_Rounds == 0)
                dval(aadj1) += 0.5;
#endif /*Check_FLT_ROUNDS*/
        }
        y = word0(rv) & Exp_mask;

        /* Check for overflow */

        if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
            dval(rv0) = dval(rv);
            word0(rv) -= P*Exp_msk1;
            adj = dval(aadj1) * ulp(dval(rv));
            dval(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 Avoid_Underflow
            if (scale && y <= 2*P*Exp_msk1) {
                if (aadj <= 0x7fffffff) {
                    if ((z = (int)aadj) <= 0)
                        z = 1;
                    aadj = z;
                    dval(aadj1) = dsign ? aadj : -aadj;
                }
                word0(aadj1) += (2*P+1)*Exp_msk1 - y;
            }
            adj = dval(aadj1) * ulp(dval(rv));
            dval(rv) += adj;
#else
#ifdef Sudden_Underflow
            if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
                dval(rv0) = dval(rv);
                word0(rv) += P*Exp_msk1;
                adj = dval(aadj1) * ulp(dval(rv));
                dval(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 = dval(aadj1) * ulp(dval(rv));
                dval(rv) += adj;
            }
#else /*Sudden_Underflow*/
            /* 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.) {
                dval(aadj1) = (double)(int)(aadj + 0.5);
                if (!dsign)
                    dval(aadj1) = -dval(aadj1);
            }
            adj = dval(aadj1) * ulp(dval(rv));
            dval(rv) += adj;
#endif /*Sudden_Underflow*/
#endif /*Avoid_Underflow*/
        }
        z = word0(rv) & Exp_mask;
#ifndef SET_INEXACT
#ifdef Avoid_Underflow
        if (!scale)
#endif
        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;
        }
#endif
cont:
        Bfree(bb);
        Bfree(bd);
        Bfree(bs);
        Bfree(delta);
    }
#ifdef SET_INEXACT
    if (inexact) {
        if (!oldinexact) {
            word0(rv0) = Exp_1 + (70 << Exp_shift);
            word1(rv0) = 0;
            dval(rv0) += 1.;
        }
    }
    else if (!oldinexact)
        clear_inexact();
#endif
#ifdef Avoid_Underflow
    if (scale) {
        word0(rv0) = Exp_1 - 2*P*Exp_msk1;
        word1(rv0) = 0;
        dval(rv) *= dval(rv0);
#ifndef NO_ERRNO
        /* try to avoid the bug of testing an 8087 register value */
        if (word0(rv) == 0 && word1(rv) == 0)
            errno = ERANGE;
#endif
    }
#endif /* Avoid_Underflow */
#ifdef SET_INEXACT
    if (inexact && !(word0(rv) & Exp_mask)) {
        /* set underflow bit */
        dval(rv0) = 1e-300;
        dval(rv0) *= dval(rv0);
    }
#endif
retfree:
    Bfree(bb);
    Bfree(bd);
    Bfree(bs);
    Bfree(bd0);
    Bfree(delta);
ret:
    if (se)
        *se = (char *)s;
    return sign ? -dval(rv) : dval(rv);