import less(1)

[unleashed/tickless.git] / usr / src / lib / libc / port / fp / double_decim.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.
 */

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

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

/*
 * Any sensible programmer would inline the following routine where
 * it is used below.  Unfortunately, the Sun SPARC compilers are not
 * consistent in generating efficient code for this, so inlining it
 * as written can cause the *_to_decimal functions to take twice as
 * long in some cases.
 *
 * We might be tempted, then, to rewrite the source to match the most
 * efficient code the compilers generate and inline that.  Alas, the
 * most efficient code on SPARC uses 32x32->64 bit multiply, which
 * can't be expressed directly in source code.  We could use long long,
 * which would imply 64x64->64 bit multiply; this would work perfectly
 * well on SPARC in v8plus mode.  But as of Solaris 10, libc for SPARC
 * is still built in v8 mode, and of course, x86 is another story.
 *
 * We could also choose to use an inline template to get the most
 * efficient code without incurring the full cost of a function call.
 * Since I expect that would not buy much performance gain, and I
 * prefer to avoid using inline templates for things that can be
 * written in a perfectly straightforward way in C, I've settled
 * for this implementation.  I hope that someday the compilers will
 * get less flaky and/or someone will come up with a better way to
 * do this.
 */
static unsigned int
__quorem10000(unsigned int x, unsigned short *pr)
{
        *pr = x % 10000;
        return (x / 10000);
}

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

        pd = *ppd;

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

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

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

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

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

/*
 * Convert a base-10^4 _big_float *pf to a decimal string in *pd,
 * rounding according to the modes in *pm and recording any exceptions
 * in *ps.  If sticky is nonzero, then additional nonzero digits are
 * assumed to follow those in *pf.  pd->sign must have already been
 * filled in, and pd->fpclass is not modified.  The resulting string
 * is stored in pd->ds, terminated by a null byte.  The length of this
 * string is stored in pd->ndigits, and the corresponding exponent
 * is stored in pd->exponent.  If the converted value is not exact,
 * the inexact flag is set in *ps.
 *
 * When pm->df == fixed_form, we may discover that the result would
 * have more than DECIMAL_STRING_LENGTH - 1 digits.  In this case,
 * we put DECIMAL_STRING_LENGTH - 1 digits into *pd, adjusting both
 * the exponent and the decimal place at which the value is rounded
 * as need be, and we set the overflow flag in *ps.  (Raising overflow
 * is a bug, but we have to do it to maintain backward compatibility.)
 *
 * *pf may be modified.
 */
static void
__big_decimal_to_string(_big_float *pf, int sticky, decimal_mode *pm,
    decimal_record *pd, fp_exception_field_type *ps)
{
        unsigned short  d;
        int             e, er, efirst, elast, i, is, j;
        char            s[4], round;

        /* set e = floor(log10(*pf)) */
        i = pf->blength - 1;
        if (i < 0) {
                e = pf->bexponent = -DECIMAL_STRING_LENGTH - 2;
        } else {
                e = pf->bexponent + (i << 2);
                d = pf->bsignificand[i];
                if (d >= 1000)
                        e += 3;
                else if (d >= 100)
                        e += 2;
                else if (d >= 10)
                        e++;
        }

        /*
         * Determine the power of ten after which to round and the
         * powers corresponding to the first and last digits desired
         * in the result.
         */
        if (pm->df == fixed_form) {
                /* F format */
                er = -pm->ndigits;
                if (er < 0) {
                        efirst = (e >= 0)? e : -1;
                        elast = er;
                } else {
                        efirst = (e >= er)? e : ((er > 0)? er - 1 : 0);
                        elast = 0;
                }

                /* check for possible overflow of pd->ds */
                if (efirst - elast >= DECIMAL_STRING_LENGTH - 1) {
                        efirst = e;
                        elast = e - DECIMAL_STRING_LENGTH + 2;
                        if (er < elast)
                                er = elast;
                        *ps |= 1 << fp_overflow;
                }
        } else {
                /* E format */
                efirst = e;
                elast = er = e - pm->ndigits + 1;
        }

        /* retrieve digits down to the (er - 1) place */
        is = 0;
        for (e = efirst; e >= pf->bexponent + (pf->blength << 2) &&
            e >= er - 1; e--)
                pd->ds[is++] = '0';

        i = pf->blength - 1;
        j = 3 - ((e - pf->bexponent) & 3);
        if (j > 0 && e >= er - 1) {
                __four_digits_quick(pf->bsignificand[i], s);
                while (j <= 3 && e >= er - 1) {
                        pd->ds[is++] = s[j++];
                        e--;
                }
                while (j <= 3)
                        sticky |= (s[j++] - '0');
                i--;
        }

        while ((i | (e - er - 2)) >= 0) {  /* i >= 0 && e >= er + 2 */
                __four_digits_quick(pf->bsignificand[i], pd->ds + is);
                is += 4;
                e -= 4;
                i--;
        }

        if (i >= 0) {
                if (e >= er - 1) {
                        __four_digits_quick(pf->bsignificand[i], s);
                        for (j = 0; e >= er - 1; j++) {
                                pd->ds[is++] = s[j];
                                e--;
                        }
                        while (j <= 3)
                                sticky |= (s[j++] - '0');
                        i--;
                }
        } else {
                while (e-- >= er - 1)
                        pd->ds[is++] = '0';
        }

        /* collect rounding information */
        round = pd->ds[--is];
        while (i >= 0)
                sticky |= pf->bsignificand[i--];

        /* add more trailing zeroes if need be */
        for (e = er - 1; e >= elast; e--)
                pd->ds[is++] = '0';

        pd->exponent = elast;
        pd->ndigits = is;
        pd->ds[is] = '\0';

        /* round */
        if (round == '0' && sticky == 0)
                return;

        *ps |= 1 << fp_inexact;

        switch (pm->rd) {
        case fp_nearest:
                if (round < '5' || (round == '5' && sticky == 0 &&
                    (is <= 0 || (pd->ds[is - 1] & 1) == 0)))
                        return;
                break;

        case fp_positive:
                if (pd->sign)
                        return;
                break;

        case fp_negative:
                if (!pd->sign)
                        return;
                break;

        default:
                return;
        }

        /* round up */
        for (i = efirst - er; i >= 0 && pd->ds[i] == '9'; i--)
                pd->ds[i] = '0';
        if (i >= 0) {
                (pd->ds[i])++;
        } else {
                /* rounding carry out has occurred */
                pd->ds[0] = '1';
                if (pm->df == floating_form) {
                        pd->exponent++;
                } else if (is == DECIMAL_STRING_LENGTH - 1) {
                        pd->exponent++;
                        *ps |= 1 << fp_overflow;
                } else {
                        if (is > 0)
                                pd->ds[is] = '0';
                        is++;
                        pd->ndigits = is;
                        pd->ds[is] = '\0';
                }
        }
}

/*
 * Convert a binary floating point value represented by *pf to a
 * decimal record *pd according to the modes in *pm.  Any exceptions
 * incurred are passed back via *ps.
 */
static void
__bigfloat_to_decimal(_big_float *bf, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
        _big_float      *pbf, *pbd, d;
        int             sticky, powten, sigbits, sigdigits, i;

        /*
         * If pm->ndigits is too large or too small, set the overflow
         * flag in *ps and do nothing.  (Raising overflow is a bug,
         * but we have to do it to maintain backward compatibility.)
         */
        if (pm->ndigits >= DECIMAL_STRING_LENGTH || pm->ndigits <=
            ((pm->df == floating_form)? 0 : -DECIMAL_STRING_LENGTH)) {
                *ps = 1 << fp_overflow;
                return;
        }

        pbf = bf;
        powten = 0;

        /* pre-scale to get the digits we want into the integer part */
        if (pm->df == fixed_form) {
                /* F format */
                if (pm->ndigits >= 0 && bf->bexponent < 0) {
                        /*
                         * Scale by 10^min(-bf->bexponent, pm->ndigits + 1).
                         */
                        powten = pm->ndigits + 1;
                        if (powten > -bf->bexponent)
                                powten = -bf->bexponent;
                        /*
                         * Take sigbits large enough to get all integral
                         * digits correct.
                         */
                        sigbits = bf->bexponent + (bf->blength << 4) +
                            (((powten * 217706) + 65535) >> 16);
                        if (sigbits < 1)
                                sigbits = 1;
                        __big_float_times_power(bf, 10, powten, sigbits, &pbf);
                }
                sigdigits = DECIMAL_STRING_LENGTH + 1;
        } else {
                /* E format */
                if (bf->bexponent < 0) {
                        /* i is a lower bound on log2(x) */
                        i = bf->bexponent + ((bf->blength - 1) << 4);
                        if (i <= 0 || ((i * 19728) >> 16) < pm->ndigits + 1) {
                                /*
                                 * Scale by 10^min(-bf->bexponent,
                                 * pm->ndigits + 1 + u) where u is
                                 * an upper bound on -log10(x).
                                 */
                                powten = pm->ndigits + 1;
                                if (i < 0)
                                        powten += ((-i * 19729) + 65535) >> 16;
                                else if (i > 0)
                                        powten -= (i * 19728) >> 16;
                                if (powten > -bf->bexponent)
                                        powten = -bf->bexponent;
                                /*
                                 * Take sigbits large enough to get
                                 * all integral digits correct.
                                 */
                                sigbits = i + 16 +
                                    (((powten * 217706) + 65535) >> 16);
                                __big_float_times_power(bf, 10, powten,
                                    sigbits, &pbf);
                        }
                }
                sigdigits = pm->ndigits + 2;
        }

        /* convert to base 10^4 */
        d.bsize = _BIG_FLOAT_SIZE;
        pbd = &d;
        __big_binary_to_big_decimal(pbf, sigdigits, &pbd, &sticky);

        /* adjust pbd->bexponent based on the scale factor above */
        pbd->bexponent -= powten;

        /* convert to ASCII */
        __big_decimal_to_string(pbd, sticky, pm, pd, ps);

        if (pbf != bf)
                (void) free((void *)pbf);
        if (pbd != &d)
                (void) free((void *)pbd);
}

/* remove trailing zeroes from the significand of p */
static void
__shorten(_big_float *p)
{
        int     length = p->blength;
        int     zeros, i;

        /* count trailing zeros */
        for (zeros = 0; zeros < length && p->bsignificand[zeros] == 0; zeros++)
                ;
        if (zeros) {
                length -= zeros;
                p->bexponent += zeros << 4;
                for (i = 0; i < length; i++)
                        p->bsignificand[i] = p->bsignificand[i + zeros];
                p->blength = length;
        }
}

/*
 * Unpack a normal or subnormal double into a _big_float.
 */
static void
__double_to_bigfloat(double *px, _big_float *pf)
{
        double_equivalence      *x;

        x = (double_equivalence *)px;
        pf->bsize = _BIG_FLOAT_SIZE;
        pf->bexponent = x->f.msw.exponent - DOUBLE_BIAS - 52;
        pf->blength = 4;
        pf->bsignificand[0] = x->f.significand2 & 0xffff;
        pf->bsignificand[1] = x->f.significand2 >> 16;
        pf->bsignificand[2] = x->f.msw.significand & 0xffff;
        pf->bsignificand[3] = x->f.msw.significand >> 16;
        if (x->f.msw.exponent == 0) {
                pf->bexponent++;
                while (pf->bsignificand[pf->blength - 1] == 0)
                        pf->blength--;
        } else {
                pf->bsignificand[3] += 0x10;
        }
        __shorten(pf);
}

/*
 * Unpack a normal or subnormal extended into a _big_float.
 */
static void
__extended_to_bigfloat(extended *px, _big_float *pf)
{
        extended_equivalence    *x;

        x = (extended_equivalence *)px;
        pf->bsize = _BIG_FLOAT_SIZE;
        pf->bexponent = x->f.msw.exponent - EXTENDED_BIAS - 63;
        pf->blength = 4;
        pf->bsignificand[0] = x->f.significand2 & 0xffff;
        pf->bsignificand[1] = x->f.significand2 >> 16;
        pf->bsignificand[2] = x->f.significand & 0xffff;
        pf->bsignificand[3] = x->f.significand >> 16;
        if (x->f.msw.exponent == 0) {
                pf->bexponent++;
                while (pf->bsignificand[pf->blength - 1] == 0)
                        pf->blength--;
        }
        __shorten(pf);
}

/*
 * Unpack a normal or subnormal quad into a _big_float.
 */
static void
__quadruple_to_bigfloat(quadruple *px, _big_float *pf)
{
        quadruple_equivalence   *x;

        x = (quadruple_equivalence *)px;
        pf->bsize = _BIG_FLOAT_SIZE;
        pf->bexponent = x->f.msw.exponent - QUAD_BIAS - 112;
        pf->bsignificand[0] = x->f.significand4 & 0xffff;
        pf->bsignificand[1] = x->f.significand4 >> 16;
        pf->bsignificand[2] = x->f.significand3 & 0xffff;
        pf->bsignificand[3] = x->f.significand3 >> 16;
        pf->bsignificand[4] = x->f.significand2 & 0xffff;
        pf->bsignificand[5] = x->f.significand2 >> 16;
        pf->bsignificand[6] = x->f.msw.significand;
        if (x->f.msw.exponent == 0) {
                pf->blength = 7;
                pf->bexponent++;
                while (pf->bsignificand[pf->blength - 1] == 0)
                        pf->blength--;
        } else {
                pf->blength = 8;
                pf->bsignificand[7] = 1;
        }
        __shorten(pf);
}

/* PUBLIC ROUTINES */

void
single_to_decimal(single *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
        single_equivalence      *kluge;
        _big_float              bf;
        fp_exception_field_type ef;
        double                  x;

        kluge = (single_equivalence *)px;
        pd->sign = kluge->f.msw.sign;

        /* decide what to do based on the class of x */
        if (kluge->f.msw.exponent == 0) {       /* 0 or subnormal */
                if (kluge->f.msw.significand == 0) {
                        pd->fpclass = fp_zero;
                        *ps = 0;
                        return;
                } else {
#if defined(__sparc) || defined(__amd64)
                        int     i;

                        pd->fpclass = fp_subnormal;
                        /*
                         * On SPARC when nonstandard mode is enabled,
                         * or on x64 when FTZ mode is enabled, simply
                         * converting *px to double can flush a sub-
                         * normal value to zero, so we have to go
                         * through all this nonsense instead.
                         */
                        i = *(int *)px;
                        x = (double)(i & ~0x80000000);
                        if (i < 0)
                                x = -x;
                        x *= 1.401298464324817070923730e-45; /* 2^-149 */
                        ef = 0;
                        if (__fast_double_to_decimal(&x, pm, pd, &ef)) {
                                __double_to_bigfloat(&x, &bf);
                                __bigfloat_to_decimal(&bf, pm, pd, &ef);
                        }
                        if (ef != 0)
                                __base_conversion_set_exception(ef);
                        *ps = ef;
                        return;
#else
                        pd->fpclass = fp_subnormal;
#endif
                }
        } else if (kluge->f.msw.exponent == 0xff) {     /* inf or nan */
                if (kluge->f.msw.significand == 0)
                        pd->fpclass = fp_infinity;
                else if (kluge->f.msw.significand >= 0x400000)
                        pd->fpclass = fp_quiet;
                else
                        pd->fpclass = fp_signaling;
                *ps = 0;
                return;
        } else {
                pd->fpclass = fp_normal;
        }

        ef = 0;
        x = *px;
        if (__fast_double_to_decimal(&x, pm, pd, &ef)) {
                __double_to_bigfloat(&x, &bf);
                __bigfloat_to_decimal(&bf, pm, pd, &ef);
        }
        if (ef != 0)
                __base_conversion_set_exception(ef);
        *ps = ef;
}

void
double_to_decimal(double *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
        double_equivalence      *kluge;
        _big_float              bf;
        fp_exception_field_type ef;

        kluge = (double_equivalence *)px;
        pd->sign = kluge->f.msw.sign;

        /* decide what to do based on the class of x */
        if (kluge->f.msw.exponent == 0) {       /* 0 or subnormal */
                if (kluge->f.msw.significand == 0 &&
                    kluge->f.significand2 == 0) {
                        pd->fpclass = fp_zero;
                        *ps = 0;
                        return;
                } else {
                        pd->fpclass = fp_subnormal;
                }
        } else if (kluge->f.msw.exponent == 0x7ff) {    /* inf or nan */
                if (kluge->f.msw.significand == 0 &&
                    kluge->f.significand2 == 0)
                        pd->fpclass = fp_infinity;
                else if (kluge->f.msw.significand >= 0x80000)
                        pd->fpclass = fp_quiet;
                else
                        pd->fpclass = fp_signaling;
                *ps = 0;
                return;
        } else {
                pd->fpclass = fp_normal;
        }

        ef = 0;
        if (__fast_double_to_decimal(px, pm, pd, &ef)) {
                __double_to_bigfloat(px, &bf);
                __bigfloat_to_decimal(&bf, pm, pd, &ef);
        }
        if (ef != 0)
                __base_conversion_set_exception(ef);
        *ps = ef;
}

void
extended_to_decimal(extended *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
        extended_equivalence    *kluge;
        _big_float              bf;
        fp_exception_field_type ef;

        kluge = (extended_equivalence *)px;
        pd->sign = kluge->f.msw.sign;

        /* decide what to do based on the class of x */
        if (kluge->f.msw.exponent == 0) {       /* 0 or subnormal */
                if ((kluge->f.significand | kluge->f.significand2) == 0) {
                        pd->fpclass = fp_zero;
                        *ps = 0;
                        return;
                } else {
                        /*
                         * x could be a pseudo-denormal, but the distinction
                         * doesn't matter
                         */
                        pd->fpclass = fp_subnormal;
                }
        } else if ((kluge->f.significand & 0x80000000) == 0) {
                /*
                 * In Intel's extended format, if the exponent is
                 * nonzero but the explicit integer bit is zero, this
                 * is an "unsupported format" bit pattern; treat it
                 * like a signaling NaN.
                 */
                pd->fpclass = fp_signaling;
                *ps = 0;
                return;
        } else if (kluge->f.msw.exponent == 0x7fff) {   /* inf or nan */
                if (((kluge->f.significand & 0x7fffffff) |
                    kluge->f.significand2) == 0)
                        pd->fpclass = fp_infinity;
                else if ((kluge->f.significand & 0x7fffffff) >= 0x40000000)
                        pd->fpclass = fp_quiet;
                else
                        pd->fpclass = fp_signaling;
                *ps = 0;
                return;
        } else {
                pd->fpclass = fp_normal;
        }

        ef = 0;
        __extended_to_bigfloat(px, &bf);
        __bigfloat_to_decimal(&bf, pm, pd, &ef);
        if (ef != 0)
                __base_conversion_set_exception(ef);
        *ps = ef;
}

void
quadruple_to_decimal(quadruple *px, decimal_mode *pm, decimal_record *pd,
    fp_exception_field_type *ps)
{
        quadruple_equivalence   *kluge;
        _big_float              bf;
        fp_exception_field_type ef;

        kluge = (quadruple_equivalence *)px;
        pd->sign = kluge->f.msw.sign;

        /* decide what to do based on the class of x */
        if (kluge->f.msw.exponent == 0) {       /* 0 or subnormal */
                if (kluge->f.msw.significand == 0 &&
                    (kluge->f.significand2 | kluge->f.significand3 |
                    kluge->f.significand4) == 0) {
                        pd->fpclass = fp_zero;
                        *ps = 0;
                        return;
                } else {
                        pd->fpclass = fp_subnormal;
                }
        } else if (kluge->f.msw.exponent == 0x7fff) {   /* inf or nan */
                if (kluge->f.msw.significand == 0 &&
                    (kluge->f.significand2 | kluge->f.significand3 |
                    kluge->f.significand4) == 0)
                        pd->fpclass = fp_infinity;
                else if (kluge->f.msw.significand >= 0x8000)
                        pd->fpclass = fp_quiet;
                else
                        pd->fpclass = fp_signaling;
                *ps = 0;
                return;
        } else {
                pd->fpclass = fp_normal;
        }

        ef = 0;
        __quadruple_to_bigfloat(px, &bf);
        __bigfloat_to_decimal(&bf, pm, pd, &ef);
        if (ef != 0)
                __base_conversion_set_exception(ef);
        *ps = ef;
}