/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License, Version 1.0 only
 * (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 2003 Sun Microsystems, Inc.  All rights reserved.
 * Use is subject to license terms.
 */

#include "quad.h"

static const double C[] = {
        0.0,
        0.5,
        1.0,
        2.0,
        68719476736.0,
        1048576.0,
        16.0,
        1.52587890625000000000e-05,
        5.96046447753906250000e-08,
        3.72529029846191406250e-09,
        8.67361737988403547206e-19,
        2.16840434497100886801e-19,
        1.32348898008484427979e-23,
        9.62964972193617926528e-35,
        4.70197740328915003187e-38
};

#define zero    C[0]
#define half    C[1]
#define one     C[2]
#define two     C[3]
#define two36   C[4]
#define two20   C[5]
#define two4    C[6]
#define twom16  C[7]
#define twom24  C[8]
#define twom28  C[9]
#define twom60  C[10]
#define twom62  C[11]
#define twom76  C[12]
#define twom113 C[13]
#define twom124 C[14]

static const unsigned fsr_rn = 0xc0000000u;

#ifdef __sparcv9

/*
 * _Qp_mul(pz, x, y) sets *pz = *x * *y.
 */
void
_Qp_mul(union longdouble *pz, const union longdouble *x,
        const union longdouble *y)

#else

/*
 * _Q_mul(x, y) returns *x * *y.
 */
union longdouble
_Q_mul(const union longdouble *x, const union longdouble *y)

#endif  /* __sparcv9 */

{
        union longdouble        z;
        union xdouble           u;
        double                  xx[5], yy[5], c, d, zz[9];
        unsigned int            xm, ym, fsr, lx, ly, wx[3], wy[3];
        unsigned int            msw, frac2, frac3, frac4, rm;
        int                     ibit, ex, ey, ez, sign;

        xm = x->l.msw & 0x7fffffff;
        ym = y->l.msw & 0x7fffffff;
        sign = (x->l.msw ^ y->l.msw) & ~0x7fffffff;

        __quad_getfsrp(&fsr);

        /* handle nan and inf cases */
        if (xm >= 0x7fff0000 || ym >= 0x7fff0000) {
                /* handle nan cases according to V9 app. B */
                if (QUAD_ISNAN(*y)) {
                        if (!(y->l.msw & 0x8000)) {
                                /* snan, signal invalid */
                                if (fsr & FSR_NVM) {
                                        __quad_fmulq(x, y, &Z);
                                } else {
                                        Z = *y;
                                        Z.l.msw |= 0x8000;
                                        fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
                                            FSR_NVC;
                                        __quad_setfsrp(&fsr);
                                }
                                QUAD_RETURN(Z);
                        } else if (QUAD_ISNAN(*x) && !(x->l.msw & 0x8000)) {
                                /* snan, signal invalid */
                                if (fsr & FSR_NVM) {
                                        __quad_fmulq(x, y, &Z);
                                } else {
                                        Z = *x;
                                        Z.l.msw |= 0x8000;
                                        fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
                                            FSR_NVC;
                                        __quad_setfsrp(&fsr);
                                }
                                QUAD_RETURN(Z);
                        }
                        Z = *y;
                        QUAD_RETURN(Z);
                }
                if (QUAD_ISNAN(*x)) {
                        if (!(x->l.msw & 0x8000)) {
                                /* snan, signal invalid */
                                if (fsr & FSR_NVM) {
                                        __quad_fmulq(x, y, &Z);
                                } else {
                                        Z = *x;
                                        Z.l.msw |= 0x8000;
                                        fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
                                            FSR_NVC;
                                        __quad_setfsrp(&fsr);
                                }
                                QUAD_RETURN(Z);
                        }
                        Z = *x;
                        QUAD_RETURN(Z);
                }
                if (xm == 0x7fff0000) {
                        /* x is inf */
                        if (QUAD_ISZERO(*y)) {
                                /* zero * inf, signal invalid */
                                if (fsr & FSR_NVM) {
                                        __quad_fmulq(x, y, &Z);
                                } else {
                                        Z.l.msw = 0x7fffffff;
                                        Z.l.frac2 = Z.l.frac3 =
                                            Z.l.frac4 = 0xffffffff;
                                        fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
                                            FSR_NVC;
                                        __quad_setfsrp(&fsr);
                                }
                                QUAD_RETURN(Z);
                        }
                        /* inf * finite, return inf */
                        Z.l.msw = sign | 0x7fff0000;
                        Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
                        QUAD_RETURN(Z);
                }
                /* y is inf */
                if (QUAD_ISZERO(*x)) {
                        /* zero * inf, signal invalid */
                        if (fsr & FSR_NVM) {
                                __quad_fmulq(x, y, &Z);
                        } else {
                                Z.l.msw = 0x7fffffff;
                                Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff;
                                fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC;
                                __quad_setfsrp(&fsr);
                        }
                        QUAD_RETURN(Z);
                }
                /* inf * finite, return inf */
                Z.l.msw = sign | 0x7fff0000;
                Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
                QUAD_RETURN(Z);
        }

        /* handle zero cases */
        if (xm == 0 || ym == 0) {
                if (QUAD_ISZERO(*x) || QUAD_ISZERO(*y)) {
                        Z.l.msw = sign;
                        Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
                        QUAD_RETURN(Z);
                }
        }

        /* now x and y are finite, nonzero */
        __quad_setfsrp(&fsr_rn);

        /* get their normalized significands and exponents */
        ex = (int)(xm >> 16);
        lx = xm & 0xffff;
        if (ex) {
                lx |= 0x10000;
                wx[0] = x->l.frac2;
                wx[1] = x->l.frac3;
                wx[2] = x->l.frac4;
        } else {
                if (lx | (x->l.frac2 & 0xfffe0000)) {
                        wx[0] = x->l.frac2;
                        wx[1] = x->l.frac3;
                        wx[2] = x->l.frac4;
                        ex = 1;
                } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
                        lx = x->l.frac2;
                        wx[0] = x->l.frac3;
                        wx[1] = x->l.frac4;
                        wx[2] = 0;
                        ex = -31;
                } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
                        lx = x->l.frac3;
                        wx[0] = x->l.frac4;
                        wx[1] = wx[2] = 0;
                        ex = -63;
                } else {
                        lx = x->l.frac4;
                        wx[0] = wx[1] = wx[2] = 0;
                        ex = -95;
                }
                while ((lx & 0x10000) == 0) {
                        lx = (lx << 1) | (wx[0] >> 31);
                        wx[0] = (wx[0] << 1) | (wx[1] >> 31);
                        wx[1] = (wx[1] << 1) | (wx[2] >> 31);
                        wx[2] <<= 1;
                        ex--;
                }
        }
        ez = ex - 0x3fff;

        ey = (int)(ym >> 16);
        ly = ym & 0xffff;
        if (ey) {
                ly |= 0x10000;
                wy[0] = y->l.frac2;
                wy[1] = y->l.frac3;
                wy[2] = y->l.frac4;
        } else {
                if (ly | (y->l.frac2 & 0xfffe0000)) {
                        wy[0] = y->l.frac2;
                        wy[1] = y->l.frac3;
                        wy[2] = y->l.frac4;
                        ey = 1;
                } else if (y->l.frac2 | (y->l.frac3 & 0xfffe0000)) {
                        ly = y->l.frac2;
                        wy[0] = y->l.frac3;
                        wy[1] = y->l.frac4;
                        wy[2] = 0;
                        ey = -31;
                } else if (y->l.frac3 | (y->l.frac4 & 0xfffe0000)) {
                        ly = y->l.frac3;
                        wy[0] = y->l.frac4;
                        wy[1] = wy[2] = 0;
                        ey = -63;
                } else {
                        ly = y->l.frac4;
                        wy[0] = wy[1] = wy[2] = 0;
                        ey = -95;
                }
                while ((ly & 0x10000) == 0) {
                        ly = (ly << 1) | (wy[0] >> 31);
                        wy[0] = (wy[0] << 1) | (wy[1] >> 31);
                        wy[1] = (wy[1] << 1) | (wy[2] >> 31);
                        wy[2] <<= 1;
                        ey--;
                }
        }
        ez += ey;

        /* extract the significand into five doubles */
        c = twom16;
        xx[0] = (double)((int)lx) * c;
        yy[0] = (double)((int)ly) * c;

        c *= twom24;
        xx[1] = (double)((int)(wx[0] >> 8)) * c;
        yy[1] = (double)((int)(wy[0] >> 8)) * c;

        c *= twom24;
        xx[2] = (double)((int)(((wx[0] << 16) | (wx[1] >> 16)) &
            0xffffff)) * c;
        yy[2] = (double)((int)(((wy[0] << 16) | (wy[1] >> 16)) &
            0xffffff)) * c;

        c *= twom24;
        xx[3] = (double)((int)(((wx[1] << 8) | (wx[2] >> 24)) &
            0xffffff)) * c;
        yy[3] = (double)((int)(((wy[1] << 8) | (wy[2] >> 24)) &
            0xffffff)) * c;

        c *= twom24;
        xx[4] = (double)((int)(wx[2] & 0xffffff)) * c;
        yy[4] = (double)((int)(wy[2] & 0xffffff)) * c;

        /* form the "digits" of the product */
        zz[0] = xx[0] * yy[0];
        zz[1] = xx[0] * yy[1] + xx[1] * yy[0];
        zz[2] = xx[0] * yy[2] + xx[1] * yy[1] + xx[2] * yy[0];
        zz[3] = xx[0] * yy[3] + xx[1] * yy[2] + xx[2] * yy[1] +
            xx[3] * yy[0];
        zz[4] = xx[0] * yy[4] + xx[1] * yy[3] + xx[2] * yy[2] +
            xx[3] * yy[1] + xx[4] * yy[0];
        zz[5] = xx[1] * yy[4] + xx[2] * yy[3] + xx[3] * yy[2] +
            xx[4] * yy[1];
        zz[6] = xx[2] * yy[4] + xx[3] * yy[3] + xx[4] * yy[2];
        zz[7] = xx[3] * yy[4] + xx[4] * yy[3];
        zz[8] = xx[4] * yy[4];

        /* collect the first few terms */
        c = (zz[1] + two20) - two20;
        zz[0] += c;
        zz[1] = zz[2] + (zz[1] - c);
        c = (zz[3] + twom28) - twom28;
        zz[1] += c;
        zz[2] = zz[4] + (zz[3] - c);

        /* propagate carries into the third term */
        d = zz[6] + (zz[7] + zz[8]);
        zz[2] += zz[5] + d;

        /* if the third term might lie on a rounding boundary, perturb it */
        /* as need be */
        if (zz[2] == (twom62 + zz[2]) - twom62)
        {
                c = (zz[5] + twom76) - twom76;
                if ((zz[5] - c) + d != zero)
                        zz[2] += twom124;
        }

        /* propagate carries to the leading term */
        c = zz[1] + zz[2];
        zz[2] = zz[2] + (zz[1] - c);
        zz[1] = c;
        c = zz[0] + zz[1];
        zz[1] = zz[1] + (zz[0] - c);
        zz[0] = c;

        /* check for carry out */
        if (c >= two) {
                /* postnormalize */
                zz[0] *= half;
                zz[1] *= half;
                zz[2] *= half;
                ez++;
        }

        /* if exponent > 0 strip off integer bit, else denormalize */
        if (ez > 0) {
                ibit = 1;
                zz[0] -= one;
        } else {
                ibit = 0;
                if (ez > -128)
                        u.l.hi = (unsigned)(0x3fe + ez) << 20;
                else
                        u.l.hi = 0x37e00000;
                u.l.lo = 0;
                zz[0] *= u.d;
                zz[1] *= u.d;
                zz[2] *= u.d;
                ez = 0;
        }

        /* the first 48 bits of fraction come from zz[0] */
        u.d = d = two36 + zz[0];
        msw = u.l.lo;
        zz[0] -= (d - two36);

        u.d = d = two4 + zz[0];
        frac2 = u.l.lo;
        zz[0] -= (d - two4);

        /* the next 32 come from zz[0] and zz[1] */
        u.d = d = twom28 + (zz[0] + zz[1]);
        frac3 = u.l.lo;
        zz[0] -= (d - twom28);

        /* condense the remaining fraction; errors here won't matter */
        c = zz[0] + zz[1];
        zz[1] = ((zz[0] - c) + zz[1]) + zz[2];
        zz[0] = c;

        /* get the last word of fraction */
        u.d = d = twom60 + (zz[0] + zz[1]);
        frac4 = u.l.lo;
        zz[0] -= (d - twom60);

        /* keep track of what's left for rounding; note that the error */
        /* in computing c will be non-negative due to rounding mode */
        c = zz[0] + zz[1];

        /* get the rounding mode, fudging directed rounding modes */
        /* as though the result were positive */
        rm = fsr >> 30;
        if (sign)
                rm ^= (rm >> 1);

        /* round and raise exceptions */
        fsr &= ~FSR_CEXC;
        if (c != zero) {
                fsr |= FSR_NXC;

                /* decide whether to round the fraction up */
                if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 ||
                    (c == twom113 && ((frac4 & 1) || (c - zz[0] != zz[1])))))) {
                        /* round up and renormalize if necessary */
                        if (++frac4 == 0)
                                if (++frac3 == 0)
                                        if (++frac2 == 0)
                                                if (++msw == 0x10000) {
                                                        msw = 0;
                                                        ez++;
                                                }
                }
        }

        /* check for under/overflow */
        if (ez >= 0x7fff) {
                if (rm == FSR_RN || rm == FSR_RP) {
                        z.l.msw = sign | 0x7fff0000;
                        z.l.frac2 = z.l.frac3 = z.l.frac4 = 0;
                } else {
                        z.l.msw = sign | 0x7ffeffff;
                        z.l.frac2 = z.l.frac3 = z.l.frac4 = 0xffffffff;
                }
                fsr |= FSR_OFC | FSR_NXC;
        } else {
                z.l.msw = sign | (ez << 16) | msw;
                z.l.frac2 = frac2;
                z.l.frac3 = frac3;
                z.l.frac4 = frac4;

                /* !ibit => exact result was tiny before rounding, */
                /* t nonzero => result delivered is inexact */
                if (!ibit) {
                        if (c != zero)
                                fsr |= FSR_UFC | FSR_NXC;
                        else if (fsr & FSR_UFM)
                                fsr |= FSR_UFC;
                }
        }

        if ((fsr & FSR_CEXC) & (fsr >> 23)) {
                __quad_setfsrp(&fsr);
                __quad_fmulq(x, y, &Z);
        } else {
                Z = z;
                fsr |= (fsr & 0x1f) << 5;
                __quad_setfsrp(&fsr);
        }
        QUAD_RETURN(Z);
}