1*0Sstevel@tonic-gate /*
2*0Sstevel@tonic-gate * CDDL HEADER START
3*0Sstevel@tonic-gate *
4*0Sstevel@tonic-gate * The contents of this file are subject to the terms of the
5*0Sstevel@tonic-gate * Common Development and Distribution License, Version 1.0 only
6*0Sstevel@tonic-gate * (the "License"). You may not use this file except in compliance
7*0Sstevel@tonic-gate * with the License.
8*0Sstevel@tonic-gate *
9*0Sstevel@tonic-gate * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10*0Sstevel@tonic-gate * or http://www.opensolaris.org/os/licensing.
11*0Sstevel@tonic-gate * See the License for the specific language governing permissions
12*0Sstevel@tonic-gate * and limitations under the License.
13*0Sstevel@tonic-gate *
14*0Sstevel@tonic-gate * When distributing Covered Code, include this CDDL HEADER in each
15*0Sstevel@tonic-gate * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16*0Sstevel@tonic-gate * If applicable, add the following below this CDDL HEADER, with the
17*0Sstevel@tonic-gate * fields enclosed by brackets "[]" replaced with your own identifying
18*0Sstevel@tonic-gate * information: Portions Copyright [yyyy] [name of copyright owner]
19*0Sstevel@tonic-gate *
20*0Sstevel@tonic-gate * CDDL HEADER END
21*0Sstevel@tonic-gate */
22*0Sstevel@tonic-gate /*
23*0Sstevel@tonic-gate * Copyright 2003 Sun Microsystems, Inc. All rights reserved.
24*0Sstevel@tonic-gate * Use is subject to license terms.
25*0Sstevel@tonic-gate */
26*0Sstevel@tonic-gate
27*0Sstevel@tonic-gate #pragma ident "%Z%%M% %I% %E% SMI"
28*0Sstevel@tonic-gate
29*0Sstevel@tonic-gate #include "quad.h"
30*0Sstevel@tonic-gate
31*0Sstevel@tonic-gate static const double C[] = {
32*0Sstevel@tonic-gate 0.0,
33*0Sstevel@tonic-gate 0.5,
34*0Sstevel@tonic-gate 1.0,
35*0Sstevel@tonic-gate 2.0,
36*0Sstevel@tonic-gate 68719476736.0,
37*0Sstevel@tonic-gate 1048576.0,
38*0Sstevel@tonic-gate 16.0,
39*0Sstevel@tonic-gate 1.52587890625000000000e-05,
40*0Sstevel@tonic-gate 5.96046447753906250000e-08,
41*0Sstevel@tonic-gate 3.72529029846191406250e-09,
42*0Sstevel@tonic-gate 8.67361737988403547206e-19,
43*0Sstevel@tonic-gate 2.16840434497100886801e-19,
44*0Sstevel@tonic-gate 1.32348898008484427979e-23,
45*0Sstevel@tonic-gate 9.62964972193617926528e-35,
46*0Sstevel@tonic-gate 4.70197740328915003187e-38
47*0Sstevel@tonic-gate };
48*0Sstevel@tonic-gate
49*0Sstevel@tonic-gate #define zero C[0]
50*0Sstevel@tonic-gate #define half C[1]
51*0Sstevel@tonic-gate #define one C[2]
52*0Sstevel@tonic-gate #define two C[3]
53*0Sstevel@tonic-gate #define two36 C[4]
54*0Sstevel@tonic-gate #define two20 C[5]
55*0Sstevel@tonic-gate #define two4 C[6]
56*0Sstevel@tonic-gate #define twom16 C[7]
57*0Sstevel@tonic-gate #define twom24 C[8]
58*0Sstevel@tonic-gate #define twom28 C[9]
59*0Sstevel@tonic-gate #define twom60 C[10]
60*0Sstevel@tonic-gate #define twom62 C[11]
61*0Sstevel@tonic-gate #define twom76 C[12]
62*0Sstevel@tonic-gate #define twom113 C[13]
63*0Sstevel@tonic-gate #define twom124 C[14]
64*0Sstevel@tonic-gate
65*0Sstevel@tonic-gate static const unsigned fsr_rn = 0xc0000000u;
66*0Sstevel@tonic-gate
67*0Sstevel@tonic-gate #ifdef __sparcv9
68*0Sstevel@tonic-gate
69*0Sstevel@tonic-gate /*
70*0Sstevel@tonic-gate * _Qp_mul(pz, x, y) sets *pz = *x * *y.
71*0Sstevel@tonic-gate */
72*0Sstevel@tonic-gate void
_Qp_mul(union longdouble * pz,const union longdouble * x,const union longdouble * y)73*0Sstevel@tonic-gate _Qp_mul(union longdouble *pz, const union longdouble *x,
74*0Sstevel@tonic-gate const union longdouble *y)
75*0Sstevel@tonic-gate
76*0Sstevel@tonic-gate #else
77*0Sstevel@tonic-gate
78*0Sstevel@tonic-gate /*
79*0Sstevel@tonic-gate * _Q_mul(x, y) returns *x * *y.
80*0Sstevel@tonic-gate */
81*0Sstevel@tonic-gate union longdouble
82*0Sstevel@tonic-gate _Q_mul(const union longdouble *x, const union longdouble *y)
83*0Sstevel@tonic-gate
84*0Sstevel@tonic-gate #endif /* __sparcv9 */
85*0Sstevel@tonic-gate
86*0Sstevel@tonic-gate {
87*0Sstevel@tonic-gate union longdouble z;
88*0Sstevel@tonic-gate union xdouble u;
89*0Sstevel@tonic-gate double xx[5], yy[5], c, d, zz[9];
90*0Sstevel@tonic-gate unsigned int xm, ym, fsr, lx, ly, wx[3], wy[3];
91*0Sstevel@tonic-gate unsigned int msw, frac2, frac3, frac4, rm;
92*0Sstevel@tonic-gate int ibit, ex, ey, ez, sign;
93*0Sstevel@tonic-gate
94*0Sstevel@tonic-gate xm = x->l.msw & 0x7fffffff;
95*0Sstevel@tonic-gate ym = y->l.msw & 0x7fffffff;
96*0Sstevel@tonic-gate sign = (x->l.msw ^ y->l.msw) & ~0x7fffffff;
97*0Sstevel@tonic-gate
98*0Sstevel@tonic-gate __quad_getfsrp(&fsr);
99*0Sstevel@tonic-gate
100*0Sstevel@tonic-gate /* handle nan and inf cases */
101*0Sstevel@tonic-gate if (xm >= 0x7fff0000 || ym >= 0x7fff0000) {
102*0Sstevel@tonic-gate /* handle nan cases according to V9 app. B */
103*0Sstevel@tonic-gate if (QUAD_ISNAN(*y)) {
104*0Sstevel@tonic-gate if (!(y->l.msw & 0x8000)) {
105*0Sstevel@tonic-gate /* snan, signal invalid */
106*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
107*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
108*0Sstevel@tonic-gate } else {
109*0Sstevel@tonic-gate Z = *y;
110*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
111*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
112*0Sstevel@tonic-gate FSR_NVC;
113*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
114*0Sstevel@tonic-gate }
115*0Sstevel@tonic-gate QUAD_RETURN(Z);
116*0Sstevel@tonic-gate } else if (QUAD_ISNAN(*x) && !(x->l.msw & 0x8000)) {
117*0Sstevel@tonic-gate /* snan, signal invalid */
118*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
119*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
120*0Sstevel@tonic-gate } else {
121*0Sstevel@tonic-gate Z = *x;
122*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
123*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
124*0Sstevel@tonic-gate FSR_NVC;
125*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
126*0Sstevel@tonic-gate }
127*0Sstevel@tonic-gate QUAD_RETURN(Z);
128*0Sstevel@tonic-gate }
129*0Sstevel@tonic-gate Z = *y;
130*0Sstevel@tonic-gate QUAD_RETURN(Z);
131*0Sstevel@tonic-gate }
132*0Sstevel@tonic-gate if (QUAD_ISNAN(*x)) {
133*0Sstevel@tonic-gate if (!(x->l.msw & 0x8000)) {
134*0Sstevel@tonic-gate /* snan, signal invalid */
135*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
136*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
137*0Sstevel@tonic-gate } else {
138*0Sstevel@tonic-gate Z = *x;
139*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
140*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
141*0Sstevel@tonic-gate FSR_NVC;
142*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
143*0Sstevel@tonic-gate }
144*0Sstevel@tonic-gate QUAD_RETURN(Z);
145*0Sstevel@tonic-gate }
146*0Sstevel@tonic-gate Z = *x;
147*0Sstevel@tonic-gate QUAD_RETURN(Z);
148*0Sstevel@tonic-gate }
149*0Sstevel@tonic-gate if (xm == 0x7fff0000) {
150*0Sstevel@tonic-gate /* x is inf */
151*0Sstevel@tonic-gate if (QUAD_ISZERO(*y)) {
152*0Sstevel@tonic-gate /* zero * inf, signal invalid */
153*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
154*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
155*0Sstevel@tonic-gate } else {
156*0Sstevel@tonic-gate Z.l.msw = 0x7fffffff;
157*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 =
158*0Sstevel@tonic-gate Z.l.frac4 = 0xffffffff;
159*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
160*0Sstevel@tonic-gate FSR_NVC;
161*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
162*0Sstevel@tonic-gate }
163*0Sstevel@tonic-gate QUAD_RETURN(Z);
164*0Sstevel@tonic-gate }
165*0Sstevel@tonic-gate /* inf * finite, return inf */
166*0Sstevel@tonic-gate Z.l.msw = sign | 0x7fff0000;
167*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
168*0Sstevel@tonic-gate QUAD_RETURN(Z);
169*0Sstevel@tonic-gate }
170*0Sstevel@tonic-gate /* y is inf */
171*0Sstevel@tonic-gate if (QUAD_ISZERO(*x)) {
172*0Sstevel@tonic-gate /* zero * inf, signal invalid */
173*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
174*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
175*0Sstevel@tonic-gate } else {
176*0Sstevel@tonic-gate Z.l.msw = 0x7fffffff;
177*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff;
178*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC;
179*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
180*0Sstevel@tonic-gate }
181*0Sstevel@tonic-gate QUAD_RETURN(Z);
182*0Sstevel@tonic-gate }
183*0Sstevel@tonic-gate /* inf * finite, return inf */
184*0Sstevel@tonic-gate Z.l.msw = sign | 0x7fff0000;
185*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
186*0Sstevel@tonic-gate QUAD_RETURN(Z);
187*0Sstevel@tonic-gate }
188*0Sstevel@tonic-gate
189*0Sstevel@tonic-gate /* handle zero cases */
190*0Sstevel@tonic-gate if (xm == 0 || ym == 0) {
191*0Sstevel@tonic-gate if (QUAD_ISZERO(*x) || QUAD_ISZERO(*y)) {
192*0Sstevel@tonic-gate Z.l.msw = sign;
193*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
194*0Sstevel@tonic-gate QUAD_RETURN(Z);
195*0Sstevel@tonic-gate }
196*0Sstevel@tonic-gate }
197*0Sstevel@tonic-gate
198*0Sstevel@tonic-gate /* now x and y are finite, nonzero */
199*0Sstevel@tonic-gate __quad_setfsrp(&fsr_rn);
200*0Sstevel@tonic-gate
201*0Sstevel@tonic-gate /* get their normalized significands and exponents */
202*0Sstevel@tonic-gate ex = (int)(xm >> 16);
203*0Sstevel@tonic-gate lx = xm & 0xffff;
204*0Sstevel@tonic-gate if (ex) {
205*0Sstevel@tonic-gate lx |= 0x10000;
206*0Sstevel@tonic-gate wx[0] = x->l.frac2;
207*0Sstevel@tonic-gate wx[1] = x->l.frac3;
208*0Sstevel@tonic-gate wx[2] = x->l.frac4;
209*0Sstevel@tonic-gate } else {
210*0Sstevel@tonic-gate if (lx | (x->l.frac2 & 0xfffe0000)) {
211*0Sstevel@tonic-gate wx[0] = x->l.frac2;
212*0Sstevel@tonic-gate wx[1] = x->l.frac3;
213*0Sstevel@tonic-gate wx[2] = x->l.frac4;
214*0Sstevel@tonic-gate ex = 1;
215*0Sstevel@tonic-gate } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
216*0Sstevel@tonic-gate lx = x->l.frac2;
217*0Sstevel@tonic-gate wx[0] = x->l.frac3;
218*0Sstevel@tonic-gate wx[1] = x->l.frac4;
219*0Sstevel@tonic-gate wx[2] = 0;
220*0Sstevel@tonic-gate ex = -31;
221*0Sstevel@tonic-gate } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
222*0Sstevel@tonic-gate lx = x->l.frac3;
223*0Sstevel@tonic-gate wx[0] = x->l.frac4;
224*0Sstevel@tonic-gate wx[1] = wx[2] = 0;
225*0Sstevel@tonic-gate ex = -63;
226*0Sstevel@tonic-gate } else {
227*0Sstevel@tonic-gate lx = x->l.frac4;
228*0Sstevel@tonic-gate wx[0] = wx[1] = wx[2] = 0;
229*0Sstevel@tonic-gate ex = -95;
230*0Sstevel@tonic-gate }
231*0Sstevel@tonic-gate while ((lx & 0x10000) == 0) {
232*0Sstevel@tonic-gate lx = (lx << 1) | (wx[0] >> 31);
233*0Sstevel@tonic-gate wx[0] = (wx[0] << 1) | (wx[1] >> 31);
234*0Sstevel@tonic-gate wx[1] = (wx[1] << 1) | (wx[2] >> 31);
235*0Sstevel@tonic-gate wx[2] <<= 1;
236*0Sstevel@tonic-gate ex--;
237*0Sstevel@tonic-gate }
238*0Sstevel@tonic-gate }
239*0Sstevel@tonic-gate ez = ex - 0x3fff;
240*0Sstevel@tonic-gate
241*0Sstevel@tonic-gate ey = (int)(ym >> 16);
242*0Sstevel@tonic-gate ly = ym & 0xffff;
243*0Sstevel@tonic-gate if (ey) {
244*0Sstevel@tonic-gate ly |= 0x10000;
245*0Sstevel@tonic-gate wy[0] = y->l.frac2;
246*0Sstevel@tonic-gate wy[1] = y->l.frac3;
247*0Sstevel@tonic-gate wy[2] = y->l.frac4;
248*0Sstevel@tonic-gate } else {
249*0Sstevel@tonic-gate if (ly | (y->l.frac2 & 0xfffe0000)) {
250*0Sstevel@tonic-gate wy[0] = y->l.frac2;
251*0Sstevel@tonic-gate wy[1] = y->l.frac3;
252*0Sstevel@tonic-gate wy[2] = y->l.frac4;
253*0Sstevel@tonic-gate ey = 1;
254*0Sstevel@tonic-gate } else if (y->l.frac2 | (y->l.frac3 & 0xfffe0000)) {
255*0Sstevel@tonic-gate ly = y->l.frac2;
256*0Sstevel@tonic-gate wy[0] = y->l.frac3;
257*0Sstevel@tonic-gate wy[1] = y->l.frac4;
258*0Sstevel@tonic-gate wy[2] = 0;
259*0Sstevel@tonic-gate ey = -31;
260*0Sstevel@tonic-gate } else if (y->l.frac3 | (y->l.frac4 & 0xfffe0000)) {
261*0Sstevel@tonic-gate ly = y->l.frac3;
262*0Sstevel@tonic-gate wy[0] = y->l.frac4;
263*0Sstevel@tonic-gate wy[1] = wy[2] = 0;
264*0Sstevel@tonic-gate ey = -63;
265*0Sstevel@tonic-gate } else {
266*0Sstevel@tonic-gate ly = y->l.frac4;
267*0Sstevel@tonic-gate wy[0] = wy[1] = wy[2] = 0;
268*0Sstevel@tonic-gate ey = -95;
269*0Sstevel@tonic-gate }
270*0Sstevel@tonic-gate while ((ly & 0x10000) == 0) {
271*0Sstevel@tonic-gate ly = (ly << 1) | (wy[0] >> 31);
272*0Sstevel@tonic-gate wy[0] = (wy[0] << 1) | (wy[1] >> 31);
273*0Sstevel@tonic-gate wy[1] = (wy[1] << 1) | (wy[2] >> 31);
274*0Sstevel@tonic-gate wy[2] <<= 1;
275*0Sstevel@tonic-gate ey--;
276*0Sstevel@tonic-gate }
277*0Sstevel@tonic-gate }
278*0Sstevel@tonic-gate ez += ey;
279*0Sstevel@tonic-gate
280*0Sstevel@tonic-gate /* extract the significand into five doubles */
281*0Sstevel@tonic-gate c = twom16;
282*0Sstevel@tonic-gate xx[0] = (double)((int)lx) * c;
283*0Sstevel@tonic-gate yy[0] = (double)((int)ly) * c;
284*0Sstevel@tonic-gate
285*0Sstevel@tonic-gate c *= twom24;
286*0Sstevel@tonic-gate xx[1] = (double)((int)(wx[0] >> 8)) * c;
287*0Sstevel@tonic-gate yy[1] = (double)((int)(wy[0] >> 8)) * c;
288*0Sstevel@tonic-gate
289*0Sstevel@tonic-gate c *= twom24;
290*0Sstevel@tonic-gate xx[2] = (double)((int)(((wx[0] << 16) | (wx[1] >> 16)) &
291*0Sstevel@tonic-gate 0xffffff)) * c;
292*0Sstevel@tonic-gate yy[2] = (double)((int)(((wy[0] << 16) | (wy[1] >> 16)) &
293*0Sstevel@tonic-gate 0xffffff)) * c;
294*0Sstevel@tonic-gate
295*0Sstevel@tonic-gate c *= twom24;
296*0Sstevel@tonic-gate xx[3] = (double)((int)(((wx[1] << 8) | (wx[2] >> 24)) &
297*0Sstevel@tonic-gate 0xffffff)) * c;
298*0Sstevel@tonic-gate yy[3] = (double)((int)(((wy[1] << 8) | (wy[2] >> 24)) &
299*0Sstevel@tonic-gate 0xffffff)) * c;
300*0Sstevel@tonic-gate
301*0Sstevel@tonic-gate c *= twom24;
302*0Sstevel@tonic-gate xx[4] = (double)((int)(wx[2] & 0xffffff)) * c;
303*0Sstevel@tonic-gate yy[4] = (double)((int)(wy[2] & 0xffffff)) * c;
304*0Sstevel@tonic-gate
305*0Sstevel@tonic-gate /* form the "digits" of the product */
306*0Sstevel@tonic-gate zz[0] = xx[0] * yy[0];
307*0Sstevel@tonic-gate zz[1] = xx[0] * yy[1] + xx[1] * yy[0];
308*0Sstevel@tonic-gate zz[2] = xx[0] * yy[2] + xx[1] * yy[1] + xx[2] * yy[0];
309*0Sstevel@tonic-gate zz[3] = xx[0] * yy[3] + xx[1] * yy[2] + xx[2] * yy[1] +
310*0Sstevel@tonic-gate xx[3] * yy[0];
311*0Sstevel@tonic-gate zz[4] = xx[0] * yy[4] + xx[1] * yy[3] + xx[2] * yy[2] +
312*0Sstevel@tonic-gate xx[3] * yy[1] + xx[4] * yy[0];
313*0Sstevel@tonic-gate zz[5] = xx[1] * yy[4] + xx[2] * yy[3] + xx[3] * yy[2] +
314*0Sstevel@tonic-gate xx[4] * yy[1];
315*0Sstevel@tonic-gate zz[6] = xx[2] * yy[4] + xx[3] * yy[3] + xx[4] * yy[2];
316*0Sstevel@tonic-gate zz[7] = xx[3] * yy[4] + xx[4] * yy[3];
317*0Sstevel@tonic-gate zz[8] = xx[4] * yy[4];
318*0Sstevel@tonic-gate
319*0Sstevel@tonic-gate /* collect the first few terms */
320*0Sstevel@tonic-gate c = (zz[1] + two20) - two20;
321*0Sstevel@tonic-gate zz[0] += c;
322*0Sstevel@tonic-gate zz[1] = zz[2] + (zz[1] - c);
323*0Sstevel@tonic-gate c = (zz[3] + twom28) - twom28;
324*0Sstevel@tonic-gate zz[1] += c;
325*0Sstevel@tonic-gate zz[2] = zz[4] + (zz[3] - c);
326*0Sstevel@tonic-gate
327*0Sstevel@tonic-gate /* propagate carries into the third term */
328*0Sstevel@tonic-gate d = zz[6] + (zz[7] + zz[8]);
329*0Sstevel@tonic-gate zz[2] += zz[5] + d;
330*0Sstevel@tonic-gate
331*0Sstevel@tonic-gate /* if the third term might lie on a rounding boundary, perturb it */
332*0Sstevel@tonic-gate /* as need be */
333*0Sstevel@tonic-gate if (zz[2] == (twom62 + zz[2]) - twom62)
334*0Sstevel@tonic-gate {
335*0Sstevel@tonic-gate c = (zz[5] + twom76) - twom76;
336*0Sstevel@tonic-gate if ((zz[5] - c) + d != zero)
337*0Sstevel@tonic-gate zz[2] += twom124;
338*0Sstevel@tonic-gate }
339*0Sstevel@tonic-gate
340*0Sstevel@tonic-gate /* propagate carries to the leading term */
341*0Sstevel@tonic-gate c = zz[1] + zz[2];
342*0Sstevel@tonic-gate zz[2] = zz[2] + (zz[1] - c);
343*0Sstevel@tonic-gate zz[1] = c;
344*0Sstevel@tonic-gate c = zz[0] + zz[1];
345*0Sstevel@tonic-gate zz[1] = zz[1] + (zz[0] - c);
346*0Sstevel@tonic-gate zz[0] = c;
347*0Sstevel@tonic-gate
348*0Sstevel@tonic-gate /* check for carry out */
349*0Sstevel@tonic-gate if (c >= two) {
350*0Sstevel@tonic-gate /* postnormalize */
351*0Sstevel@tonic-gate zz[0] *= half;
352*0Sstevel@tonic-gate zz[1] *= half;
353*0Sstevel@tonic-gate zz[2] *= half;
354*0Sstevel@tonic-gate ez++;
355*0Sstevel@tonic-gate }
356*0Sstevel@tonic-gate
357*0Sstevel@tonic-gate /* if exponent > 0 strip off integer bit, else denormalize */
358*0Sstevel@tonic-gate if (ez > 0) {
359*0Sstevel@tonic-gate ibit = 1;
360*0Sstevel@tonic-gate zz[0] -= one;
361*0Sstevel@tonic-gate } else {
362*0Sstevel@tonic-gate ibit = 0;
363*0Sstevel@tonic-gate if (ez > -128)
364*0Sstevel@tonic-gate u.l.hi = (unsigned)(0x3fe + ez) << 20;
365*0Sstevel@tonic-gate else
366*0Sstevel@tonic-gate u.l.hi = 0x37e00000;
367*0Sstevel@tonic-gate u.l.lo = 0;
368*0Sstevel@tonic-gate zz[0] *= u.d;
369*0Sstevel@tonic-gate zz[1] *= u.d;
370*0Sstevel@tonic-gate zz[2] *= u.d;
371*0Sstevel@tonic-gate ez = 0;
372*0Sstevel@tonic-gate }
373*0Sstevel@tonic-gate
374*0Sstevel@tonic-gate /* the first 48 bits of fraction come from zz[0] */
375*0Sstevel@tonic-gate u.d = d = two36 + zz[0];
376*0Sstevel@tonic-gate msw = u.l.lo;
377*0Sstevel@tonic-gate zz[0] -= (d - two36);
378*0Sstevel@tonic-gate
379*0Sstevel@tonic-gate u.d = d = two4 + zz[0];
380*0Sstevel@tonic-gate frac2 = u.l.lo;
381*0Sstevel@tonic-gate zz[0] -= (d - two4);
382*0Sstevel@tonic-gate
383*0Sstevel@tonic-gate /* the next 32 come from zz[0] and zz[1] */
384*0Sstevel@tonic-gate u.d = d = twom28 + (zz[0] + zz[1]);
385*0Sstevel@tonic-gate frac3 = u.l.lo;
386*0Sstevel@tonic-gate zz[0] -= (d - twom28);
387*0Sstevel@tonic-gate
388*0Sstevel@tonic-gate /* condense the remaining fraction; errors here won't matter */
389*0Sstevel@tonic-gate c = zz[0] + zz[1];
390*0Sstevel@tonic-gate zz[1] = ((zz[0] - c) + zz[1]) + zz[2];
391*0Sstevel@tonic-gate zz[0] = c;
392*0Sstevel@tonic-gate
393*0Sstevel@tonic-gate /* get the last word of fraction */
394*0Sstevel@tonic-gate u.d = d = twom60 + (zz[0] + zz[1]);
395*0Sstevel@tonic-gate frac4 = u.l.lo;
396*0Sstevel@tonic-gate zz[0] -= (d - twom60);
397*0Sstevel@tonic-gate
398*0Sstevel@tonic-gate /* keep track of what's left for rounding; note that the error */
399*0Sstevel@tonic-gate /* in computing c will be non-negative due to rounding mode */
400*0Sstevel@tonic-gate c = zz[0] + zz[1];
401*0Sstevel@tonic-gate
402*0Sstevel@tonic-gate /* get the rounding mode, fudging directed rounding modes */
403*0Sstevel@tonic-gate /* as though the result were positive */
404*0Sstevel@tonic-gate rm = fsr >> 30;
405*0Sstevel@tonic-gate if (sign)
406*0Sstevel@tonic-gate rm ^= (rm >> 1);
407*0Sstevel@tonic-gate
408*0Sstevel@tonic-gate /* round and raise exceptions */
409*0Sstevel@tonic-gate fsr &= ~FSR_CEXC;
410*0Sstevel@tonic-gate if (c != zero) {
411*0Sstevel@tonic-gate fsr |= FSR_NXC;
412*0Sstevel@tonic-gate
413*0Sstevel@tonic-gate /* decide whether to round the fraction up */
414*0Sstevel@tonic-gate if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 ||
415*0Sstevel@tonic-gate (c == twom113 && ((frac4 & 1) || (c - zz[0] != zz[1])))))) {
416*0Sstevel@tonic-gate /* round up and renormalize if necessary */
417*0Sstevel@tonic-gate if (++frac4 == 0)
418*0Sstevel@tonic-gate if (++frac3 == 0)
419*0Sstevel@tonic-gate if (++frac2 == 0)
420*0Sstevel@tonic-gate if (++msw == 0x10000) {
421*0Sstevel@tonic-gate msw = 0;
422*0Sstevel@tonic-gate ez++;
423*0Sstevel@tonic-gate }
424*0Sstevel@tonic-gate }
425*0Sstevel@tonic-gate }
426*0Sstevel@tonic-gate
427*0Sstevel@tonic-gate /* check for under/overflow */
428*0Sstevel@tonic-gate if (ez >= 0x7fff) {
429*0Sstevel@tonic-gate if (rm == FSR_RN || rm == FSR_RP) {
430*0Sstevel@tonic-gate z.l.msw = sign | 0x7fff0000;
431*0Sstevel@tonic-gate z.l.frac2 = z.l.frac3 = z.l.frac4 = 0;
432*0Sstevel@tonic-gate } else {
433*0Sstevel@tonic-gate z.l.msw = sign | 0x7ffeffff;
434*0Sstevel@tonic-gate z.l.frac2 = z.l.frac3 = z.l.frac4 = 0xffffffff;
435*0Sstevel@tonic-gate }
436*0Sstevel@tonic-gate fsr |= FSR_OFC | FSR_NXC;
437*0Sstevel@tonic-gate } else {
438*0Sstevel@tonic-gate z.l.msw = sign | (ez << 16) | msw;
439*0Sstevel@tonic-gate z.l.frac2 = frac2;
440*0Sstevel@tonic-gate z.l.frac3 = frac3;
441*0Sstevel@tonic-gate z.l.frac4 = frac4;
442*0Sstevel@tonic-gate
443*0Sstevel@tonic-gate /* !ibit => exact result was tiny before rounding, */
444*0Sstevel@tonic-gate /* t nonzero => result delivered is inexact */
445*0Sstevel@tonic-gate if (!ibit) {
446*0Sstevel@tonic-gate if (c != zero)
447*0Sstevel@tonic-gate fsr |= FSR_UFC | FSR_NXC;
448*0Sstevel@tonic-gate else if (fsr & FSR_UFM)
449*0Sstevel@tonic-gate fsr |= FSR_UFC;
450*0Sstevel@tonic-gate }
451*0Sstevel@tonic-gate }
452*0Sstevel@tonic-gate
453*0Sstevel@tonic-gate if ((fsr & FSR_CEXC) & (fsr >> 23)) {
454*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
455*0Sstevel@tonic-gate __quad_fmulq(x, y, &Z);
456*0Sstevel@tonic-gate } else {
457*0Sstevel@tonic-gate Z = z;
458*0Sstevel@tonic-gate fsr |= (fsr & 0x1f) << 5;
459*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
460*0Sstevel@tonic-gate }
461*0Sstevel@tonic-gate QUAD_RETURN(Z);
462*0Sstevel@tonic-gate }
463