port/fp/decimal_bin.c

/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License (the "License").
 * You may not use this file except in compliance with the License.
 *
 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
 * or http://www.opensolaris.org/os/licensing.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 *
 * When distributing Covered Code, include this CDDL HEADER in each
 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
 * If applicable, add the following below this CDDL HEADER, with the
 * fields enclosed by brackets "[]" replaced with your own identifying
 * information: Portions Copyright [yyyy] [name of copyright owner]
 *
 * CDDL HEADER END
 */

/*
 * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
 * Use is subject to license terms.
 */

#pragma ident   "%Z%%M% %I% %E% SMI"

/*
 * Conversion from decimal to binary floating point
 */

#include "lint.h"
#include <stdlib.h>
#include "base_conversion.h"

/*
 * Convert the integer part of a nonzero base-10^4 _big_float *pd
 * to base 2^16 in **ppb.  The converted value is accurate to nsig
 * significant bits.  On exit, *sticky is nonzero if *pd had a
 * nonzero fractional part.  If pd->exponent > 0 and **ppb is not
 * large enough to hold the final converted value (i.e., the con-
 * verted significand scaled by 10^pd->exponent), then on exit,
 * *ppb will point to a newly allocated _big_float, which must be
 * freed by the caller.  (The number of significant bits we need
 * should fit in pb, but __big_float_times_power may allocate new
 * storage anyway because the exact product could require more than
 * 16000 bits.)
 *
 * This routine does not check that **ppb is large enough to hold
 * the result of converting the significand of *pd.
 */
static void
__big_decimal_to_big_binary(_big_float *pd, int nsig, _big_float **ppb,
    int *sticky)
{
    _big_float  *pb;
    int     i, j, len, s;
    unsigned int    carry;

    pb = *ppb;

    /* convert pd a digit at a time, most significant first */
    if (pd->bexponent + ((pd->blength - 1) << 2) >= 0) {
        pb->bsignificand[0] = pd->bsignificand[pd->blength - 1];
        len = 1;
        for (i = pd->blength - 2; i >= 0 &&
            pd->bexponent + (i << 2) >= 0; i--) {
            /* multiply pb by 10^4 and add next digit */
            carry = pd->bsignificand[i];
            for (j = 0; j < len; j++) {
                carry += (unsigned int)pb->bsignificand[j]
                    * 10000;
                pb->bsignificand[j] = carry & 0xffff;
                carry >>= 16;
            }
            if (carry)
                pb->bsignificand[j++] = carry;
            len = j;
        }
    } else {
        i = pd->blength - 1;
        len = 0;
    }

    /* convert any partial digit */
    if (i >= 0 && pd->bexponent + (i << 2) > -4) {
        s = pd->bexponent + (i << 2) + 4;
        /* multiply pb by 10^s and add partial digit */
        carry = pd->bsignificand[i];
        if (s == 1) {
            s = carry % 1000;
            carry = carry / 1000;
            for (j = 0; j < len; j++) {
                carry += (unsigned int)pb->bsignificand[j]
                    * 10;
                pb->bsignificand[j] = carry & 0xffff;
                carry >>= 16;
            }
        } else if (s == 2) {
            s = carry % 100;
            carry = carry / 100;
            for (j = 0; j < len; j++) {
                carry += (unsigned int)pb->bsignificand[j]
                    * 100;
                pb->bsignificand[j] = carry & 0xffff;
                carry >>= 16;
            }
        } else {
            s = carry % 10;
            carry = carry / 10;
            for (j = 0; j < len; j++) {
                carry += (unsigned int)pb->bsignificand[j]
                    * 1000;
                pb->bsignificand[j] = carry & 0xffff;
                carry >>= 16;
            }
        }
        if (carry)
            pb->bsignificand[j++] = carry;
        len = j;
        i--;
    } else {
        s = 0;
    }

    pb->blength = len;
    pb->bexponent = 0;

    /* continue accumulating sticky flag */
    while (i >= 0)
        s |= pd->bsignificand[i--];
    *sticky = s;

    if (pd->bexponent > 0) {
        /* scale pb by 10^pd->exponent */
        __big_float_times_power(pb, 10, pd->bexponent, nsig, ppb);
    }
}

/*
 * Convert the decimal_record *pd to an unpacked datum *px accurately
 * enough that *px can be rounded correctly to sigbits significant bits.
 * (We may assume sigbits <= 113.)
 */
static void
__decimal_to_unpacked(unpacked *px, decimal_record *pd, int sigbits)
{
    _big_float  d, b, *pbd, *pbb;
    char        *ds;
    int     ids, i, ix, exp, ndigs;
    int     sticky, powtwo, sigdigits;

    px->sign = pd->sign;
    px->fpclass = pd->fpclass;
    ds = pd->ds;
    ndigs = pd->ndigits;
    exp = pd->exponent;

    /* remove trailing zeroes */
    while (ndigs > 0 && ds[ndigs - 1] == '0') {
        exp++;
        ndigs--;
    }
    if (ndigs < 1) {
        /* nothing left */
        px->fpclass = fp_zero;
        return;
    }

    /* convert remaining digits to a base-10^4 _big_float */
    d.bsize = _BIG_FLOAT_SIZE;
    d.bexponent = exp;
    d.blength = (ndigs + 3) >> 2;
    i = d.blength - 1;
    ids = ndigs - (d.blength << 2);
    switch (ids) {
    case -1:
        d.bsignificand[i] = 100 * ds[ids + 1] +
            10 * ds[ids + 2] + ds[ids + 3] - 111 * '0';
        i--;
        ids += 4;
        break;

    case -2:
        d.bsignificand[i] = 10 * ds[ids + 2] + ds[ids + 3] - 11 * '0';
        i--;
        ids += 4;
        break;

    case -3:
        d.bsignificand[i] = ds[ids + 3] - '0';
        i--;
        ids += 4;
        break;
    }
    while (i >= 0) {
        d.bsignificand[i] = 1000 * ds[ids] + 100 * ds[ids + 1] +
            10 * ds[ids + 2] + ds[ids + 3] - 1111 * '0';
        i--;
        ids += 4;
    }

    pbd = &d;
    powtwo = 0;

    /* pre-scale to get the bits we want into the integer part */
    if (exp < 0) {
        /* i is a lower bound on log10(x) */
        i = exp + ndigs - 1;
        if (i <= 0 || ((i * 217705) >> 16) < sigbits + 2) {
            /*
             * Scale by 2^(sigbits + 2 + u) where
             * u is an upper bound on -log2(x).
             */
            powtwo = sigbits + 2;
            if (i < 0)
                powtwo += ((-i * 217706) + 65535) >> 16;
            else if (i > 0)
                powtwo -= (i * 217705) >> 16;
            /*
             * Take sigdigits large enough to get
             * all integral digits correct.
             */
            sigdigits = i + 1 + (((powtwo * 19729) + 65535) >> 16);
            __big_float_times_power(&d, 2, powtwo, sigdigits, &pbd);
        }
    }

    /* convert to base 2^16 */
    b.bsize = _BIG_FLOAT_SIZE;
    pbb = &b;
    __big_decimal_to_big_binary(pbd, sigbits + 2, &pbb, &sticky);

    /* adjust pbb->bexponent based on the scale factor above */
    pbb->bexponent -= powtwo;

    /* convert to unpacked */
    ix = 0;
    for (i = pbb->blength - 1; i > 0 && ix < 5; i -= 2) {
        px->significand[ix++] = (pbb->bsignificand[i] << 16) |
            pbb->bsignificand[i - 1];
    }
    if (ix < 5) {
        /* pad with zeroes */
        if (i == 0)
            px->significand[ix++] = pbb->bsignificand[i] << 16;
        while (ix < 5)
            px->significand[ix++] = 0;
    } else {
        /* truncate and set a sticky bit if necessary */
        while (i >= 0 && pbb->bsignificand[i] == 0)
            i--;
        if (i >= 0)
            px->significand[4] |= 1;
    }
    if (sticky | pd->more)
        px->significand[4] |= 1;
    px->exponent = pbb->bexponent + (pbb->blength << 4) - 1;

    /* normalize so the most significant bit is set */
    while (px->significand[0] < 0x80000000u) {
        px->significand[0] = (px->significand[0] << 1) |
            (px->significand[1] >> 31);
        px->significand[1] = (px->significand[1] << 1) |
            (px->significand[2] >> 31);
        px->significand[2] = (px->significand[2] << 1) |
            (px->significand[3] >> 31);
        px->significand[3] = (px->significand[3] << 1) |
            (px->significand[4] >> 31);
        px->significand[4] <<= 1;
        px->exponent--;
    }

    if (pbd != &d)
        (void) free((void *)pbd);
    if (pbb != &b)
        (void) free((void *)pbb);
}

/*
 * Convert a string s consisting of n <= 18 ASCII decimal digits
 * to an integer value in double precision format, and set *pe
 * to the number of rounding errors incurred (0 or 1).
 */
static double
__digits_to_double(char *s, int n, int *pe)
{
    int i, acc;
    double  t, th, tl;

    if (n <= 9) {
        acc = s[0] - '0';
        for (i = 1; i < n; i++) {
            /* acc <- 10 * acc + next digit */
            acc = (acc << 1) + (acc << 3) + s[i] - '0';
        }
        t = (double)acc;
        *pe = 0;
    } else {
        acc = s[0] - '0';
        for (i = 1; i < (n - 9); i++) {
            /* acc <- 10 * acc + next digit */
            acc = (acc << 1) + (acc << 3) + s[i] - '0';
        }
        th = 1.0e9 * (double)acc;   /* this will be exact */
        acc = s[n - 9] - '0';
        for (i = n - 8; i < n; i++) {
            /* acc <- 10 * acc + next digit */
            acc = (acc << 1) + (acc << 3) + s[i] - '0';
        }
        tl = (double)acc;
        /* add and indicate whether or not the sum is exact */
        t = th + tl;
        *pe = ((t - th) == tl)? 0 : 1;
    }
    return (t);
}

static union {
    int i[2];
    double  d;
} C[] = {
#ifdef _LITTLE_ENDIAN
    { 0x00000000, 0x3cc40000 },
#else
    { 0x3cc40000, 0x00000000 }, /* 5 * 2^-53 */
#endif
};

#define five2m53    C[0].d

static int
__fast_decimal_to_single(single *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    double          dds, delta, ddsplus, ddsminus, df1;
    int         n, exp, rounded, e;
    float           f1, f2;
    __ieee_flags_type   fb;

    if (pm->rd != fp_nearest)
        return (0);

    exp = pd->exponent;
    if (pd->ndigits <= 18) {
        rounded = 0;
        n = pd->ndigits;
    } else {
        rounded = 1;
        n = 18;
        exp += pd->ndigits - 18;
    }
    /*
     * exp must be in the range of the table, and the result
     * must not underflow or overflow.
     */
    if (exp < -__TBL_TENS_MAX || exp + n < -36 || exp + n > 38)
        return (0);

    __get_ieee_flags(&fb);
    dds = __digits_to_double(pd->ds, n, &e);
    if (e != 0)
        rounded = 1;
    if (exp > 0) {
        /* small positive exponent */
        if (exp > __TBL_TENS_EXACT)
            rounded = 1;
        if (rounded) {
            dds *= __tbl_tens[exp];
        } else {
            dds = __mul_set(dds, __tbl_tens[exp], &e);
            if (e)
                rounded = 1;
        }
    } else if (exp < 0) {
        /* small negative exponent */
        if (-exp > __TBL_TENS_EXACT)
            rounded = 1;
        if (rounded) {
            dds /= __tbl_tens[-exp];
        } else {
            dds = __div_set(dds, __tbl_tens[-exp], &e);
            if (e)
                rounded = 1;
        }
    }

    /*
     * At this point dds may have four rounding errors due to
     * (i) truncation of pd->ds to 18 digits, (ii) inexact con-
     * version of pd->ds to binary, (iii) scaling by a power of
     * ten that is not exactly representable, and (iv) roundoff
     * error in the multiplication.  Below we will incur one more
     * rounding error when we add or subtract delta and dds.  We
     * construct delta so that even after this last rounding error,
     * ddsplus is an upper bound on the exact value and ddsminus
     * is a lower bound.  Then as long as both of these quantities
     * round to the same single precision number, that number
     * will be the correctly rounded single precision result.
     * (If any rounding errors have been committed, then we must
     * also be certain that the result can't be exact.)
     */
    delta = five2m53 * dds;
    ddsplus = dds + delta;
    ddsminus = dds - delta;
    f1 = (float)(ddsplus);
    f2 = (float)(ddsminus);
    df1 = f1;
    __set_ieee_flags(&fb);
    if (f1 != f2)
        return (0);
    if (rounded) {
        /*
         * If ddsminus <= f1 <= ddsplus, the result might be
         * exact; we have to convert the long way to be sure.
         */
        if (ddsminus <= df1 && df1 <= ddsplus)
            return (0);
        *ps = (1 << fp_inexact);
    } else {
        *ps = (df1 == dds)? 0 : (1 << fp_inexact);
    }
    *px = (pd->sign)? -f1 : f1;
    return (1);
}

/*
 * Attempt conversion to double using floating-point arithmetic.
 * Return 1 if it works (at most one rounding error), 0 if it doesn't.
 */
static int
__fast_decimal_to_double(double *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    double          dds;
    int         e;
    __ieee_flags_type   fb;

    if (pm->rd != fp_nearest || pd->ndigits > 18 || pd->exponent
        > __TBL_TENS_EXACT || pd->exponent < -__TBL_TENS_EXACT)
        return (0);

    __get_ieee_flags(&fb);
    dds = __digits_to_double(pd->ds, pd->ndigits, &e);
    if (e != 0) {
        __set_ieee_flags(&fb);
        return (0);
    }
    if (pd->exponent > 0)
        dds = __mul_set(dds, __tbl_tens[pd->exponent], &e);
    else if (pd->exponent < 0)
        dds = __div_set(dds, __tbl_tens[-pd->exponent], &e);
    *px = (pd->sign)? -dds : dds;
    *ps = (e)? (1 << fp_inexact) : 0;
    __set_ieee_flags(&fb);
    return (1);
}

/* PUBLIC FUNCTIONS */

/*
 * The following routines convert the decimal record *pd to a floating
 * point value *px observing the rounding mode specified in pm->rd and
 * passing back any floating point exceptions that occur in *ps.
 *
 * pd->sign and pd->fpclass are always taken into account.  pd->exponent
 * and pd->ds are used when pd->fpclass is fp_normal or fp_subnormal.
 * In these cases, pd->ds must contain a null-terminated string of one
 * or more ASCII digits, the first of which is not zero, and pd->ndigits
 * must equal the length of this string.  If m is the integer represented
 * by the string pd->ds, then *px will be set to a correctly rounded
 * approximation to
 *
 *   -1**(pd->sign) * m * 10**(pd->exponent)
 *
 * (If pd->more != 0 then additional nonzero digits are assumed to follow
 * those in pd->ds, so m is effectively replaced by m + epsilon in the
 * expression above.)
 *
 * For example, if pd->exponent == -2 and pd->ds holds "1234", then *px
 * will be a correctly rounded approximation to 12.34.
 *
 * Note that the only mode that matters is the rounding direction pm->rd;
 * pm->df and pm->ndigits are never used.
 */

/* maximum decimal exponent we need to consider */
#define SINGLE_MAXE   47
#define DOUBLE_MAXE  326
#define EXTENDED_MAXE   4968
#define QUAD_MAXE   4968

void
decimal_to_single(single *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    single_equivalence  *kluge;
    unpacked        u;
    fp_exception_field_type ef;
    int         i;

    /* special values */
    kluge = (single_equivalence *)px;
    switch (pd->fpclass) {
    case fp_zero:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0;
        kluge->f.msw.significand = 0;
        *ps = 0;
        return;

    case fp_infinity:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0xff;
        kluge->f.msw.significand = 0;
        *ps = 0;
        return;

    case fp_quiet:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0xff;
        kluge->f.msw.significand = 0x7fffff;
        *ps = 0;
        return;

    case fp_signaling:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0xff;
        kluge->f.msw.significand = 0x3fffff;
        *ps = 0;
        return;
    }

    /* numeric values */
    ef = 0;
    if (pd->exponent + pd->ndigits > SINGLE_MAXE) {
        /* result must overflow */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = 0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else if (pd->exponent + pd->ndigits < -SINGLE_MAXE) {
        /* result must underflow completely */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = -0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else {
        /* result may be in range */
        if (__fast_decimal_to_single(px, pm, pd, &ef) == 1) {
            *ps = ef;
            if (ef != 0)
                __base_conversion_set_exception(ef);
            return;
        }
        __decimal_to_unpacked(&u, pd, 24);
    }
    __pack_single(&u, px, pm->rd, &ef);
    *ps = ef;
    if (ef != 0)
        __base_conversion_set_exception(ef);
}

void
decimal_to_double(double *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    double_equivalence  *kluge;
    unpacked        u;
    fp_exception_field_type ef;
    int         i;

    /* special values */
    kluge = (double_equivalence *)px;
    switch (pd->fpclass) {
    case fp_zero:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0;
        kluge->f.msw.significand = 0;
        kluge->f.significand2 = 0;
        *ps = 0;
        return;

    case fp_infinity:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7ff;
        kluge->f.msw.significand = 0;
        kluge->f.significand2 = 0;
        *ps = 0;
        return;

    case fp_quiet:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7ff;
        kluge->f.msw.significand = 0xfffff;
        kluge->f.significand2 = 0xffffffff;
        *ps = 0;
        return;

    case fp_signaling:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7ff;
        kluge->f.msw.significand = 0x7ffff;
        kluge->f.significand2 = 0xffffffff;
        *ps = 0;
        return;
    }

    /* numeric values */
    ef = 0;
    if (pd->exponent + pd->ndigits > DOUBLE_MAXE) {
        /* result must overflow */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = 0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else if (pd->exponent + pd->ndigits < -DOUBLE_MAXE) {
        /* result must underflow completely */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = -0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else {
        /* result may be in range */
        if (__fast_decimal_to_double(px, pm, pd, &ef) == 1) {
            *ps = ef;
            if (ef != 0)
                __base_conversion_set_exception(ef);
            return;
        }
        __decimal_to_unpacked(&u, pd, 53);
    }
    __pack_double(&u, px, pm->rd, &ef);
    *ps = ef;
    if (ef != 0)
        __base_conversion_set_exception(ef);
}

void
decimal_to_extended(extended *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    extended_equivalence    *kluge;
    unpacked        u;
    double_equivalence  dd;
    fp_exception_field_type ef;
    int         i;

    /* special values */
    kluge = (extended_equivalence *)px;
    switch (pd->fpclass) {
    case fp_zero:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0;
        kluge->f.significand = 0;
        kluge->f.significand2 = 0;
        *ps = 0;
        return;

    case fp_infinity:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.significand = 0x80000000;
        kluge->f.significand2 = 0;
        *ps = 0;
        return;

    case fp_quiet:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.significand = 0xffffffff;
        kluge->f.significand2 = 0xffffffff;
        *ps = 0;
        return;

    case fp_signaling:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.significand = 0xbfffffff;
        kluge->f.significand2 = 0xffffffff;
        *ps = 0;
        return;
    }

    /* numeric values */
    ef = 0;
    if (pd->exponent + pd->ndigits > EXTENDED_MAXE) {
        /* result must overflow */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = 0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else if (pd->exponent + pd->ndigits < -EXTENDED_MAXE) {
        /* result must underflow completely */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = -0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else {
        /* result may be in range */
        if (__fast_decimal_to_double(&dd.x, pm, pd, &ef) == 1 &&
            ef == 0) {
            u.sign = dd.f.msw.sign;
            u.fpclass = fp_normal;
            u.exponent = dd.f.msw.exponent - DOUBLE_BIAS;
            u.significand[0] = ((0x100000 |
                dd.f.msw.significand) << 11) |
                (dd.f.significand2 >> 21);
            u.significand[1] = dd.f.significand2 << 11;
            for (i = 2; i < UNPACKED_SIZE; i++)
                u.significand[i] = 0;
        } else {
            __decimal_to_unpacked(&u, pd, 64);
        }
    }
    __pack_extended(&u, px, pm->rd, &ef);
    *ps = ef;
    if (ef != 0)
        __base_conversion_set_exception(ef);
}

void
decimal_to_quadruple(quadruple *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
    quadruple_equivalence   *kluge;
    unpacked        u;
    double_equivalence  dd;
    fp_exception_field_type ef;
    int         i;

    /* special values */
    kluge = (quadruple_equivalence *)px;
    switch (pd->fpclass) {
    case fp_zero:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0;
        kluge->f.msw.significand = 0;
        kluge->f.significand2 = 0;
        kluge->f.significand3 = 0;
        kluge->f.significand4 = 0;
        *ps = 0;
        return;

    case fp_infinity:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.msw.significand = 0;
        kluge->f.significand2 = 0;
        kluge->f.significand3 = 0;
        kluge->f.significand4 = 0;
        *ps = 0;
        return;

    case fp_quiet:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.msw.significand = 0xffff;
        kluge->f.significand2 = 0xffffffff;
        kluge->f.significand3 = 0xffffffff;
        kluge->f.significand4 = 0xffffffff;
        *ps = 0;
        return;

    case fp_signaling:
        kluge->f.msw.sign = (pd->sign)? 1 : 0;
        kluge->f.msw.exponent = 0x7fff;
        kluge->f.msw.significand = 0x7fff;
        kluge->f.significand2 = 0xffffffff;
        kluge->f.significand3 = 0xffffffff;
        kluge->f.significand4 = 0xffffffff;
        *ps = 0;
        return;
    }

    /* numeric values */
    ef = 0;
    if (pd->exponent + pd->ndigits > QUAD_MAXE) {
        /* result must overflow */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = 0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else if (pd->exponent + pd->ndigits < -QUAD_MAXE) {
        /* result must underflow completely */
        u.sign = (pd->sign != 0);
        u.fpclass = fp_normal;
        u.exponent = -0x000fffff;
        u.significand[0] = 0x80000000;
        for (i = 1; i < UNPACKED_SIZE; i++)
            u.significand[i] = 0;
    } else {
        /* result may be in range */
        if (__fast_decimal_to_double(&dd.x, pm, pd, &ef) == 1 &&
            ef == 0) {
            u.sign = dd.f.msw.sign;
            u.fpclass = fp_normal;
            u.exponent = dd.f.msw.exponent - DOUBLE_BIAS;
            u.significand[0] = ((0x100000 |
                dd.f.msw.significand) << 11) |
                (dd.f.significand2 >> 21);
            u.significand[1] = dd.f.significand2 << 11;
            for (i = 2; i < UNPACKED_SIZE; i++)
                u.significand[i] = 0;
        } else {
            __decimal_to_unpacked(&u, pd, 113);
        }
    }
    __pack_quadruple(&u, px, pm->rd, &ef);
    *ps = ef;
    if (ef != 0)
        __base_conversion_set_exception(ef);
}