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