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 1.0,
34*0Sstevel@tonic-gate 68719476736.0,
35*0Sstevel@tonic-gate 402653184.0,
36*0Sstevel@tonic-gate 24.0,
37*0Sstevel@tonic-gate 16.0,
38*0Sstevel@tonic-gate 3.66210937500000000000e-04,
39*0Sstevel@tonic-gate 1.52587890625000000000e-05,
40*0Sstevel@tonic-gate 1.43051147460937500000e-06,
41*0Sstevel@tonic-gate 5.96046447753906250000e-08,
42*0Sstevel@tonic-gate 3.72529029846191406250e-09,
43*0Sstevel@tonic-gate 2.18278728425502777100e-11,
44*0Sstevel@tonic-gate 8.52651282912120223045e-14,
45*0Sstevel@tonic-gate 3.55271367880050092936e-15,
46*0Sstevel@tonic-gate 1.30104260698260532081e-18,
47*0Sstevel@tonic-gate 8.67361737988403547206e-19,
48*0Sstevel@tonic-gate 2.16840434497100886801e-19,
49*0Sstevel@tonic-gate 3.17637355220362627151e-22,
50*0Sstevel@tonic-gate 7.75481824268463445192e-26,
51*0Sstevel@tonic-gate 4.62223186652936604733e-33,
52*0Sstevel@tonic-gate 9.62964972193617926528e-35,
53*0Sstevel@tonic-gate 4.70197740328915003187e-38,
54*0Sstevel@tonic-gate 2.75506488473973634680e-40,
55*0Sstevel@tonic-gate };
56*0Sstevel@tonic-gate
57*0Sstevel@tonic-gate #define zero C[0]
58*0Sstevel@tonic-gate #define one C[1]
59*0Sstevel@tonic-gate #define two36 C[2]
60*0Sstevel@tonic-gate #define three2p27 C[3]
61*0Sstevel@tonic-gate #define three2p3 C[4]
62*0Sstevel@tonic-gate #define two4 C[5]
63*0Sstevel@tonic-gate #define three2m13 C[6]
64*0Sstevel@tonic-gate #define twom16 C[7]
65*0Sstevel@tonic-gate #define three2m21 C[8]
66*0Sstevel@tonic-gate #define twom24 C[9]
67*0Sstevel@tonic-gate #define twom28 C[10]
68*0Sstevel@tonic-gate #define three2m37 C[11]
69*0Sstevel@tonic-gate #define three2m45 C[12]
70*0Sstevel@tonic-gate #define twom48 C[13]
71*0Sstevel@tonic-gate #define three2m61 C[14]
72*0Sstevel@tonic-gate #define twom60 C[15]
73*0Sstevel@tonic-gate #define twom62 C[16]
74*0Sstevel@tonic-gate #define three2m73 C[17]
75*0Sstevel@tonic-gate #define three2m85 C[18]
76*0Sstevel@tonic-gate #define three2m109 C[19]
77*0Sstevel@tonic-gate #define twom113 C[20]
78*0Sstevel@tonic-gate #define twom124 C[21]
79*0Sstevel@tonic-gate #define three2m133 C[22]
80*0Sstevel@tonic-gate
81*0Sstevel@tonic-gate static const unsigned int
82*0Sstevel@tonic-gate fsr_re = 0x00000000u,
83*0Sstevel@tonic-gate fsr_rn = 0xc0000000u;
84*0Sstevel@tonic-gate
85*0Sstevel@tonic-gate #ifdef __sparcv9
86*0Sstevel@tonic-gate
87*0Sstevel@tonic-gate /*
88*0Sstevel@tonic-gate * _Qp_div(pz, x, y) sets *pz = *x / *y.
89*0Sstevel@tonic-gate */
90*0Sstevel@tonic-gate void
_Qp_div(union longdouble * pz,const union longdouble * x,const union longdouble * y)91*0Sstevel@tonic-gate _Qp_div(union longdouble *pz, const union longdouble *x,
92*0Sstevel@tonic-gate const union longdouble *y)
93*0Sstevel@tonic-gate
94*0Sstevel@tonic-gate #else
95*0Sstevel@tonic-gate
96*0Sstevel@tonic-gate /*
97*0Sstevel@tonic-gate * _Q_div(x, y) returns *x / *y.
98*0Sstevel@tonic-gate */
99*0Sstevel@tonic-gate union longdouble
100*0Sstevel@tonic-gate _Q_div(const union longdouble *x, const union longdouble *y)
101*0Sstevel@tonic-gate
102*0Sstevel@tonic-gate #endif /* __sparcv9 */
103*0Sstevel@tonic-gate
104*0Sstevel@tonic-gate {
105*0Sstevel@tonic-gate union longdouble z;
106*0Sstevel@tonic-gate union xdouble u;
107*0Sstevel@tonic-gate double c, d, ry, xx[4], yy[5], zz[5];
108*0Sstevel@tonic-gate unsigned int xm, ym, fsr, lx, ly, wx[3], wy[3];
109*0Sstevel@tonic-gate unsigned int msw, frac2, frac3, frac4, rm;
110*0Sstevel@tonic-gate int ibit, ex, ey, ez, sign;
111*0Sstevel@tonic-gate
112*0Sstevel@tonic-gate xm = x->l.msw & 0x7fffffff;
113*0Sstevel@tonic-gate ym = y->l.msw & 0x7fffffff;
114*0Sstevel@tonic-gate sign = (x->l.msw ^ y->l.msw) & ~0x7fffffff;
115*0Sstevel@tonic-gate
116*0Sstevel@tonic-gate __quad_getfsrp(&fsr);
117*0Sstevel@tonic-gate
118*0Sstevel@tonic-gate /* handle nan and inf cases */
119*0Sstevel@tonic-gate if (xm >= 0x7fff0000 || ym >= 0x7fff0000) {
120*0Sstevel@tonic-gate /* handle nan cases according to V9 app. B */
121*0Sstevel@tonic-gate if (QUAD_ISNAN(*y)) {
122*0Sstevel@tonic-gate if (!(y->l.msw & 0x8000)) {
123*0Sstevel@tonic-gate /* snan, signal invalid */
124*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
125*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
126*0Sstevel@tonic-gate } else {
127*0Sstevel@tonic-gate Z = *y;
128*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
129*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
130*0Sstevel@tonic-gate FSR_NVC;
131*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
132*0Sstevel@tonic-gate }
133*0Sstevel@tonic-gate QUAD_RETURN(Z);
134*0Sstevel@tonic-gate } else if (QUAD_ISNAN(*x) && !(x->l.msw & 0x8000)) {
135*0Sstevel@tonic-gate /* snan, signal invalid */
136*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
137*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
138*0Sstevel@tonic-gate } else {
139*0Sstevel@tonic-gate Z = *x;
140*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
141*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
142*0Sstevel@tonic-gate FSR_NVC;
143*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
144*0Sstevel@tonic-gate }
145*0Sstevel@tonic-gate QUAD_RETURN(Z);
146*0Sstevel@tonic-gate }
147*0Sstevel@tonic-gate Z = *y;
148*0Sstevel@tonic-gate QUAD_RETURN(Z);
149*0Sstevel@tonic-gate }
150*0Sstevel@tonic-gate if (QUAD_ISNAN(*x)) {
151*0Sstevel@tonic-gate if (!(x->l.msw & 0x8000)) {
152*0Sstevel@tonic-gate /* snan, signal invalid */
153*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
154*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
155*0Sstevel@tonic-gate } else {
156*0Sstevel@tonic-gate Z = *x;
157*0Sstevel@tonic-gate Z.l.msw |= 0x8000;
158*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
159*0Sstevel@tonic-gate FSR_NVC;
160*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
161*0Sstevel@tonic-gate }
162*0Sstevel@tonic-gate QUAD_RETURN(Z);
163*0Sstevel@tonic-gate }
164*0Sstevel@tonic-gate Z = *x;
165*0Sstevel@tonic-gate QUAD_RETURN(Z);
166*0Sstevel@tonic-gate }
167*0Sstevel@tonic-gate if (xm == 0x7fff0000) {
168*0Sstevel@tonic-gate /* x is inf */
169*0Sstevel@tonic-gate if (ym == 0x7fff0000) {
170*0Sstevel@tonic-gate /* inf / inf, signal invalid */
171*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
172*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
173*0Sstevel@tonic-gate } else {
174*0Sstevel@tonic-gate Z.l.msw = 0x7fffffff;
175*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 =
176*0Sstevel@tonic-gate Z.l.frac4 = 0xffffffff;
177*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
178*0Sstevel@tonic-gate 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 /* y is inf */
189*0Sstevel@tonic-gate Z.l.msw = sign;
190*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
191*0Sstevel@tonic-gate QUAD_RETURN(Z);
192*0Sstevel@tonic-gate }
193*0Sstevel@tonic-gate
194*0Sstevel@tonic-gate /* handle zero cases */
195*0Sstevel@tonic-gate if (xm == 0 || ym == 0) {
196*0Sstevel@tonic-gate if (QUAD_ISZERO(*x)) {
197*0Sstevel@tonic-gate if (QUAD_ISZERO(*y)) {
198*0Sstevel@tonic-gate /* zero / zero, signal invalid */
199*0Sstevel@tonic-gate if (fsr & FSR_NVM) {
200*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
201*0Sstevel@tonic-gate } else {
202*0Sstevel@tonic-gate Z.l.msw = 0x7fffffff;
203*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 =
204*0Sstevel@tonic-gate Z.l.frac4 = 0xffffffff;
205*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
206*0Sstevel@tonic-gate FSR_NVC;
207*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
208*0Sstevel@tonic-gate }
209*0Sstevel@tonic-gate QUAD_RETURN(Z);
210*0Sstevel@tonic-gate }
211*0Sstevel@tonic-gate /* zero / nonzero, return zero */
212*0Sstevel@tonic-gate Z.l.msw = sign;
213*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
214*0Sstevel@tonic-gate QUAD_RETURN(Z);
215*0Sstevel@tonic-gate }
216*0Sstevel@tonic-gate if (QUAD_ISZERO(*y)) {
217*0Sstevel@tonic-gate /* nonzero / zero, signal zero divide */
218*0Sstevel@tonic-gate if (fsr & FSR_DZM) {
219*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
220*0Sstevel@tonic-gate } else {
221*0Sstevel@tonic-gate Z.l.msw = sign | 0x7fff0000;
222*0Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0;
223*0Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_DZA | FSR_DZC;
224*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
225*0Sstevel@tonic-gate }
226*0Sstevel@tonic-gate QUAD_RETURN(Z);
227*0Sstevel@tonic-gate }
228*0Sstevel@tonic-gate }
229*0Sstevel@tonic-gate
230*0Sstevel@tonic-gate /* now x and y are finite, nonzero */
231*0Sstevel@tonic-gate __quad_setfsrp(&fsr_re);
232*0Sstevel@tonic-gate
233*0Sstevel@tonic-gate /* get their normalized significands and exponents */
234*0Sstevel@tonic-gate ex = (int)(xm >> 16);
235*0Sstevel@tonic-gate lx = xm & 0xffff;
236*0Sstevel@tonic-gate if (ex) {
237*0Sstevel@tonic-gate lx |= 0x10000;
238*0Sstevel@tonic-gate wx[0] = x->l.frac2;
239*0Sstevel@tonic-gate wx[1] = x->l.frac3;
240*0Sstevel@tonic-gate wx[2] = x->l.frac4;
241*0Sstevel@tonic-gate } else {
242*0Sstevel@tonic-gate if (lx | (x->l.frac2 & 0xfffe0000)) {
243*0Sstevel@tonic-gate wx[0] = x->l.frac2;
244*0Sstevel@tonic-gate wx[1] = x->l.frac3;
245*0Sstevel@tonic-gate wx[2] = x->l.frac4;
246*0Sstevel@tonic-gate ex = 1;
247*0Sstevel@tonic-gate } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
248*0Sstevel@tonic-gate lx = x->l.frac2;
249*0Sstevel@tonic-gate wx[0] = x->l.frac3;
250*0Sstevel@tonic-gate wx[1] = x->l.frac4;
251*0Sstevel@tonic-gate wx[2] = 0;
252*0Sstevel@tonic-gate ex = -31;
253*0Sstevel@tonic-gate } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
254*0Sstevel@tonic-gate lx = x->l.frac3;
255*0Sstevel@tonic-gate wx[0] = x->l.frac4;
256*0Sstevel@tonic-gate wx[1] = wx[2] = 0;
257*0Sstevel@tonic-gate ex = -63;
258*0Sstevel@tonic-gate } else {
259*0Sstevel@tonic-gate lx = x->l.frac4;
260*0Sstevel@tonic-gate wx[0] = wx[1] = wx[2] = 0;
261*0Sstevel@tonic-gate ex = -95;
262*0Sstevel@tonic-gate }
263*0Sstevel@tonic-gate while ((lx & 0x10000) == 0) {
264*0Sstevel@tonic-gate lx = (lx << 1) | (wx[0] >> 31);
265*0Sstevel@tonic-gate wx[0] = (wx[0] << 1) | (wx[1] >> 31);
266*0Sstevel@tonic-gate wx[1] = (wx[1] << 1) | (wx[2] >> 31);
267*0Sstevel@tonic-gate wx[2] <<= 1;
268*0Sstevel@tonic-gate ex--;
269*0Sstevel@tonic-gate }
270*0Sstevel@tonic-gate }
271*0Sstevel@tonic-gate ez = ex;
272*0Sstevel@tonic-gate
273*0Sstevel@tonic-gate ey = (int)(ym >> 16);
274*0Sstevel@tonic-gate ly = ym & 0xffff;
275*0Sstevel@tonic-gate if (ey) {
276*0Sstevel@tonic-gate ly |= 0x10000;
277*0Sstevel@tonic-gate wy[0] = y->l.frac2;
278*0Sstevel@tonic-gate wy[1] = y->l.frac3;
279*0Sstevel@tonic-gate wy[2] = y->l.frac4;
280*0Sstevel@tonic-gate } else {
281*0Sstevel@tonic-gate if (ly | (y->l.frac2 & 0xfffe0000)) {
282*0Sstevel@tonic-gate wy[0] = y->l.frac2;
283*0Sstevel@tonic-gate wy[1] = y->l.frac3;
284*0Sstevel@tonic-gate wy[2] = y->l.frac4;
285*0Sstevel@tonic-gate ey = 1;
286*0Sstevel@tonic-gate } else if (y->l.frac2 | (y->l.frac3 & 0xfffe0000)) {
287*0Sstevel@tonic-gate ly = y->l.frac2;
288*0Sstevel@tonic-gate wy[0] = y->l.frac3;
289*0Sstevel@tonic-gate wy[1] = y->l.frac4;
290*0Sstevel@tonic-gate wy[2] = 0;
291*0Sstevel@tonic-gate ey = -31;
292*0Sstevel@tonic-gate } else if (y->l.frac3 | (y->l.frac4 & 0xfffe0000)) {
293*0Sstevel@tonic-gate ly = y->l.frac3;
294*0Sstevel@tonic-gate wy[0] = y->l.frac4;
295*0Sstevel@tonic-gate wy[1] = wy[2] = 0;
296*0Sstevel@tonic-gate ey = -63;
297*0Sstevel@tonic-gate } else {
298*0Sstevel@tonic-gate ly = y->l.frac4;
299*0Sstevel@tonic-gate wy[0] = wy[1] = wy[2] = 0;
300*0Sstevel@tonic-gate ey = -95;
301*0Sstevel@tonic-gate }
302*0Sstevel@tonic-gate while ((ly & 0x10000) == 0) {
303*0Sstevel@tonic-gate ly = (ly << 1) | (wy[0] >> 31);
304*0Sstevel@tonic-gate wy[0] = (wy[0] << 1) | (wy[1] >> 31);
305*0Sstevel@tonic-gate wy[1] = (wy[1] << 1) | (wy[2] >> 31);
306*0Sstevel@tonic-gate wy[2] <<= 1;
307*0Sstevel@tonic-gate ey--;
308*0Sstevel@tonic-gate }
309*0Sstevel@tonic-gate }
310*0Sstevel@tonic-gate ez -= ey - 0x3fff;
311*0Sstevel@tonic-gate
312*0Sstevel@tonic-gate /* extract the significands into doubles */
313*0Sstevel@tonic-gate c = twom16;
314*0Sstevel@tonic-gate xx[0] = (double)((int)lx) * c;
315*0Sstevel@tonic-gate yy[0] = (double)((int)ly) * c;
316*0Sstevel@tonic-gate
317*0Sstevel@tonic-gate c *= twom24;
318*0Sstevel@tonic-gate xx[0] += (double)((int)(wx[0] >> 8)) * c;
319*0Sstevel@tonic-gate yy[1] = (double)((int)(wy[0] >> 8)) * c;
320*0Sstevel@tonic-gate
321*0Sstevel@tonic-gate c *= twom24;
322*0Sstevel@tonic-gate xx[1] = (double)((int)(((wx[0] << 16) |
323*0Sstevel@tonic-gate (wx[1] >> 16)) & 0xffffff)) * c;
324*0Sstevel@tonic-gate yy[2] = (double)((int)(((wy[0] << 16) |
325*0Sstevel@tonic-gate (wy[1] >> 16)) & 0xffffff)) * c;
326*0Sstevel@tonic-gate
327*0Sstevel@tonic-gate c *= twom24;
328*0Sstevel@tonic-gate xx[2] = (double)((int)(((wx[1] << 8) |
329*0Sstevel@tonic-gate (wx[2] >> 24)) & 0xffffff)) * c;
330*0Sstevel@tonic-gate yy[3] = (double)((int)(((wy[1] << 8) |
331*0Sstevel@tonic-gate (wy[2] >> 24)) & 0xffffff)) * c;
332*0Sstevel@tonic-gate
333*0Sstevel@tonic-gate c *= twom24;
334*0Sstevel@tonic-gate xx[3] = (double)((int)(wx[2] & 0xffffff)) * c;
335*0Sstevel@tonic-gate yy[4] = (double)((int)(wy[2] & 0xffffff)) * c;
336*0Sstevel@tonic-gate
337*0Sstevel@tonic-gate /* approximate the reciprocal of y */
338*0Sstevel@tonic-gate ry = one / ((yy[0] + yy[1]) + yy[2]);
339*0Sstevel@tonic-gate
340*0Sstevel@tonic-gate /* compute the first five "digits" of the quotient */
341*0Sstevel@tonic-gate zz[0] = (ry * (xx[0] + xx[1]) + three2p27) - three2p27;
342*0Sstevel@tonic-gate xx[0] = ((xx[0] - zz[0] * yy[0]) - zz[0] * yy[1]) + xx[1];
343*0Sstevel@tonic-gate d = zz[0] * yy[2];
344*0Sstevel@tonic-gate c = (d + three2m13) - three2m13;
345*0Sstevel@tonic-gate xx[0] -= c;
346*0Sstevel@tonic-gate xx[1] = xx[2] - (d - c);
347*0Sstevel@tonic-gate d = zz[0] * yy[3];
348*0Sstevel@tonic-gate c = (d + three2m37) - three2m37;
349*0Sstevel@tonic-gate xx[1] -= c;
350*0Sstevel@tonic-gate xx[2] = xx[3] - (d - c);
351*0Sstevel@tonic-gate d = zz[0] * yy[4];
352*0Sstevel@tonic-gate c = (d + three2m61) - three2m61;
353*0Sstevel@tonic-gate xx[2] -= c;
354*0Sstevel@tonic-gate xx[3] = c - d;
355*0Sstevel@tonic-gate
356*0Sstevel@tonic-gate zz[1] = (ry * (xx[0] + xx[1]) + three2p3) - three2p3;
357*0Sstevel@tonic-gate xx[0] = ((xx[0] - zz[1] * yy[0]) - zz[1] * yy[1]) + xx[1];
358*0Sstevel@tonic-gate d = zz[1] * yy[2];
359*0Sstevel@tonic-gate c = (d + three2m37) - three2m37;
360*0Sstevel@tonic-gate xx[0] -= c;
361*0Sstevel@tonic-gate xx[1] = xx[2] - (d - c);
362*0Sstevel@tonic-gate d = zz[1] * yy[3];
363*0Sstevel@tonic-gate c = (d + three2m61) - three2m61;
364*0Sstevel@tonic-gate xx[1] -= c;
365*0Sstevel@tonic-gate xx[2] = xx[3] - (d - c);
366*0Sstevel@tonic-gate d = zz[1] * yy[4];
367*0Sstevel@tonic-gate c = (d + three2m85) - three2m85;
368*0Sstevel@tonic-gate xx[2] -= c;
369*0Sstevel@tonic-gate xx[3] = c - d;
370*0Sstevel@tonic-gate
371*0Sstevel@tonic-gate zz[2] = (ry * (xx[0] + xx[1]) + three2m21) - three2m21;
372*0Sstevel@tonic-gate xx[0] = ((xx[0] - zz[2] * yy[0]) - zz[2] * yy[1]) + xx[1];
373*0Sstevel@tonic-gate d = zz[2] * yy[2];
374*0Sstevel@tonic-gate c = (d + three2m61) - three2m61;
375*0Sstevel@tonic-gate xx[0] -= c;
376*0Sstevel@tonic-gate xx[1] = xx[2] - (d - c);
377*0Sstevel@tonic-gate d = zz[2] * yy[3];
378*0Sstevel@tonic-gate c = (d + three2m85) - three2m85;
379*0Sstevel@tonic-gate xx[1] -= c;
380*0Sstevel@tonic-gate xx[2] = xx[3] - (d - c);
381*0Sstevel@tonic-gate d = zz[2] * yy[4];
382*0Sstevel@tonic-gate c = (d + three2m109) - three2m109;
383*0Sstevel@tonic-gate xx[2] -= c;
384*0Sstevel@tonic-gate xx[3] = c - d;
385*0Sstevel@tonic-gate
386*0Sstevel@tonic-gate zz[3] = (ry * (xx[0] + xx[1]) + three2m45) - three2m45;
387*0Sstevel@tonic-gate xx[0] = ((xx[0] - zz[3] * yy[0]) - zz[3] * yy[1]) + xx[1];
388*0Sstevel@tonic-gate d = zz[3] * yy[2];
389*0Sstevel@tonic-gate c = (d + three2m85) - three2m85;
390*0Sstevel@tonic-gate xx[0] -= c;
391*0Sstevel@tonic-gate xx[1] = xx[2] - (d - c);
392*0Sstevel@tonic-gate d = zz[3] * yy[3];
393*0Sstevel@tonic-gate c = (d + three2m109) - three2m109;
394*0Sstevel@tonic-gate xx[1] -= c;
395*0Sstevel@tonic-gate xx[2] = xx[3] - (d - c);
396*0Sstevel@tonic-gate d = zz[3] * yy[4];
397*0Sstevel@tonic-gate c = (d + three2m133) - three2m133;
398*0Sstevel@tonic-gate xx[2] -= c;
399*0Sstevel@tonic-gate xx[3] = c - d;
400*0Sstevel@tonic-gate
401*0Sstevel@tonic-gate zz[4] = (ry * (xx[0] + xx[1]) + three2m73) - three2m73;
402*0Sstevel@tonic-gate
403*0Sstevel@tonic-gate /* reduce to three doubles, making sure zz[1] is positive */
404*0Sstevel@tonic-gate zz[0] += zz[1] - twom48;
405*0Sstevel@tonic-gate zz[1] = twom48 + zz[2] + zz[3];
406*0Sstevel@tonic-gate zz[2] = zz[4];
407*0Sstevel@tonic-gate
408*0Sstevel@tonic-gate /* if the third term might lie on a rounding boundary, perturb it */
409*0Sstevel@tonic-gate if (zz[2] == (twom62 + zz[2]) - twom62) {
410*0Sstevel@tonic-gate /* here we just need to get the sign of the remainder */
411*0Sstevel@tonic-gate c = (((((xx[0] - zz[4] * yy[0]) - zz[4] * yy[1]) + xx[1]) +
412*0Sstevel@tonic-gate (xx[2] - zz[4] * yy[2])) + (xx[3] - zz[4] * yy[3]))
413*0Sstevel@tonic-gate - zz[4] * yy[4];
414*0Sstevel@tonic-gate if (c < zero)
415*0Sstevel@tonic-gate zz[2] -= twom124;
416*0Sstevel@tonic-gate else if (c > zero)
417*0Sstevel@tonic-gate zz[2] += twom124;
418*0Sstevel@tonic-gate }
419*0Sstevel@tonic-gate
420*0Sstevel@tonic-gate /*
421*0Sstevel@tonic-gate * propagate carries/borrows, using round-to-negative-infinity mode
422*0Sstevel@tonic-gate * to make all terms nonnegative (note that we can't encounter a
423*0Sstevel@tonic-gate * borrow so large that the roundoff is unrepresentable because
424*0Sstevel@tonic-gate * we took care to make zz[1] positive above)
425*0Sstevel@tonic-gate */
426*0Sstevel@tonic-gate __quad_setfsrp(&fsr_rn);
427*0Sstevel@tonic-gate c = zz[1] + zz[2];
428*0Sstevel@tonic-gate zz[2] += (zz[1] - c);
429*0Sstevel@tonic-gate zz[1] = c;
430*0Sstevel@tonic-gate c = zz[0] + zz[1];
431*0Sstevel@tonic-gate zz[1] += (zz[0] - c);
432*0Sstevel@tonic-gate zz[0] = c;
433*0Sstevel@tonic-gate
434*0Sstevel@tonic-gate /* check for borrow */
435*0Sstevel@tonic-gate if (c < one) {
436*0Sstevel@tonic-gate /* postnormalize */
437*0Sstevel@tonic-gate zz[0] += zz[0];
438*0Sstevel@tonic-gate zz[1] += zz[1];
439*0Sstevel@tonic-gate zz[2] += zz[2];
440*0Sstevel@tonic-gate ez--;
441*0Sstevel@tonic-gate }
442*0Sstevel@tonic-gate
443*0Sstevel@tonic-gate /* if exponent > 0 strip off integer bit, else denormalize */
444*0Sstevel@tonic-gate if (ez > 0) {
445*0Sstevel@tonic-gate ibit = 1;
446*0Sstevel@tonic-gate zz[0] -= one;
447*0Sstevel@tonic-gate } else {
448*0Sstevel@tonic-gate ibit = 0;
449*0Sstevel@tonic-gate if (ez > -128)
450*0Sstevel@tonic-gate u.l.hi = (unsigned int)(0x3fe + ez) << 20;
451*0Sstevel@tonic-gate else
452*0Sstevel@tonic-gate u.l.hi = 0x37e00000;
453*0Sstevel@tonic-gate u.l.lo = 0;
454*0Sstevel@tonic-gate zz[0] *= u.d;
455*0Sstevel@tonic-gate zz[1] *= u.d;
456*0Sstevel@tonic-gate zz[2] *= u.d;
457*0Sstevel@tonic-gate ez = 0;
458*0Sstevel@tonic-gate }
459*0Sstevel@tonic-gate
460*0Sstevel@tonic-gate /* the first 48 bits of fraction come from zz[0] */
461*0Sstevel@tonic-gate u.d = d = two36 + zz[0];
462*0Sstevel@tonic-gate msw = u.l.lo;
463*0Sstevel@tonic-gate zz[0] -= (d - two36);
464*0Sstevel@tonic-gate
465*0Sstevel@tonic-gate u.d = d = two4 + zz[0];
466*0Sstevel@tonic-gate frac2 = u.l.lo;
467*0Sstevel@tonic-gate zz[0] -= (d - two4);
468*0Sstevel@tonic-gate
469*0Sstevel@tonic-gate /* the next 32 come from zz[0] and zz[1] */
470*0Sstevel@tonic-gate u.d = d = twom28 + (zz[0] + zz[1]);
471*0Sstevel@tonic-gate frac3 = u.l.lo;
472*0Sstevel@tonic-gate zz[0] -= (d - twom28);
473*0Sstevel@tonic-gate
474*0Sstevel@tonic-gate /* condense the remaining fraction; errors here won't matter */
475*0Sstevel@tonic-gate c = zz[0] + zz[1];
476*0Sstevel@tonic-gate zz[1] = ((zz[0] - c) + zz[1]) + zz[2];
477*0Sstevel@tonic-gate zz[0] = c;
478*0Sstevel@tonic-gate
479*0Sstevel@tonic-gate /* get the last word of fraction */
480*0Sstevel@tonic-gate u.d = d = twom60 + (zz[0] + zz[1]);
481*0Sstevel@tonic-gate frac4 = u.l.lo;
482*0Sstevel@tonic-gate zz[0] -= (d - twom60);
483*0Sstevel@tonic-gate
484*0Sstevel@tonic-gate /* keep track of what's left for rounding; note that the error */
485*0Sstevel@tonic-gate /* in computing c will be non-negative due to rounding mode */
486*0Sstevel@tonic-gate c = zz[0] + zz[1];
487*0Sstevel@tonic-gate
488*0Sstevel@tonic-gate /* get the rounding mode, fudging directed rounding modes */
489*0Sstevel@tonic-gate /* as though the result were positive */
490*0Sstevel@tonic-gate rm = fsr >> 30;
491*0Sstevel@tonic-gate if (sign)
492*0Sstevel@tonic-gate rm ^= (rm >> 1);
493*0Sstevel@tonic-gate
494*0Sstevel@tonic-gate /* round and raise exceptions */
495*0Sstevel@tonic-gate fsr &= ~FSR_CEXC;
496*0Sstevel@tonic-gate if (c != zero) {
497*0Sstevel@tonic-gate fsr |= FSR_NXC;
498*0Sstevel@tonic-gate
499*0Sstevel@tonic-gate /* decide whether to round the fraction up */
500*0Sstevel@tonic-gate if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 ||
501*0Sstevel@tonic-gate (c == twom113 && ((frac4 & 1) || (c - zz[0] !=
502*0Sstevel@tonic-gate zz[1])))))) {
503*0Sstevel@tonic-gate /* round up and renormalize if necessary */
504*0Sstevel@tonic-gate if (++frac4 == 0)
505*0Sstevel@tonic-gate if (++frac3 == 0)
506*0Sstevel@tonic-gate if (++frac2 == 0)
507*0Sstevel@tonic-gate if (++msw == 0x10000) {
508*0Sstevel@tonic-gate msw = 0;
509*0Sstevel@tonic-gate ez++;
510*0Sstevel@tonic-gate }
511*0Sstevel@tonic-gate }
512*0Sstevel@tonic-gate }
513*0Sstevel@tonic-gate
514*0Sstevel@tonic-gate /* check for under/overflow */
515*0Sstevel@tonic-gate if (ez >= 0x7fff) {
516*0Sstevel@tonic-gate if (rm == FSR_RN || rm == FSR_RP) {
517*0Sstevel@tonic-gate z.l.msw = sign | 0x7fff0000;
518*0Sstevel@tonic-gate z.l.frac2 = z.l.frac3 = z.l.frac4 = 0;
519*0Sstevel@tonic-gate } else {
520*0Sstevel@tonic-gate z.l.msw = sign | 0x7ffeffff;
521*0Sstevel@tonic-gate z.l.frac2 = z.l.frac3 = z.l.frac4 = 0xffffffff;
522*0Sstevel@tonic-gate }
523*0Sstevel@tonic-gate fsr |= FSR_OFC | FSR_NXC;
524*0Sstevel@tonic-gate } else {
525*0Sstevel@tonic-gate z.l.msw = sign | (ez << 16) | msw;
526*0Sstevel@tonic-gate z.l.frac2 = frac2;
527*0Sstevel@tonic-gate z.l.frac3 = frac3;
528*0Sstevel@tonic-gate z.l.frac4 = frac4;
529*0Sstevel@tonic-gate
530*0Sstevel@tonic-gate /* !ibit => exact result was tiny before rounding, */
531*0Sstevel@tonic-gate /* t nonzero => result delivered is inexact */
532*0Sstevel@tonic-gate if (!ibit) {
533*0Sstevel@tonic-gate if (c != zero)
534*0Sstevel@tonic-gate fsr |= FSR_UFC | FSR_NXC;
535*0Sstevel@tonic-gate else if (fsr & FSR_UFM)
536*0Sstevel@tonic-gate fsr |= FSR_UFC;
537*0Sstevel@tonic-gate }
538*0Sstevel@tonic-gate }
539*0Sstevel@tonic-gate
540*0Sstevel@tonic-gate if ((fsr & FSR_CEXC) & (fsr >> 23)) {
541*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
542*0Sstevel@tonic-gate __quad_fdivq(x, y, &Z);
543*0Sstevel@tonic-gate } else {
544*0Sstevel@tonic-gate Z = z;
545*0Sstevel@tonic-gate fsr |= (fsr & 0x1f) << 5;
546*0Sstevel@tonic-gate __quad_setfsrp(&fsr);
547*0Sstevel@tonic-gate }
548*0Sstevel@tonic-gate QUAD_RETURN(Z);
549*0Sstevel@tonic-gate }
550