/*
* 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
* 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.
*/
#pragma ident "%Z%%M% %I% %E% SMI"
#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]
#ifdef __sparcv9
/*
* _Qp_mul(pz, x, y) sets *pz = *x * *y.
*/
void
const union longdouble *y)
#else
/*
* _Q_mul(x, y) returns *x * *y.
*/
union longdouble
#endif /* __sparcv9 */
{
union longdouble z;
union xdouble u;
/* handle nan and inf cases */
/* handle nan cases according to V9 app. B */
if (QUAD_ISNAN(*y)) {
if (!(y->l.msw & 0x8000)) {
/* snan, signal invalid */
__quad_fmulq(x, y, &Z);
} else {
Z = *y;
Z.l.msw |= 0x8000;
}
QUAD_RETURN(Z);
/* snan, signal invalid */
__quad_fmulq(x, y, &Z);
} else {
Z = *x;
Z.l.msw |= 0x8000;
}
QUAD_RETURN(Z);
}
Z = *y;
QUAD_RETURN(Z);
}
if (QUAD_ISNAN(*x)) {
if (!(x->l.msw & 0x8000)) {
/* snan, signal invalid */
__quad_fmulq(x, y, &Z);
} else {
Z = *x;
Z.l.msw |= 0x8000;
}
QUAD_RETURN(Z);
}
Z = *x;
QUAD_RETURN(Z);
}
if (xm == 0x7fff0000) {
/* x is inf */
if (QUAD_ISZERO(*y)) {
/* zero * inf, signal invalid */
__quad_fmulq(x, y, &Z);
} else {
Z.l.msw = 0x7fffffff;
Z.l.frac4 = 0xffffffff;
}
QUAD_RETURN(Z);
}
/* inf * finite, return inf */
QUAD_RETURN(Z);
}
/* y is inf */
if (QUAD_ISZERO(*x)) {
/* zero * inf, signal invalid */
__quad_fmulq(x, y, &Z);
} else {
Z.l.msw = 0x7fffffff;
}
QUAD_RETURN(Z);
}
/* inf * finite, return inf */
QUAD_RETURN(Z);
}
/* handle zero cases */
if (QUAD_ISZERO(*x) || QUAD_ISZERO(*y)) {
QUAD_RETURN(Z);
}
}
/* now x and y are finite, nonzero */
/* get their normalized significands and exponents */
if (ex) {
lx |= 0x10000;
} else {
ex = 1;
wx[2] = 0;
ex = -31;
ex = -63;
} else {
ex = -95;
}
while ((lx & 0x10000) == 0) {
ex--;
}
}
if (ey) {
ly |= 0x10000;
} else {
ey = 1;
wy[2] = 0;
ey = -31;
ey = -63;
} else {
ey = -95;
}
while ((ly & 0x10000) == 0) {
ey--;
}
}
/* extract the significand into five doubles */
c = twom16;
c *= twom24;
c *= twom24;
0xffffff)) * c;
0xffffff)) * c;
c *= twom24;
0xffffff)) * c;
0xffffff)) * c;
c *= twom24;
/* form the "digits" of the product */
/* collect the first few terms */
zz[0] += c;
zz[1] += c;
/* propagate carries into the third term */
/* if the third term might lie on a rounding boundary, perturb it */
/* as need be */
{
}
/* propagate carries to the leading term */
zz[1] = c;
zz[0] = c;
/* check for carry out */
if (c >= two) {
/* postnormalize */
ez++;
}
/* if exponent > 0 strip off integer bit, else denormalize */
if (ez > 0) {
ibit = 1;
} else {
ibit = 0;
if (ez > -128)
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] */
/* the next 32 come from zz[0] and zz[1] */
/* condense the remaining fraction; errors here won't matter */
zz[0] = c;
/* get the last word of fraction */
/* keep track of what's left for rounding; note that the error */
/* in computing c will be non-negative due to rounding mode */
/* get the rounding mode, fudging directed rounding modes */
/* as though the result were positive */
if (sign)
/* round and raise exceptions */
if (c != zero) {
/* decide whether to round the fraction up */
/* round up and renormalize if necessary */
if (++frac4 == 0)
if (++frac3 == 0)
if (++frac2 == 0)
if (++msw == 0x10000) {
msw = 0;
ez++;
}
}
}
if (ez >= 0x7fff) {
} else {
}
} else {
/* !ibit => exact result was tiny before rounding, */
/* t nonzero => result delivered is inexact */
if (!ibit) {
if (c != zero)
}
}
__quad_fmulq(x, y, &Z);
} else {
Z = z;
}
QUAD_RETURN(Z);
}