strtod.c File Reference

#include "gdtoaimp.h"
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <fenv.h>

Include dependency graph for strtod.c:

Go to the source code of this file.

Macros
#define	Rounding   Flt_Rounds
#define	fpi1   fpi

Functions
double	strtod (CONST char *s00, char **se)

Macro Definition Documentation

fpi1

#define fpi1   fpi

Rounding

#define Rounding   Flt_Rounds

Definition at line 59 of file strtod.c.

Function Documentation

strtod()

double strtod	(	CONST char *	s00,
		char **	se
	)

Definition at line 67 of file strtod.c.

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

    sign = nz0 = nz = decpt = 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')
    {
#ifndef NO_HEX_FP
        {
            static FPI fpi = {53, 1 - 1023 - 53 + 1, 2046 - 1023 - 53 + 1, 1, SI};
            Long exp;
            ULong bits[2];
            switch(s[1])
            {
                case 'x':
                case 'X':
                {
#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
                    FPI fpi1 = fpi;
                    switch(fegetround())
                    {
                        case FE_TOWARDZERO:
                            fpi1.rounding = 0;
                            break;
                        case FE_UPWARD:
                            fpi1.rounding = 2;
                            break;
                        case FE_DOWNWARD:
                            fpi1.rounding = 3;
                    }
#else
#define fpi1 fpi
#endif
                    switch((i = gethex(&s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask)
                    {
                        case STRTOG_NoNumber:
                            s = s00;
                            sign = 0;
                        case STRTOG_Zero:
                            break;
                        default:
                            if(bb)
                            {
                                copybits(bits, fpi.nbits, bb);
                                Bfree(bb);
                            }
                            ULtod(((U*)&rv)->L, bits, exp, i);
                    }
                }
                    goto ret;
            }
        }
#endif
        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 + (unsigned)c - '0';
        }
        else if(nd < 16)
        {
            z = 10 * z + (unsigned)c - '0';
        }
    }
    nd0 = nd;
#ifdef USE_LOCALE
    if(c == *localeconv()->decimal_point)
#else
    if(c == '.')
#endif
    {
        decpt = 1;
        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 + (unsigned)c;
                }
                else if(nd <= DBL_DIG + 1)
                {
                    z = 10 * z + (unsigned)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 */
            ULong bits[2];
            static FPI fpinan = /* only 52 explicit bits */
                {52, 1 - 1023 - 53 + 1, 2046 - 1023 - 53 + 1, 1, SI};
            if(!decpt)
                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"))
                        {
#ifndef No_Hex_NaN
                            if(*s == '(' /*)*/
                               && hexnan(&s, &fpinan, bits) == STRTOG_NaNbits)
                            {
                                word0(rv) = 0x7ff00000 | bits[1];
                                word1(rv) = bits[0];
                            }
                            else
                            {
#endif
                                word0(rv) = NAN_WORD0;
                                word1(rv) = NAN_WORD1;
#ifndef No_Hex_NaN
                            }
#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 = 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)
                {
                    rv = -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)
                {
                    rv = -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)
            {
                rv = -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(fabs(dval(rv)) <= DBL_EPSILON)
            {
            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 = (ULong)adj;
            if(fabs(y - adj) > DBL_EPSILON)
            {
                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) = (ULong)(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(fabs(dval(rv)) <= DBL_EPSILON)
            {
                goto undfl;
            }
#endif
        }
#ifdef Avoid_Underflow
        dsign = 1 - dsign;
#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(Rounding)
        {
            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 /*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 = 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)
            {
                z = (ULong)aadj;
                if(z <= 0)
                {
                    z = 1;
                }
                aadj = z;
                aadj1 = dsign ? aadj : -aadj;
            }
            word0(aadj1) += (2 * P + 1) * Exp_msk1 - y;
        }
        adj = 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 = 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 = 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.)
                {
                    aadj1 = (double)(int)(aadj + 0.5);
                    if(!dsign)
                    {
                        aadj1 = -aadj1;
                    }
                }
                adj = 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;
            }
        }
#ifdef Avoid_Underflow
    }
#endif
#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 : if(bb)
{
    Bfree(bb);
}

if(bd)
{
    Bfree(bd);
}

if(bs)
{
    Bfree(bs);
}

if(bd0)
{
    Bfree(bd0);
}

if(delta)
{
    Bfree(delta);
}

ret : if(se)
{
    *se = (char*)(uintptr_t)s;
}
return sign ? -dval(rv) : dval(rv);
}

References Balloc(), Bcopy, Bfree(), Bias, Big0, Big1, bigtens, Bndry_mask, Bndry_mask1, cmp(), CONST, copybits(), d2b(), DBL_DIG, DBL_EPSILON, DBL_MAX_10_EXP, DBL_MAX_EXP, lconv::decimal_point, diff(), dval, ERANGE, errno, Exp_1, Exp_mask, Exp_msk1, Exp_shift, fabs(), FLT_RADIX, Flt_Rounds, fpi1, Frac_mask, gethex(), hexnan(), Honor_FLT_ROUNDS, i2b(), IBM, localeconv(), Log2P, Long, LSB, lshift(), match, mult(), FPI::nbits, NULL, P, pow5mult(), ratio(), rounded_product, rounded_quotient, Rounding, s2b(), SI, STRTOG_NaNbits, STRTOG_NoNumber, STRTOG_Retmask, STRTOG_Zero, Ten_pmax, tens, Tiny0, Tiny1, tinytens, ulp(), ULtod(), word0, and word1.