xref: /illumos-gate/usr/src/lib/libm/common/m9x/fmal.c (revision 25c28e83)
1*25c28e83SPiotr Jasiukajtis /*
2*25c28e83SPiotr Jasiukajtis  * CDDL HEADER START
3*25c28e83SPiotr Jasiukajtis  *
4*25c28e83SPiotr Jasiukajtis  * The contents of this file are subject to the terms of the
5*25c28e83SPiotr Jasiukajtis  * Common Development and Distribution License (the "License").
6*25c28e83SPiotr Jasiukajtis  * You may not use this file except in compliance with the License.
7*25c28e83SPiotr Jasiukajtis  *
8*25c28e83SPiotr Jasiukajtis  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9*25c28e83SPiotr Jasiukajtis  * or http://www.opensolaris.org/os/licensing.
10*25c28e83SPiotr Jasiukajtis  * See the License for the specific language governing permissions
11*25c28e83SPiotr Jasiukajtis  * and limitations under the License.
12*25c28e83SPiotr Jasiukajtis  *
13*25c28e83SPiotr Jasiukajtis  * When distributing Covered Code, include this CDDL HEADER in each
14*25c28e83SPiotr Jasiukajtis  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15*25c28e83SPiotr Jasiukajtis  * If applicable, add the following below this CDDL HEADER, with the
16*25c28e83SPiotr Jasiukajtis  * fields enclosed by brackets "[]" replaced with your own identifying
17*25c28e83SPiotr Jasiukajtis  * information: Portions Copyright [yyyy] [name of copyright owner]
18*25c28e83SPiotr Jasiukajtis  *
19*25c28e83SPiotr Jasiukajtis  * CDDL HEADER END
20*25c28e83SPiotr Jasiukajtis  */
21*25c28e83SPiotr Jasiukajtis 
22*25c28e83SPiotr Jasiukajtis /*
23*25c28e83SPiotr Jasiukajtis  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24*25c28e83SPiotr Jasiukajtis  */
25*25c28e83SPiotr Jasiukajtis /*
26*25c28e83SPiotr Jasiukajtis  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27*25c28e83SPiotr Jasiukajtis  * Use is subject to license terms.
28*25c28e83SPiotr Jasiukajtis  */
29*25c28e83SPiotr Jasiukajtis 
30*25c28e83SPiotr Jasiukajtis #if defined(ELFOBJ)
31*25c28e83SPiotr Jasiukajtis #pragma weak fmal = __fmal
32*25c28e83SPiotr Jasiukajtis #endif
33*25c28e83SPiotr Jasiukajtis 
34*25c28e83SPiotr Jasiukajtis #include "libm.h"
35*25c28e83SPiotr Jasiukajtis #include "fma.h"
36*25c28e83SPiotr Jasiukajtis #include "fenv_inlines.h"
37*25c28e83SPiotr Jasiukajtis 
38*25c28e83SPiotr Jasiukajtis #if defined(__sparc)
39*25c28e83SPiotr Jasiukajtis 
40*25c28e83SPiotr Jasiukajtis static const union {
41*25c28e83SPiotr Jasiukajtis 	unsigned i[2];
42*25c28e83SPiotr Jasiukajtis 	double d;
43*25c28e83SPiotr Jasiukajtis } C[] = {
44*25c28e83SPiotr Jasiukajtis 	{ 0x3fe00000u, 0 },
45*25c28e83SPiotr Jasiukajtis 	{ 0x40000000u, 0 },
46*25c28e83SPiotr Jasiukajtis 	{ 0x3ef00000u, 0 },
47*25c28e83SPiotr Jasiukajtis 	{ 0x3e700000u, 0 },
48*25c28e83SPiotr Jasiukajtis 	{ 0x41300000u, 0 },
49*25c28e83SPiotr Jasiukajtis 	{ 0x3e300000u, 0 },
50*25c28e83SPiotr Jasiukajtis 	{ 0x3b300000u, 0 },
51*25c28e83SPiotr Jasiukajtis 	{ 0x38300000u, 0 },
52*25c28e83SPiotr Jasiukajtis 	{ 0x42300000u, 0 },
53*25c28e83SPiotr Jasiukajtis 	{ 0x3df00000u, 0 },
54*25c28e83SPiotr Jasiukajtis 	{ 0x7fe00000u, 0 },
55*25c28e83SPiotr Jasiukajtis 	{ 0x00100000u, 0 },
56*25c28e83SPiotr Jasiukajtis 	{ 0x00100001u, 0 },
57*25c28e83SPiotr Jasiukajtis 	{ 0, 0 },
58*25c28e83SPiotr Jasiukajtis 	{ 0x7ff00000u, 0 },
59*25c28e83SPiotr Jasiukajtis 	{ 0x7ff00001u, 0 }
60*25c28e83SPiotr Jasiukajtis };
61*25c28e83SPiotr Jasiukajtis 
62*25c28e83SPiotr Jasiukajtis #define	half	C[0].d
63*25c28e83SPiotr Jasiukajtis #define	two	C[1].d
64*25c28e83SPiotr Jasiukajtis #define	twom16	C[2].d
65*25c28e83SPiotr Jasiukajtis #define	twom24	C[3].d
66*25c28e83SPiotr Jasiukajtis #define	two20	C[4].d
67*25c28e83SPiotr Jasiukajtis #define	twom28	C[5].d
68*25c28e83SPiotr Jasiukajtis #define	twom76	C[6].d
69*25c28e83SPiotr Jasiukajtis #define	twom124	C[7].d
70*25c28e83SPiotr Jasiukajtis #define	two36	C[8].d
71*25c28e83SPiotr Jasiukajtis #define	twom32	C[9].d
72*25c28e83SPiotr Jasiukajtis #define	huge	C[10].d
73*25c28e83SPiotr Jasiukajtis #define	tiny	C[11].d
74*25c28e83SPiotr Jasiukajtis #define	tiny2	C[12].d
75*25c28e83SPiotr Jasiukajtis #define	zero	C[13].d
76*25c28e83SPiotr Jasiukajtis #define	inf	C[14].d
77*25c28e83SPiotr Jasiukajtis #define	snan	C[15].d
78*25c28e83SPiotr Jasiukajtis 
79*25c28e83SPiotr Jasiukajtis static const unsigned int fsr_rm = 0xc0000000u;
80*25c28e83SPiotr Jasiukajtis 
81*25c28e83SPiotr Jasiukajtis /*
82*25c28e83SPiotr Jasiukajtis  * fmal for SPARC: 128-bit quad precision, big-endian
83*25c28e83SPiotr Jasiukajtis  */
84*25c28e83SPiotr Jasiukajtis long double
85*25c28e83SPiotr Jasiukajtis __fmal(long double x, long double y, long double z) {
86*25c28e83SPiotr Jasiukajtis 	union {
87*25c28e83SPiotr Jasiukajtis 		unsigned int i[4];
88*25c28e83SPiotr Jasiukajtis 		long double q;
89*25c28e83SPiotr Jasiukajtis 	} xx, yy, zz;
90*25c28e83SPiotr Jasiukajtis 	union {
91*25c28e83SPiotr Jasiukajtis 		unsigned int i[2];
92*25c28e83SPiotr Jasiukajtis 		double d;
93*25c28e83SPiotr Jasiukajtis 	} u;
94*25c28e83SPiotr Jasiukajtis 	double dx[5], dy[5], dxy[9], c, s;
95*25c28e83SPiotr Jasiukajtis 	unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
96*25c28e83SPiotr Jasiukajtis 	unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
97*25c28e83SPiotr Jasiukajtis 	unsigned int rm, sticky;
98*25c28e83SPiotr Jasiukajtis 	unsigned int fsr;
99*25c28e83SPiotr Jasiukajtis 	int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit;
100*25c28e83SPiotr Jasiukajtis 	int cx, cy, cz;
101*25c28e83SPiotr Jasiukajtis 	volatile double	dummy;
102*25c28e83SPiotr Jasiukajtis 
103*25c28e83SPiotr Jasiukajtis 	/* extract the high order words of the arguments */
104*25c28e83SPiotr Jasiukajtis 	xx.q = x;
105*25c28e83SPiotr Jasiukajtis 	yy.q = y;
106*25c28e83SPiotr Jasiukajtis 	zz.q = z;
107*25c28e83SPiotr Jasiukajtis 	hx = xx.i[0] & ~0x80000000;
108*25c28e83SPiotr Jasiukajtis 	hy = yy.i[0] & ~0x80000000;
109*25c28e83SPiotr Jasiukajtis 	hz = zz.i[0] & ~0x80000000;
110*25c28e83SPiotr Jasiukajtis 
111*25c28e83SPiotr Jasiukajtis 	/*
112*25c28e83SPiotr Jasiukajtis 	 * distinguish zero, finite nonzero, infinite, and quiet nan
113*25c28e83SPiotr Jasiukajtis 	 * arguments; raise invalid and return for signaling nans
114*25c28e83SPiotr Jasiukajtis 	 */
115*25c28e83SPiotr Jasiukajtis 	if (hx >= 0x7fff0000) {
116*25c28e83SPiotr Jasiukajtis 		if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) {
117*25c28e83SPiotr Jasiukajtis 			if (!(hx & 0x8000)) {
118*25c28e83SPiotr Jasiukajtis 				/* signaling nan, raise invalid */
119*25c28e83SPiotr Jasiukajtis 				dummy = snan;
120*25c28e83SPiotr Jasiukajtis 				dummy += snan;
121*25c28e83SPiotr Jasiukajtis 				xx.i[0] |= 0x8000;
122*25c28e83SPiotr Jasiukajtis 				return (xx.q);
123*25c28e83SPiotr Jasiukajtis 			}
124*25c28e83SPiotr Jasiukajtis 			cx = 3;	/* quiet nan */
125*25c28e83SPiotr Jasiukajtis 		} else
126*25c28e83SPiotr Jasiukajtis 			cx = 2;	/* inf */
127*25c28e83SPiotr Jasiukajtis 	} else if (hx == 0) {
128*25c28e83SPiotr Jasiukajtis 		cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
129*25c28e83SPiotr Jasiukajtis 				/* subnormal or zero */
130*25c28e83SPiotr Jasiukajtis 	} else
131*25c28e83SPiotr Jasiukajtis 		cx = 1;		/* finite nonzero */
132*25c28e83SPiotr Jasiukajtis 
133*25c28e83SPiotr Jasiukajtis 	if (hy >= 0x7fff0000) {
134*25c28e83SPiotr Jasiukajtis 		if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
135*25c28e83SPiotr Jasiukajtis 			if (!(hy & 0x8000)) {
136*25c28e83SPiotr Jasiukajtis 				dummy = snan;
137*25c28e83SPiotr Jasiukajtis 				dummy += snan;
138*25c28e83SPiotr Jasiukajtis 				yy.i[0] |= 0x8000;
139*25c28e83SPiotr Jasiukajtis 				return (yy.q);
140*25c28e83SPiotr Jasiukajtis 			}
141*25c28e83SPiotr Jasiukajtis 			cy = 3;
142*25c28e83SPiotr Jasiukajtis 		} else
143*25c28e83SPiotr Jasiukajtis 			cy = 2;
144*25c28e83SPiotr Jasiukajtis 	} else if (hy == 0) {
145*25c28e83SPiotr Jasiukajtis 		cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
146*25c28e83SPiotr Jasiukajtis 	} else
147*25c28e83SPiotr Jasiukajtis 		cy = 1;
148*25c28e83SPiotr Jasiukajtis 
149*25c28e83SPiotr Jasiukajtis 	if (hz >= 0x7fff0000) {
150*25c28e83SPiotr Jasiukajtis 		if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
151*25c28e83SPiotr Jasiukajtis 			if (!(hz & 0x8000)) {
152*25c28e83SPiotr Jasiukajtis 				dummy = snan;
153*25c28e83SPiotr Jasiukajtis 				dummy += snan;
154*25c28e83SPiotr Jasiukajtis 				zz.i[0] |= 0x8000;
155*25c28e83SPiotr Jasiukajtis 				return (zz.q);
156*25c28e83SPiotr Jasiukajtis 			}
157*25c28e83SPiotr Jasiukajtis 			cz = 3;
158*25c28e83SPiotr Jasiukajtis 		} else
159*25c28e83SPiotr Jasiukajtis 			cz = 2;
160*25c28e83SPiotr Jasiukajtis 	} else if (hz == 0) {
161*25c28e83SPiotr Jasiukajtis 		cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
162*25c28e83SPiotr Jasiukajtis 	} else
163*25c28e83SPiotr Jasiukajtis 		cz = 1;
164*25c28e83SPiotr Jasiukajtis 
165*25c28e83SPiotr Jasiukajtis 	/* get the fsr and clear current exceptions */
166*25c28e83SPiotr Jasiukajtis 	__fenv_getfsr32(&fsr);
167*25c28e83SPiotr Jasiukajtis 	fsr &= ~FSR_CEXC;
168*25c28e83SPiotr Jasiukajtis 
169*25c28e83SPiotr Jasiukajtis 	/* handle all other zero, inf, and nan cases */
170*25c28e83SPiotr Jasiukajtis 	if (cx != 1 || cy != 1 || cz != 1) {
171*25c28e83SPiotr Jasiukajtis 		/* if x or y is a quiet nan, return it */
172*25c28e83SPiotr Jasiukajtis 		if (cx == 3) {
173*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
174*25c28e83SPiotr Jasiukajtis 			return (x);
175*25c28e83SPiotr Jasiukajtis 		}
176*25c28e83SPiotr Jasiukajtis 		if (cy == 3) {
177*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
178*25c28e83SPiotr Jasiukajtis 			return (y);
179*25c28e83SPiotr Jasiukajtis 		}
180*25c28e83SPiotr Jasiukajtis 
181*25c28e83SPiotr Jasiukajtis 		/* if x*y is 0*inf, raise invalid and return the default nan */
182*25c28e83SPiotr Jasiukajtis 		if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) {
183*25c28e83SPiotr Jasiukajtis 			dummy = zero;
184*25c28e83SPiotr Jasiukajtis 			dummy *= inf;
185*25c28e83SPiotr Jasiukajtis 			zz.i[0] = 0x7fffffff;
186*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
187*25c28e83SPiotr Jasiukajtis 			return (zz.q);
188*25c28e83SPiotr Jasiukajtis 		}
189*25c28e83SPiotr Jasiukajtis 
190*25c28e83SPiotr Jasiukajtis 		/* if z is a quiet nan, return it */
191*25c28e83SPiotr Jasiukajtis 		if (cz == 3) {
192*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
193*25c28e83SPiotr Jasiukajtis 			return (z);
194*25c28e83SPiotr Jasiukajtis 		}
195*25c28e83SPiotr Jasiukajtis 
196*25c28e83SPiotr Jasiukajtis 		/*
197*25c28e83SPiotr Jasiukajtis 		 * now none of x, y, or z is nan; handle cases where x or y
198*25c28e83SPiotr Jasiukajtis 		 * is inf
199*25c28e83SPiotr Jasiukajtis 		 */
200*25c28e83SPiotr Jasiukajtis 		if (cx == 2 || cy == 2) {
201*25c28e83SPiotr Jasiukajtis 			/*
202*25c28e83SPiotr Jasiukajtis 			 * if z is also inf, either we have inf-inf or
203*25c28e83SPiotr Jasiukajtis 			 * the result is the same as z depending on signs
204*25c28e83SPiotr Jasiukajtis 			 */
205*25c28e83SPiotr Jasiukajtis 			if (cz == 2) {
206*25c28e83SPiotr Jasiukajtis 				if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
207*25c28e83SPiotr Jasiukajtis 					dummy = inf;
208*25c28e83SPiotr Jasiukajtis 					dummy -= inf;
209*25c28e83SPiotr Jasiukajtis 					zz.i[0] = 0x7fffffff;
210*25c28e83SPiotr Jasiukajtis 					zz.i[1] = zz.i[2] = zz.i[3] =
211*25c28e83SPiotr Jasiukajtis 						0xffffffff;
212*25c28e83SPiotr Jasiukajtis 					return (zz.q);
213*25c28e83SPiotr Jasiukajtis 				}
214*25c28e83SPiotr Jasiukajtis 				__fenv_setfsr32(&fsr);
215*25c28e83SPiotr Jasiukajtis 				return (z);
216*25c28e83SPiotr Jasiukajtis 			}
217*25c28e83SPiotr Jasiukajtis 
218*25c28e83SPiotr Jasiukajtis 			/* otherwise the result is inf with appropriate sign */
219*25c28e83SPiotr Jasiukajtis 			zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) |
220*25c28e83SPiotr Jasiukajtis 				0x7fff0000;
221*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[2] = zz.i[3] = 0;
222*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
223*25c28e83SPiotr Jasiukajtis 			return (zz.q);
224*25c28e83SPiotr Jasiukajtis 		}
225*25c28e83SPiotr Jasiukajtis 
226*25c28e83SPiotr Jasiukajtis 		/* if z is inf, return it */
227*25c28e83SPiotr Jasiukajtis 		if (cz == 2) {
228*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
229*25c28e83SPiotr Jasiukajtis 			return (z);
230*25c28e83SPiotr Jasiukajtis 		}
231*25c28e83SPiotr Jasiukajtis 
232*25c28e83SPiotr Jasiukajtis 		/*
233*25c28e83SPiotr Jasiukajtis 		 * now x, y, and z are all finite; handle cases where x or y
234*25c28e83SPiotr Jasiukajtis 		 * is zero
235*25c28e83SPiotr Jasiukajtis 		 */
236*25c28e83SPiotr Jasiukajtis 		if (cx == 0 || cy == 0) {
237*25c28e83SPiotr Jasiukajtis 			/* either we have 0-0 or the result is the same as z */
238*25c28e83SPiotr Jasiukajtis 			if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
239*25c28e83SPiotr Jasiukajtis 				0) {
240*25c28e83SPiotr Jasiukajtis 				zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
241*25c28e83SPiotr Jasiukajtis 					0;
242*25c28e83SPiotr Jasiukajtis 				__fenv_setfsr32(&fsr);
243*25c28e83SPiotr Jasiukajtis 				return (zz.q);
244*25c28e83SPiotr Jasiukajtis 			}
245*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
246*25c28e83SPiotr Jasiukajtis 			return (z);
247*25c28e83SPiotr Jasiukajtis 		}
248*25c28e83SPiotr Jasiukajtis 
249*25c28e83SPiotr Jasiukajtis 		/* if we get here, x and y are nonzero finite, z must be zero */
250*25c28e83SPiotr Jasiukajtis 		return (x * y);
251*25c28e83SPiotr Jasiukajtis 	}
252*25c28e83SPiotr Jasiukajtis 
253*25c28e83SPiotr Jasiukajtis 	/*
254*25c28e83SPiotr Jasiukajtis 	 * now x, y, and z are all finite and nonzero; set round-to-
255*25c28e83SPiotr Jasiukajtis 	 * negative-infinity mode
256*25c28e83SPiotr Jasiukajtis 	 */
257*25c28e83SPiotr Jasiukajtis 	__fenv_setfsr32(&fsr_rm);
258*25c28e83SPiotr Jasiukajtis 
259*25c28e83SPiotr Jasiukajtis 	/*
260*25c28e83SPiotr Jasiukajtis 	 * get the signs and exponents and normalize the significands
261*25c28e83SPiotr Jasiukajtis 	 * of x and y
262*25c28e83SPiotr Jasiukajtis 	 */
263*25c28e83SPiotr Jasiukajtis 	sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000;
264*25c28e83SPiotr Jasiukajtis 	ex = hx >> 16;
265*25c28e83SPiotr Jasiukajtis 	hx &= 0xffff;
266*25c28e83SPiotr Jasiukajtis 	if (!ex) {
267*25c28e83SPiotr Jasiukajtis 		if (hx | (xx.i[1] & 0xfffe0000)) {
268*25c28e83SPiotr Jasiukajtis 			ex = 1;
269*25c28e83SPiotr Jasiukajtis 		} else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) {
270*25c28e83SPiotr Jasiukajtis 			hx = xx.i[1];
271*25c28e83SPiotr Jasiukajtis 			xx.i[1] = xx.i[2];
272*25c28e83SPiotr Jasiukajtis 			xx.i[2] = xx.i[3];
273*25c28e83SPiotr Jasiukajtis 			xx.i[3] = 0;
274*25c28e83SPiotr Jasiukajtis 			ex = -31;
275*25c28e83SPiotr Jasiukajtis 		} else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) {
276*25c28e83SPiotr Jasiukajtis 			hx = xx.i[2];
277*25c28e83SPiotr Jasiukajtis 			xx.i[1] = xx.i[3];
278*25c28e83SPiotr Jasiukajtis 			xx.i[2] = xx.i[3] = 0;
279*25c28e83SPiotr Jasiukajtis 			ex = -63;
280*25c28e83SPiotr Jasiukajtis 		} else {
281*25c28e83SPiotr Jasiukajtis 			hx = xx.i[3];
282*25c28e83SPiotr Jasiukajtis 			xx.i[1] = xx.i[2] = xx.i[3] = 0;
283*25c28e83SPiotr Jasiukajtis 			ex = -95;
284*25c28e83SPiotr Jasiukajtis 		}
285*25c28e83SPiotr Jasiukajtis 		while ((hx & 0x10000) == 0) {
286*25c28e83SPiotr Jasiukajtis 			hx = (hx << 1) | (xx.i[1] >> 31);
287*25c28e83SPiotr Jasiukajtis 			xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
288*25c28e83SPiotr Jasiukajtis 			xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
289*25c28e83SPiotr Jasiukajtis 			xx.i[3] <<= 1;
290*25c28e83SPiotr Jasiukajtis 			ex--;
291*25c28e83SPiotr Jasiukajtis 		}
292*25c28e83SPiotr Jasiukajtis 	} else
293*25c28e83SPiotr Jasiukajtis 		hx |= 0x10000;
294*25c28e83SPiotr Jasiukajtis 	ey = hy >> 16;
295*25c28e83SPiotr Jasiukajtis 	hy &= 0xffff;
296*25c28e83SPiotr Jasiukajtis 	if (!ey) {
297*25c28e83SPiotr Jasiukajtis 		if (hy | (yy.i[1] & 0xfffe0000)) {
298*25c28e83SPiotr Jasiukajtis 			ey = 1;
299*25c28e83SPiotr Jasiukajtis 		} else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
300*25c28e83SPiotr Jasiukajtis 			hy = yy.i[1];
301*25c28e83SPiotr Jasiukajtis 			yy.i[1] = yy.i[2];
302*25c28e83SPiotr Jasiukajtis 			yy.i[2] = yy.i[3];
303*25c28e83SPiotr Jasiukajtis 			yy.i[3] = 0;
304*25c28e83SPiotr Jasiukajtis 			ey = -31;
305*25c28e83SPiotr Jasiukajtis 		} else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) {
306*25c28e83SPiotr Jasiukajtis 			hy = yy.i[2];
307*25c28e83SPiotr Jasiukajtis 			yy.i[1] = yy.i[3];
308*25c28e83SPiotr Jasiukajtis 			yy.i[2] = yy.i[3] = 0;
309*25c28e83SPiotr Jasiukajtis 			ey = -63;
310*25c28e83SPiotr Jasiukajtis 		} else {
311*25c28e83SPiotr Jasiukajtis 			hy = yy.i[3];
312*25c28e83SPiotr Jasiukajtis 			yy.i[1] = yy.i[2] = yy.i[3] = 0;
313*25c28e83SPiotr Jasiukajtis 			ey = -95;
314*25c28e83SPiotr Jasiukajtis 		}
315*25c28e83SPiotr Jasiukajtis 		while ((hy & 0x10000) == 0) {
316*25c28e83SPiotr Jasiukajtis 			hy = (hy << 1) | (yy.i[1] >> 31);
317*25c28e83SPiotr Jasiukajtis 			yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
318*25c28e83SPiotr Jasiukajtis 			yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
319*25c28e83SPiotr Jasiukajtis 			yy.i[3] <<= 1;
320*25c28e83SPiotr Jasiukajtis 			ey--;
321*25c28e83SPiotr Jasiukajtis 		}
322*25c28e83SPiotr Jasiukajtis 	} else
323*25c28e83SPiotr Jasiukajtis 		hy |= 0x10000;
324*25c28e83SPiotr Jasiukajtis 	exy = ex + ey - 0x3fff;
325*25c28e83SPiotr Jasiukajtis 
326*25c28e83SPiotr Jasiukajtis 	/* convert the significands of x and y to doubles */
327*25c28e83SPiotr Jasiukajtis 	c = twom16;
328*25c28e83SPiotr Jasiukajtis 	dx[0] = (double) ((int) hx) * c;
329*25c28e83SPiotr Jasiukajtis 	dy[0] = (double) ((int) hy) * c;
330*25c28e83SPiotr Jasiukajtis 
331*25c28e83SPiotr Jasiukajtis 	c *= twom24;
332*25c28e83SPiotr Jasiukajtis 	dx[1] = (double) ((int) (xx.i[1] >> 8)) * c;
333*25c28e83SPiotr Jasiukajtis 	dy[1] = (double) ((int) (yy.i[1] >> 8)) * c;
334*25c28e83SPiotr Jasiukajtis 
335*25c28e83SPiotr Jasiukajtis 	c *= twom24;
336*25c28e83SPiotr Jasiukajtis 	dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) &
337*25c28e83SPiotr Jasiukajtis 	    0xffffff)) * c;
338*25c28e83SPiotr Jasiukajtis 	dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) &
339*25c28e83SPiotr Jasiukajtis 	    0xffffff)) * c;
340*25c28e83SPiotr Jasiukajtis 
341*25c28e83SPiotr Jasiukajtis 	c *= twom24;
342*25c28e83SPiotr Jasiukajtis 	dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) &
343*25c28e83SPiotr Jasiukajtis 	    0xffffff)) * c;
344*25c28e83SPiotr Jasiukajtis 	dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) &
345*25c28e83SPiotr Jasiukajtis 	    0xffffff)) * c;
346*25c28e83SPiotr Jasiukajtis 
347*25c28e83SPiotr Jasiukajtis 	c *= twom24;
348*25c28e83SPiotr Jasiukajtis 	dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c;
349*25c28e83SPiotr Jasiukajtis 	dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c;
350*25c28e83SPiotr Jasiukajtis 
351*25c28e83SPiotr Jasiukajtis 	/* form the "digits" of the product */
352*25c28e83SPiotr Jasiukajtis 	dxy[0] = dx[0] * dy[0];
353*25c28e83SPiotr Jasiukajtis 	dxy[1] = dx[0] * dy[1] + dx[1] * dy[0];
354*25c28e83SPiotr Jasiukajtis 	dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
355*25c28e83SPiotr Jasiukajtis 	dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] +
356*25c28e83SPiotr Jasiukajtis 	    dx[3] * dy[0];
357*25c28e83SPiotr Jasiukajtis 	dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] +
358*25c28e83SPiotr Jasiukajtis 	    dx[3] * dy[1] + dx[4] * dy[0];
359*25c28e83SPiotr Jasiukajtis 	dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] +
360*25c28e83SPiotr Jasiukajtis 	    dx[4] * dy[1];
361*25c28e83SPiotr Jasiukajtis 	dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
362*25c28e83SPiotr Jasiukajtis 	dxy[7] = dx[3] * dy[4] + dx[4] * dy[3];
363*25c28e83SPiotr Jasiukajtis 	dxy[8] = dx[4] * dy[4];
364*25c28e83SPiotr Jasiukajtis 
365*25c28e83SPiotr Jasiukajtis 	/* split odd-numbered terms and combine into even-numbered terms */
366*25c28e83SPiotr Jasiukajtis 	c = (dxy[1] + two20) - two20;
367*25c28e83SPiotr Jasiukajtis 	dxy[0] += c;
368*25c28e83SPiotr Jasiukajtis 	dxy[1] -= c;
369*25c28e83SPiotr Jasiukajtis 	c = (dxy[3] + twom28) - twom28;
370*25c28e83SPiotr Jasiukajtis 	dxy[2] += c + dxy[1];
371*25c28e83SPiotr Jasiukajtis 	dxy[3] -= c;
372*25c28e83SPiotr Jasiukajtis 	c = (dxy[5] + twom76) - twom76;
373*25c28e83SPiotr Jasiukajtis 	dxy[4] += c + dxy[3];
374*25c28e83SPiotr Jasiukajtis 	dxy[5] -= c;
375*25c28e83SPiotr Jasiukajtis 	c = (dxy[7] + twom124) - twom124;
376*25c28e83SPiotr Jasiukajtis 	dxy[6] += c + dxy[5];
377*25c28e83SPiotr Jasiukajtis 	dxy[8] += (dxy[7] - c);
378*25c28e83SPiotr Jasiukajtis 
379*25c28e83SPiotr Jasiukajtis 	/* propagate carries, adjusting the exponent if need be */
380*25c28e83SPiotr Jasiukajtis 	dxy[7] = dxy[6] + dxy[8];
381*25c28e83SPiotr Jasiukajtis 	dxy[5] = dxy[4] + dxy[7];
382*25c28e83SPiotr Jasiukajtis 	dxy[3] = dxy[2] + dxy[5];
383*25c28e83SPiotr Jasiukajtis 	dxy[1] = dxy[0] + dxy[3];
384*25c28e83SPiotr Jasiukajtis 	if (dxy[1] >= two) {
385*25c28e83SPiotr Jasiukajtis 		dxy[0] *= half;
386*25c28e83SPiotr Jasiukajtis 		dxy[1] *= half;
387*25c28e83SPiotr Jasiukajtis 		dxy[2] *= half;
388*25c28e83SPiotr Jasiukajtis 		dxy[3] *= half;
389*25c28e83SPiotr Jasiukajtis 		dxy[4] *= half;
390*25c28e83SPiotr Jasiukajtis 		dxy[5] *= half;
391*25c28e83SPiotr Jasiukajtis 		dxy[6] *= half;
392*25c28e83SPiotr Jasiukajtis 		dxy[7] *= half;
393*25c28e83SPiotr Jasiukajtis 		dxy[8] *= half;
394*25c28e83SPiotr Jasiukajtis 		exy++;
395*25c28e83SPiotr Jasiukajtis 	}
396*25c28e83SPiotr Jasiukajtis 
397*25c28e83SPiotr Jasiukajtis 	/* extract the significand of x*y */
398*25c28e83SPiotr Jasiukajtis 	s = two36;
399*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[1] + s;
400*25c28e83SPiotr Jasiukajtis 	xy0 = u.i[1];
401*25c28e83SPiotr Jasiukajtis 	c -= s;
402*25c28e83SPiotr Jasiukajtis 	dxy[1] -= c;
403*25c28e83SPiotr Jasiukajtis 	dxy[0] -= c;
404*25c28e83SPiotr Jasiukajtis 
405*25c28e83SPiotr Jasiukajtis 	s *= twom32;
406*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[1] + s;
407*25c28e83SPiotr Jasiukajtis 	xy1 = u.i[1];
408*25c28e83SPiotr Jasiukajtis 	c -= s;
409*25c28e83SPiotr Jasiukajtis 	dxy[2] += (dxy[0] - c);
410*25c28e83SPiotr Jasiukajtis 	dxy[3] = dxy[2] + dxy[5];
411*25c28e83SPiotr Jasiukajtis 
412*25c28e83SPiotr Jasiukajtis 	s *= twom32;
413*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[3] + s;
414*25c28e83SPiotr Jasiukajtis 	xy2 = u.i[1];
415*25c28e83SPiotr Jasiukajtis 	c -= s;
416*25c28e83SPiotr Jasiukajtis 	dxy[4] += (dxy[2] - c);
417*25c28e83SPiotr Jasiukajtis 	dxy[5] = dxy[4] + dxy[7];
418*25c28e83SPiotr Jasiukajtis 
419*25c28e83SPiotr Jasiukajtis 	s *= twom32;
420*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[5] + s;
421*25c28e83SPiotr Jasiukajtis 	xy3 = u.i[1];
422*25c28e83SPiotr Jasiukajtis 	c -= s;
423*25c28e83SPiotr Jasiukajtis 	dxy[4] -= c;
424*25c28e83SPiotr Jasiukajtis 	dxy[5] = dxy[4] + dxy[7];
425*25c28e83SPiotr Jasiukajtis 
426*25c28e83SPiotr Jasiukajtis 	s *= twom32;
427*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[5] + s;
428*25c28e83SPiotr Jasiukajtis 	xy4 = u.i[1];
429*25c28e83SPiotr Jasiukajtis 	c -= s;
430*25c28e83SPiotr Jasiukajtis 	dxy[6] += (dxy[4] - c);
431*25c28e83SPiotr Jasiukajtis 	dxy[7] = dxy[6] + dxy[8];
432*25c28e83SPiotr Jasiukajtis 
433*25c28e83SPiotr Jasiukajtis 	s *= twom32;
434*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[7] + s;
435*25c28e83SPiotr Jasiukajtis 	xy5 = u.i[1];
436*25c28e83SPiotr Jasiukajtis 	c -= s;
437*25c28e83SPiotr Jasiukajtis 	dxy[8] += (dxy[6] - c);
438*25c28e83SPiotr Jasiukajtis 
439*25c28e83SPiotr Jasiukajtis 	s *= twom32;
440*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[8] + s;
441*25c28e83SPiotr Jasiukajtis 	xy6 = u.i[1];
442*25c28e83SPiotr Jasiukajtis 	c -= s;
443*25c28e83SPiotr Jasiukajtis 	dxy[8] -= c;
444*25c28e83SPiotr Jasiukajtis 
445*25c28e83SPiotr Jasiukajtis 	s *= twom32;
446*25c28e83SPiotr Jasiukajtis 	u.d = c = dxy[8] + s;
447*25c28e83SPiotr Jasiukajtis 	xy7 = u.i[1];
448*25c28e83SPiotr Jasiukajtis 
449*25c28e83SPiotr Jasiukajtis 	/* extract the sign, exponent, and significand of z */
450*25c28e83SPiotr Jasiukajtis 	sz = zz.i[0] & 0x80000000;
451*25c28e83SPiotr Jasiukajtis 	ez = hz >> 16;
452*25c28e83SPiotr Jasiukajtis 	z0 = hz & 0xffff;
453*25c28e83SPiotr Jasiukajtis 	if (!ez) {
454*25c28e83SPiotr Jasiukajtis 		if (z0 | (zz.i[1] & 0xfffe0000)) {
455*25c28e83SPiotr Jasiukajtis 			z1 = zz.i[1];
456*25c28e83SPiotr Jasiukajtis 			z2 = zz.i[2];
457*25c28e83SPiotr Jasiukajtis 			z3 = zz.i[3];
458*25c28e83SPiotr Jasiukajtis 			ez = 1;
459*25c28e83SPiotr Jasiukajtis 		} else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) {
460*25c28e83SPiotr Jasiukajtis 			z0 = zz.i[1];
461*25c28e83SPiotr Jasiukajtis 			z1 = zz.i[2];
462*25c28e83SPiotr Jasiukajtis 			z2 = zz.i[3];
463*25c28e83SPiotr Jasiukajtis 			z3 = 0;
464*25c28e83SPiotr Jasiukajtis 			ez = -31;
465*25c28e83SPiotr Jasiukajtis 		} else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) {
466*25c28e83SPiotr Jasiukajtis 			z0 = zz.i[2];
467*25c28e83SPiotr Jasiukajtis 			z1 = zz.i[3];
468*25c28e83SPiotr Jasiukajtis 			z2 = z3 = 0;
469*25c28e83SPiotr Jasiukajtis 			ez = -63;
470*25c28e83SPiotr Jasiukajtis 		} else {
471*25c28e83SPiotr Jasiukajtis 			z0 = zz.i[3];
472*25c28e83SPiotr Jasiukajtis 			z1 = z2 = z3 = 0;
473*25c28e83SPiotr Jasiukajtis 			ez = -95;
474*25c28e83SPiotr Jasiukajtis 		}
475*25c28e83SPiotr Jasiukajtis 		while ((z0 & 0x10000) == 0) {
476*25c28e83SPiotr Jasiukajtis 			z0 = (z0 << 1) | (z1 >> 31);
477*25c28e83SPiotr Jasiukajtis 			z1 = (z1 << 1) | (z2 >> 31);
478*25c28e83SPiotr Jasiukajtis 			z2 = (z2 << 1) | (z3 >> 31);
479*25c28e83SPiotr Jasiukajtis 			z3 <<= 1;
480*25c28e83SPiotr Jasiukajtis 			ez--;
481*25c28e83SPiotr Jasiukajtis 		}
482*25c28e83SPiotr Jasiukajtis 	} else {
483*25c28e83SPiotr Jasiukajtis 		z0 |= 0x10000;
484*25c28e83SPiotr Jasiukajtis 		z1 = zz.i[1];
485*25c28e83SPiotr Jasiukajtis 		z2 = zz.i[2];
486*25c28e83SPiotr Jasiukajtis 		z3 = zz.i[3];
487*25c28e83SPiotr Jasiukajtis 	}
488*25c28e83SPiotr Jasiukajtis 	z4 = z5 = z6 = z7 = 0;
489*25c28e83SPiotr Jasiukajtis 
490*25c28e83SPiotr Jasiukajtis 	/*
491*25c28e83SPiotr Jasiukajtis 	 * now x*y is represented by sxy, exy, and xy[0-7], and z is
492*25c28e83SPiotr Jasiukajtis 	 * represented likewise; swap if need be so |xy| <= |z|
493*25c28e83SPiotr Jasiukajtis 	 */
494*25c28e83SPiotr Jasiukajtis 	if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
495*25c28e83SPiotr Jasiukajtis 		(xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 ||
496*25c28e83SPiotr Jasiukajtis 		(xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
497*25c28e83SPiotr Jasiukajtis 		e = sxy; sxy = sz; sz = e;
498*25c28e83SPiotr Jasiukajtis 		e = exy; exy = ez; ez = e;
499*25c28e83SPiotr Jasiukajtis 		e = xy0; xy0 = z0; z0 = e;
500*25c28e83SPiotr Jasiukajtis 		e = xy1; xy1 = z1; z1 = e;
501*25c28e83SPiotr Jasiukajtis 		e = xy2; xy2 = z2; z2 = e;
502*25c28e83SPiotr Jasiukajtis 		e = xy3; xy3 = z3; z3 = e;
503*25c28e83SPiotr Jasiukajtis 		z4 = xy4; xy4 = 0;
504*25c28e83SPiotr Jasiukajtis 		z5 = xy5; xy5 = 0;
505*25c28e83SPiotr Jasiukajtis 		z6 = xy6; xy6 = 0;
506*25c28e83SPiotr Jasiukajtis 		z7 = xy7; xy7 = 0;
507*25c28e83SPiotr Jasiukajtis 	}
508*25c28e83SPiotr Jasiukajtis 
509*25c28e83SPiotr Jasiukajtis 	/* shift the significand of xy keeping a sticky bit */
510*25c28e83SPiotr Jasiukajtis 	e = ez - exy;
511*25c28e83SPiotr Jasiukajtis 	if (e > 236) {
512*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
513*25c28e83SPiotr Jasiukajtis 		xy7 = 1;
514*25c28e83SPiotr Jasiukajtis 	} else if (e >= 224) {
515*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 |
516*25c28e83SPiotr Jasiukajtis 			((xy0 << 1) << (255 - e));
517*25c28e83SPiotr Jasiukajtis 		xy7 = xy0 >> (e - 224);
518*25c28e83SPiotr Jasiukajtis 		if (sticky)
519*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
520*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
521*25c28e83SPiotr Jasiukajtis 	} else if (e >= 192) {
522*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 |
523*25c28e83SPiotr Jasiukajtis 			((xy1 << 1) << (223 - e));
524*25c28e83SPiotr Jasiukajtis 		xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
525*25c28e83SPiotr Jasiukajtis 		if (sticky)
526*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
527*25c28e83SPiotr Jasiukajtis 		xy6 = xy0 >> (e - 192);
528*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0;
529*25c28e83SPiotr Jasiukajtis 	} else if (e >= 160) {
530*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | xy5 | xy4 | xy3 |
531*25c28e83SPiotr Jasiukajtis 			((xy2 << 1) << (191 - e));
532*25c28e83SPiotr Jasiukajtis 		xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
533*25c28e83SPiotr Jasiukajtis 		if (sticky)
534*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
535*25c28e83SPiotr Jasiukajtis 		xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
536*25c28e83SPiotr Jasiukajtis 		xy5 = xy0 >> (e - 160);
537*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = xy4 = 0;
538*25c28e83SPiotr Jasiukajtis 	} else if (e >= 128) {
539*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
540*25c28e83SPiotr Jasiukajtis 		xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
541*25c28e83SPiotr Jasiukajtis 		if (sticky)
542*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
543*25c28e83SPiotr Jasiukajtis 		xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
544*25c28e83SPiotr Jasiukajtis 		xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
545*25c28e83SPiotr Jasiukajtis 		xy4 = xy0 >> (e - 128);
546*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = 0;
547*25c28e83SPiotr Jasiukajtis 	} else if (e >= 96) {
548*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
549*25c28e83SPiotr Jasiukajtis 		xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
550*25c28e83SPiotr Jasiukajtis 		if (sticky)
551*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
552*25c28e83SPiotr Jasiukajtis 		xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
553*25c28e83SPiotr Jasiukajtis 		xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
554*25c28e83SPiotr Jasiukajtis 		xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
555*25c28e83SPiotr Jasiukajtis 		xy3 = xy0 >> (e - 96);
556*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = 0;
557*25c28e83SPiotr Jasiukajtis 	} else if (e >= 64) {
558*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
559*25c28e83SPiotr Jasiukajtis 		xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
560*25c28e83SPiotr Jasiukajtis 		if (sticky)
561*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
562*25c28e83SPiotr Jasiukajtis 		xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
563*25c28e83SPiotr Jasiukajtis 		xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
564*25c28e83SPiotr Jasiukajtis 		xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
565*25c28e83SPiotr Jasiukajtis 		xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
566*25c28e83SPiotr Jasiukajtis 		xy2 = xy0 >> (e - 64);
567*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = 0;
568*25c28e83SPiotr Jasiukajtis 	} else if (e >= 32) {
569*25c28e83SPiotr Jasiukajtis 		sticky = xy7 | ((xy6 << 1) << (63 - e));
570*25c28e83SPiotr Jasiukajtis 		xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
571*25c28e83SPiotr Jasiukajtis 		if (sticky)
572*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
573*25c28e83SPiotr Jasiukajtis 		xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
574*25c28e83SPiotr Jasiukajtis 		xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
575*25c28e83SPiotr Jasiukajtis 		xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
576*25c28e83SPiotr Jasiukajtis 		xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
577*25c28e83SPiotr Jasiukajtis 		xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
578*25c28e83SPiotr Jasiukajtis 		xy1 = xy0 >> (e - 32);
579*25c28e83SPiotr Jasiukajtis 		xy0 = 0;
580*25c28e83SPiotr Jasiukajtis 	} else if (e) {
581*25c28e83SPiotr Jasiukajtis 		sticky = (xy7 << 1) << (31 - e);
582*25c28e83SPiotr Jasiukajtis 		xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
583*25c28e83SPiotr Jasiukajtis 		if (sticky)
584*25c28e83SPiotr Jasiukajtis 			xy7 |= 1;
585*25c28e83SPiotr Jasiukajtis 		xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
586*25c28e83SPiotr Jasiukajtis 		xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
587*25c28e83SPiotr Jasiukajtis 		xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
588*25c28e83SPiotr Jasiukajtis 		xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
589*25c28e83SPiotr Jasiukajtis 		xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
590*25c28e83SPiotr Jasiukajtis 		xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
591*25c28e83SPiotr Jasiukajtis 		xy0 >>= e;
592*25c28e83SPiotr Jasiukajtis 	}
593*25c28e83SPiotr Jasiukajtis 
594*25c28e83SPiotr Jasiukajtis 	/* if this is a magnitude subtract, negate the significand of xy */
595*25c28e83SPiotr Jasiukajtis 	if (sxy ^ sz) {
596*25c28e83SPiotr Jasiukajtis 		xy0 = ~xy0;
597*25c28e83SPiotr Jasiukajtis 		xy1 = ~xy1;
598*25c28e83SPiotr Jasiukajtis 		xy2 = ~xy2;
599*25c28e83SPiotr Jasiukajtis 		xy3 = ~xy3;
600*25c28e83SPiotr Jasiukajtis 		xy4 = ~xy4;
601*25c28e83SPiotr Jasiukajtis 		xy5 = ~xy5;
602*25c28e83SPiotr Jasiukajtis 		xy6 = ~xy6;
603*25c28e83SPiotr Jasiukajtis 		xy7 = -xy7;
604*25c28e83SPiotr Jasiukajtis 		if (xy7 == 0)
605*25c28e83SPiotr Jasiukajtis 			if (++xy6 == 0)
606*25c28e83SPiotr Jasiukajtis 				if (++xy5 == 0)
607*25c28e83SPiotr Jasiukajtis 					if (++xy4 == 0)
608*25c28e83SPiotr Jasiukajtis 						if (++xy3 == 0)
609*25c28e83SPiotr Jasiukajtis 							if (++xy2 == 0)
610*25c28e83SPiotr Jasiukajtis 								if (++xy1 == 0)
611*25c28e83SPiotr Jasiukajtis 									xy0++;
612*25c28e83SPiotr Jasiukajtis 	}
613*25c28e83SPiotr Jasiukajtis 
614*25c28e83SPiotr Jasiukajtis 	/* add, propagating carries */
615*25c28e83SPiotr Jasiukajtis 	z7 += xy7;
616*25c28e83SPiotr Jasiukajtis 	e = (z7 < xy7);
617*25c28e83SPiotr Jasiukajtis 	z6 += xy6;
618*25c28e83SPiotr Jasiukajtis 	if (e) {
619*25c28e83SPiotr Jasiukajtis 		z6++;
620*25c28e83SPiotr Jasiukajtis 		e = (z6 <= xy6);
621*25c28e83SPiotr Jasiukajtis 	} else
622*25c28e83SPiotr Jasiukajtis 		e = (z6 < xy6);
623*25c28e83SPiotr Jasiukajtis 	z5 += xy5;
624*25c28e83SPiotr Jasiukajtis 	if (e) {
625*25c28e83SPiotr Jasiukajtis 		z5++;
626*25c28e83SPiotr Jasiukajtis 		e = (z5 <= xy5);
627*25c28e83SPiotr Jasiukajtis 	} else
628*25c28e83SPiotr Jasiukajtis 		e = (z5 < xy5);
629*25c28e83SPiotr Jasiukajtis 	z4 += xy4;
630*25c28e83SPiotr Jasiukajtis 	if (e) {
631*25c28e83SPiotr Jasiukajtis 		z4++;
632*25c28e83SPiotr Jasiukajtis 		e = (z4 <= xy4);
633*25c28e83SPiotr Jasiukajtis 	} else
634*25c28e83SPiotr Jasiukajtis 		e = (z4 < xy4);
635*25c28e83SPiotr Jasiukajtis 	z3 += xy3;
636*25c28e83SPiotr Jasiukajtis 	if (e) {
637*25c28e83SPiotr Jasiukajtis 		z3++;
638*25c28e83SPiotr Jasiukajtis 		e = (z3 <= xy3);
639*25c28e83SPiotr Jasiukajtis 	} else
640*25c28e83SPiotr Jasiukajtis 		e = (z3 < xy3);
641*25c28e83SPiotr Jasiukajtis 	z2 += xy2;
642*25c28e83SPiotr Jasiukajtis 	if (e) {
643*25c28e83SPiotr Jasiukajtis 		z2++;
644*25c28e83SPiotr Jasiukajtis 		e = (z2 <= xy2);
645*25c28e83SPiotr Jasiukajtis 	} else
646*25c28e83SPiotr Jasiukajtis 		e = (z2 < xy2);
647*25c28e83SPiotr Jasiukajtis 	z1 += xy1;
648*25c28e83SPiotr Jasiukajtis 	if (e) {
649*25c28e83SPiotr Jasiukajtis 		z1++;
650*25c28e83SPiotr Jasiukajtis 		e = (z1 <= xy1);
651*25c28e83SPiotr Jasiukajtis 	} else
652*25c28e83SPiotr Jasiukajtis 		e = (z1 < xy1);
653*25c28e83SPiotr Jasiukajtis 	z0 += xy0;
654*25c28e83SPiotr Jasiukajtis 	if (e)
655*25c28e83SPiotr Jasiukajtis 		z0++;
656*25c28e83SPiotr Jasiukajtis 
657*25c28e83SPiotr Jasiukajtis 	/* postnormalize and collect rounding information into z4 */
658*25c28e83SPiotr Jasiukajtis 	if (ez < 1) {
659*25c28e83SPiotr Jasiukajtis 		/* result is tiny; shift right until exponent is within range */
660*25c28e83SPiotr Jasiukajtis 		e = 1 - ez;
661*25c28e83SPiotr Jasiukajtis 		if (e > 116) {
662*25c28e83SPiotr Jasiukajtis 			z4 = 1; /* result can't be exactly zero */
663*25c28e83SPiotr Jasiukajtis 			z0 = z1 = z2 = z3 = 0;
664*25c28e83SPiotr Jasiukajtis 		} else if (e >= 96) {
665*25c28e83SPiotr Jasiukajtis 			sticky = z7 | z6 | z5 | z4 | z3 | z2 |
666*25c28e83SPiotr Jasiukajtis 				((z1 << 1) << (127 - e));
667*25c28e83SPiotr Jasiukajtis 			z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
668*25c28e83SPiotr Jasiukajtis 			if (sticky)
669*25c28e83SPiotr Jasiukajtis 				z4 |= 1;
670*25c28e83SPiotr Jasiukajtis 			z3 = z0 >> (e - 96);
671*25c28e83SPiotr Jasiukajtis 			z0 = z1 = z2 = 0;
672*25c28e83SPiotr Jasiukajtis 		} else if (e >= 64) {
673*25c28e83SPiotr Jasiukajtis 			sticky = z7 | z6 | z5 | z4 | z3 |
674*25c28e83SPiotr Jasiukajtis 				((z2 << 1) << (95 - e));
675*25c28e83SPiotr Jasiukajtis 			z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
676*25c28e83SPiotr Jasiukajtis 			if (sticky)
677*25c28e83SPiotr Jasiukajtis 				z4 |= 1;
678*25c28e83SPiotr Jasiukajtis 			z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
679*25c28e83SPiotr Jasiukajtis 			z2 = z0 >> (e - 64);
680*25c28e83SPiotr Jasiukajtis 			z0 = z1 = 0;
681*25c28e83SPiotr Jasiukajtis 		} else if (e >= 32) {
682*25c28e83SPiotr Jasiukajtis 			sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
683*25c28e83SPiotr Jasiukajtis 			z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
684*25c28e83SPiotr Jasiukajtis 			if (sticky)
685*25c28e83SPiotr Jasiukajtis 				z4 |= 1;
686*25c28e83SPiotr Jasiukajtis 			z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
687*25c28e83SPiotr Jasiukajtis 			z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
688*25c28e83SPiotr Jasiukajtis 			z1 = z0 >> (e - 32);
689*25c28e83SPiotr Jasiukajtis 			z0 = 0;
690*25c28e83SPiotr Jasiukajtis 		} else {
691*25c28e83SPiotr Jasiukajtis 			sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
692*25c28e83SPiotr Jasiukajtis 			z4 = (z4 >> e) | ((z3 << 1) << (31 - e));
693*25c28e83SPiotr Jasiukajtis 			if (sticky)
694*25c28e83SPiotr Jasiukajtis 				z4 |= 1;
695*25c28e83SPiotr Jasiukajtis 			z3 = (z3 >> e) | ((z2 << 1) << (31 - e));
696*25c28e83SPiotr Jasiukajtis 			z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
697*25c28e83SPiotr Jasiukajtis 			z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
698*25c28e83SPiotr Jasiukajtis 			z0 >>= e;
699*25c28e83SPiotr Jasiukajtis 		}
700*25c28e83SPiotr Jasiukajtis 		ez = 1;
701*25c28e83SPiotr Jasiukajtis 	} else if (z0 >= 0x20000) {
702*25c28e83SPiotr Jasiukajtis 		/* carry out; shift right by one */
703*25c28e83SPiotr Jasiukajtis 		sticky = (z4 & 1) | z5 | z6 | z7;
704*25c28e83SPiotr Jasiukajtis 		z4 = (z4 >> 1) | (z3 << 31);
705*25c28e83SPiotr Jasiukajtis 		if (sticky)
706*25c28e83SPiotr Jasiukajtis 			z4 |= 1;
707*25c28e83SPiotr Jasiukajtis 		z3 = (z3 >> 1) | (z2 << 31);
708*25c28e83SPiotr Jasiukajtis 		z2 = (z2 >> 1) | (z1 << 31);
709*25c28e83SPiotr Jasiukajtis 		z1 = (z1 >> 1) | (z0 << 31);
710*25c28e83SPiotr Jasiukajtis 		z0 >>= 1;
711*25c28e83SPiotr Jasiukajtis 		ez++;
712*25c28e83SPiotr Jasiukajtis 	} else {
713*25c28e83SPiotr Jasiukajtis 		if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7)
714*25c28e83SPiotr Jasiukajtis 			!= 0) {
715*25c28e83SPiotr Jasiukajtis 			/*
716*25c28e83SPiotr Jasiukajtis 			 * borrow/cancellation; shift left as much as
717*25c28e83SPiotr Jasiukajtis 			 * exponent allows
718*25c28e83SPiotr Jasiukajtis 			 */
719*25c28e83SPiotr Jasiukajtis 			while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
720*25c28e83SPiotr Jasiukajtis 				z0 = z1;
721*25c28e83SPiotr Jasiukajtis 				z1 = z2;
722*25c28e83SPiotr Jasiukajtis 				z2 = z3;
723*25c28e83SPiotr Jasiukajtis 				z3 = z4;
724*25c28e83SPiotr Jasiukajtis 				z4 = z5;
725*25c28e83SPiotr Jasiukajtis 				z5 = z6;
726*25c28e83SPiotr Jasiukajtis 				z6 = z7;
727*25c28e83SPiotr Jasiukajtis 				z7 = 0;
728*25c28e83SPiotr Jasiukajtis 				ez -= 32;
729*25c28e83SPiotr Jasiukajtis 			}
730*25c28e83SPiotr Jasiukajtis 			while (z0 < 0x10000 && ez > 1) {
731*25c28e83SPiotr Jasiukajtis 				z0 = (z0 << 1) | (z1 >> 31);
732*25c28e83SPiotr Jasiukajtis 				z1 = (z1 << 1) | (z2 >> 31);
733*25c28e83SPiotr Jasiukajtis 				z2 = (z2 << 1) | (z3 >> 31);
734*25c28e83SPiotr Jasiukajtis 				z3 = (z3 << 1) | (z4 >> 31);
735*25c28e83SPiotr Jasiukajtis 				z4 = (z4 << 1) | (z5 >> 31);
736*25c28e83SPiotr Jasiukajtis 				z5 = (z5 << 1) | (z6 >> 31);
737*25c28e83SPiotr Jasiukajtis 				z6 = (z6 << 1) | (z7 >> 31);
738*25c28e83SPiotr Jasiukajtis 				z7 <<= 1;
739*25c28e83SPiotr Jasiukajtis 				ez--;
740*25c28e83SPiotr Jasiukajtis 			}
741*25c28e83SPiotr Jasiukajtis 		}
742*25c28e83SPiotr Jasiukajtis 		if (z5 | z6 | z7)
743*25c28e83SPiotr Jasiukajtis 			z4 |= 1;
744*25c28e83SPiotr Jasiukajtis 	}
745*25c28e83SPiotr Jasiukajtis 
746*25c28e83SPiotr Jasiukajtis 	/* get the rounding mode */
747*25c28e83SPiotr Jasiukajtis 	rm = fsr >> 30;
748*25c28e83SPiotr Jasiukajtis 
749*25c28e83SPiotr Jasiukajtis 	/* strip off the integer bit, if there is one */
750*25c28e83SPiotr Jasiukajtis 	ibit = z0 & 0x10000;
751*25c28e83SPiotr Jasiukajtis 	if (ibit)
752*25c28e83SPiotr Jasiukajtis 		z0 -= 0x10000;
753*25c28e83SPiotr Jasiukajtis 	else {
754*25c28e83SPiotr Jasiukajtis 		ez = 0;
755*25c28e83SPiotr Jasiukajtis 		if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
756*25c28e83SPiotr Jasiukajtis 			zz.i[0] = rm == FSR_RM ? 0x80000000 : 0;
757*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[2] = zz.i[3] = 0;
758*25c28e83SPiotr Jasiukajtis 			__fenv_setfsr32(&fsr);
759*25c28e83SPiotr Jasiukajtis 			return (zz.q);
760*25c28e83SPiotr Jasiukajtis 		}
761*25c28e83SPiotr Jasiukajtis 	}
762*25c28e83SPiotr Jasiukajtis 
763*25c28e83SPiotr Jasiukajtis 	/*
764*25c28e83SPiotr Jasiukajtis 	 * flip the sense of directed roundings if the result is negative;
765*25c28e83SPiotr Jasiukajtis 	 * the logic below applies to a positive result
766*25c28e83SPiotr Jasiukajtis 	 */
767*25c28e83SPiotr Jasiukajtis 	if (sz)
768*25c28e83SPiotr Jasiukajtis 		rm ^= rm >> 1;
769*25c28e83SPiotr Jasiukajtis 
770*25c28e83SPiotr Jasiukajtis 	/* round and raise exceptions */
771*25c28e83SPiotr Jasiukajtis 	if (z4) {
772*25c28e83SPiotr Jasiukajtis 		fsr |= FSR_NXC;
773*25c28e83SPiotr Jasiukajtis 
774*25c28e83SPiotr Jasiukajtis 		/* decide whether to round the fraction up */
775*25c28e83SPiotr Jasiukajtis 		if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u ||
776*25c28e83SPiotr Jasiukajtis 			(z4 == 0x80000000u && (z3 & 1))))) {
777*25c28e83SPiotr Jasiukajtis 			/* round up and renormalize if necessary */
778*25c28e83SPiotr Jasiukajtis 			if (++z3 == 0)
779*25c28e83SPiotr Jasiukajtis 				if (++z2 == 0)
780*25c28e83SPiotr Jasiukajtis 					if (++z1 == 0)
781*25c28e83SPiotr Jasiukajtis 						if (++z0 == 0x10000) {
782*25c28e83SPiotr Jasiukajtis 							z0 = 0;
783*25c28e83SPiotr Jasiukajtis 							ez++;
784*25c28e83SPiotr Jasiukajtis 						}
785*25c28e83SPiotr Jasiukajtis 		}
786*25c28e83SPiotr Jasiukajtis 	}
787*25c28e83SPiotr Jasiukajtis 
788*25c28e83SPiotr Jasiukajtis 	/* check for under/overflow */
789*25c28e83SPiotr Jasiukajtis 	if (ez >= 0x7fff) {
790*25c28e83SPiotr Jasiukajtis 		if (rm == FSR_RN || rm == FSR_RP) {
791*25c28e83SPiotr Jasiukajtis 			zz.i[0] = sz | 0x7fff0000;
792*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[2] = zz.i[3] = 0;
793*25c28e83SPiotr Jasiukajtis 		} else {
794*25c28e83SPiotr Jasiukajtis 			zz.i[0] = sz | 0x7ffeffff;
795*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
796*25c28e83SPiotr Jasiukajtis 		}
797*25c28e83SPiotr Jasiukajtis 		fsr |= FSR_OFC | FSR_NXC;
798*25c28e83SPiotr Jasiukajtis 	} else {
799*25c28e83SPiotr Jasiukajtis 		zz.i[0] = sz | (ez << 16) | z0;
800*25c28e83SPiotr Jasiukajtis 		zz.i[1] = z1;
801*25c28e83SPiotr Jasiukajtis 		zz.i[2] = z2;
802*25c28e83SPiotr Jasiukajtis 		zz.i[3] = z3;
803*25c28e83SPiotr Jasiukajtis 
804*25c28e83SPiotr Jasiukajtis 		/*
805*25c28e83SPiotr Jasiukajtis 		 * !ibit => exact result was tiny before rounding,
806*25c28e83SPiotr Jasiukajtis 		 * z4 nonzero => result delivered is inexact
807*25c28e83SPiotr Jasiukajtis 		 */
808*25c28e83SPiotr Jasiukajtis 		if (!ibit) {
809*25c28e83SPiotr Jasiukajtis 			if (z4)
810*25c28e83SPiotr Jasiukajtis 				fsr |= FSR_UFC | FSR_NXC;
811*25c28e83SPiotr Jasiukajtis 			else if (fsr & FSR_UFM)
812*25c28e83SPiotr Jasiukajtis 				fsr |= FSR_UFC;
813*25c28e83SPiotr Jasiukajtis 		}
814*25c28e83SPiotr Jasiukajtis 	}
815*25c28e83SPiotr Jasiukajtis 
816*25c28e83SPiotr Jasiukajtis 	/* restore the fsr and emulate exceptions as needed */
817*25c28e83SPiotr Jasiukajtis 	if ((fsr & FSR_CEXC) & (fsr >> 23)) {
818*25c28e83SPiotr Jasiukajtis 		__fenv_setfsr32(&fsr);
819*25c28e83SPiotr Jasiukajtis 		if (fsr & FSR_OFC) {
820*25c28e83SPiotr Jasiukajtis 			dummy = huge;
821*25c28e83SPiotr Jasiukajtis 			dummy *= huge;
822*25c28e83SPiotr Jasiukajtis 		} else if (fsr & FSR_UFC) {
823*25c28e83SPiotr Jasiukajtis 			dummy = tiny;
824*25c28e83SPiotr Jasiukajtis 			if (fsr & FSR_NXC)
825*25c28e83SPiotr Jasiukajtis 				dummy *= tiny;
826*25c28e83SPiotr Jasiukajtis 			else
827*25c28e83SPiotr Jasiukajtis 				dummy -= tiny2;
828*25c28e83SPiotr Jasiukajtis 		} else {
829*25c28e83SPiotr Jasiukajtis 			dummy = huge;
830*25c28e83SPiotr Jasiukajtis 			dummy += tiny;
831*25c28e83SPiotr Jasiukajtis 		}
832*25c28e83SPiotr Jasiukajtis 	} else {
833*25c28e83SPiotr Jasiukajtis 		fsr |= (fsr & 0x1f) << 5;
834*25c28e83SPiotr Jasiukajtis 		__fenv_setfsr32(&fsr);
835*25c28e83SPiotr Jasiukajtis 	}
836*25c28e83SPiotr Jasiukajtis 	return (zz.q);
837*25c28e83SPiotr Jasiukajtis }
838*25c28e83SPiotr Jasiukajtis 
839*25c28e83SPiotr Jasiukajtis #elif defined(__x86)
840*25c28e83SPiotr Jasiukajtis 
841*25c28e83SPiotr Jasiukajtis static const union {
842*25c28e83SPiotr Jasiukajtis 	unsigned i[2];
843*25c28e83SPiotr Jasiukajtis 	double d;
844*25c28e83SPiotr Jasiukajtis } C[] = {
845*25c28e83SPiotr Jasiukajtis 	{ 0, 0x3fe00000u },
846*25c28e83SPiotr Jasiukajtis 	{ 0, 0x40000000u },
847*25c28e83SPiotr Jasiukajtis 	{ 0, 0x3df00000u },
848*25c28e83SPiotr Jasiukajtis 	{ 0, 0x3bf00000u },
849*25c28e83SPiotr Jasiukajtis 	{ 0, 0x41f00000u },
850*25c28e83SPiotr Jasiukajtis 	{ 0, 0x43e00000u },
851*25c28e83SPiotr Jasiukajtis 	{ 0, 0x7fe00000u },
852*25c28e83SPiotr Jasiukajtis 	{ 0, 0x00100000u },
853*25c28e83SPiotr Jasiukajtis 	{ 0, 0x00100001u }
854*25c28e83SPiotr Jasiukajtis };
855*25c28e83SPiotr Jasiukajtis 
856*25c28e83SPiotr Jasiukajtis #define	half	C[0].d
857*25c28e83SPiotr Jasiukajtis #define	two	C[1].d
858*25c28e83SPiotr Jasiukajtis #define	twom32	C[2].d
859*25c28e83SPiotr Jasiukajtis #define	twom64	C[3].d
860*25c28e83SPiotr Jasiukajtis #define	two32	C[4].d
861*25c28e83SPiotr Jasiukajtis #define	two63	C[5].d
862*25c28e83SPiotr Jasiukajtis #define	huge	C[6].d
863*25c28e83SPiotr Jasiukajtis #define	tiny	C[7].d
864*25c28e83SPiotr Jasiukajtis #define	tiny2	C[8].d
865*25c28e83SPiotr Jasiukajtis 
866*25c28e83SPiotr Jasiukajtis #if defined(__amd64)
867*25c28e83SPiotr Jasiukajtis #define	NI	4
868*25c28e83SPiotr Jasiukajtis #else
869*25c28e83SPiotr Jasiukajtis #define	NI	3
870*25c28e83SPiotr Jasiukajtis #endif
871*25c28e83SPiotr Jasiukajtis 
872*25c28e83SPiotr Jasiukajtis /*
873*25c28e83SPiotr Jasiukajtis  * fmal for x86: 80-bit extended double precision, little-endian
874*25c28e83SPiotr Jasiukajtis  */
875*25c28e83SPiotr Jasiukajtis long double
876*25c28e83SPiotr Jasiukajtis __fmal(long double x, long double y, long double z) {
877*25c28e83SPiotr Jasiukajtis 	union {
878*25c28e83SPiotr Jasiukajtis 		unsigned i[NI];
879*25c28e83SPiotr Jasiukajtis 		long double e;
880*25c28e83SPiotr Jasiukajtis 	} xx, yy, zz;
881*25c28e83SPiotr Jasiukajtis 	long double xhi, yhi, xlo, ylo, t;
882*25c28e83SPiotr Jasiukajtis 	unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
883*25c28e83SPiotr Jasiukajtis 	unsigned oldcwsw, cwsw, rm, sticky, carry;
884*25c28e83SPiotr Jasiukajtis 	int ex, ey, ez, exy, sxy, sz, e, tinyafter;
885*25c28e83SPiotr Jasiukajtis 	volatile double	dummy;
886*25c28e83SPiotr Jasiukajtis 
887*25c28e83SPiotr Jasiukajtis 	/* extract the exponents of the arguments */
888*25c28e83SPiotr Jasiukajtis 	xx.e = x;
889*25c28e83SPiotr Jasiukajtis 	yy.e = y;
890*25c28e83SPiotr Jasiukajtis 	zz.e = z;
891*25c28e83SPiotr Jasiukajtis 	ex = xx.i[2] & 0x7fff;
892*25c28e83SPiotr Jasiukajtis 	ey = yy.i[2] & 0x7fff;
893*25c28e83SPiotr Jasiukajtis 	ez = zz.i[2] & 0x7fff;
894*25c28e83SPiotr Jasiukajtis 
895*25c28e83SPiotr Jasiukajtis 	/* dispense with inf, nan, and zero cases */
896*25c28e83SPiotr Jasiukajtis 	if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 ||
897*25c28e83SPiotr Jasiukajtis 		(ey | yy.i[1] | yy.i[0]) == 0)	/* x or y is inf, nan, or 0 */
898*25c28e83SPiotr Jasiukajtis 		return (x * y + z);
899*25c28e83SPiotr Jasiukajtis 
900*25c28e83SPiotr Jasiukajtis 	if (ez == 0x7fff)			/* z is inf or nan */
901*25c28e83SPiotr Jasiukajtis 		return (x + z);	/* avoid spurious under/overflow in x * y */
902*25c28e83SPiotr Jasiukajtis 
903*25c28e83SPiotr Jasiukajtis 	if ((ez | zz.i[1] | zz.i[0]) == 0)	/* z is zero */
904*25c28e83SPiotr Jasiukajtis 		/*
905*25c28e83SPiotr Jasiukajtis 		 * x * y isn't zero but could underflow to zero,
906*25c28e83SPiotr Jasiukajtis 		 * so don't add z, lest we perturb the sign
907*25c28e83SPiotr Jasiukajtis 		 */
908*25c28e83SPiotr Jasiukajtis 		return (x * y);
909*25c28e83SPiotr Jasiukajtis 
910*25c28e83SPiotr Jasiukajtis 	/*
911*25c28e83SPiotr Jasiukajtis 	 * now x, y, and z are all finite and nonzero; extract signs and
912*25c28e83SPiotr Jasiukajtis 	 * normalize the significands (this will raise the denormal operand
913*25c28e83SPiotr Jasiukajtis 	 * exception if need be)
914*25c28e83SPiotr Jasiukajtis 	 */
915*25c28e83SPiotr Jasiukajtis 	sxy = (xx.i[2] ^ yy.i[2]) & 0x8000;
916*25c28e83SPiotr Jasiukajtis 	sz = zz.i[2] & 0x8000;
917*25c28e83SPiotr Jasiukajtis 	if (!ex) {
918*25c28e83SPiotr Jasiukajtis 		xx.e = x * two63;
919*25c28e83SPiotr Jasiukajtis 		ex = (xx.i[2] & 0x7fff) - 63;
920*25c28e83SPiotr Jasiukajtis 	}
921*25c28e83SPiotr Jasiukajtis 	if (!ey) {
922*25c28e83SPiotr Jasiukajtis 		yy.e = y * two63;
923*25c28e83SPiotr Jasiukajtis 		ey = (yy.i[2] & 0x7fff) - 63;
924*25c28e83SPiotr Jasiukajtis 	}
925*25c28e83SPiotr Jasiukajtis 	if (!ez) {
926*25c28e83SPiotr Jasiukajtis 		zz.e = z * two63;
927*25c28e83SPiotr Jasiukajtis 		ez = (zz.i[2] & 0x7fff) - 63;
928*25c28e83SPiotr Jasiukajtis 	}
929*25c28e83SPiotr Jasiukajtis 
930*25c28e83SPiotr Jasiukajtis 	/*
931*25c28e83SPiotr Jasiukajtis 	 * save the control and status words, mask all exceptions, and
932*25c28e83SPiotr Jasiukajtis 	 * set rounding to 64-bit precision and toward-zero
933*25c28e83SPiotr Jasiukajtis 	 */
934*25c28e83SPiotr Jasiukajtis 	__fenv_getcwsw(&oldcwsw);
935*25c28e83SPiotr Jasiukajtis 	cwsw = (oldcwsw & 0xf0c0ffff) | 0x0f3f0000;
936*25c28e83SPiotr Jasiukajtis 	__fenv_setcwsw(&cwsw);
937*25c28e83SPiotr Jasiukajtis 
938*25c28e83SPiotr Jasiukajtis 	/* multiply x*y to 128 bits */
939*25c28e83SPiotr Jasiukajtis 	exy = ex + ey - 0x3fff;
940*25c28e83SPiotr Jasiukajtis 	xx.i[2] = 0x3fff;
941*25c28e83SPiotr Jasiukajtis 	yy.i[2] = 0x3fff;
942*25c28e83SPiotr Jasiukajtis 	x = xx.e;
943*25c28e83SPiotr Jasiukajtis 	y = yy.e;
944*25c28e83SPiotr Jasiukajtis 	xhi = ((x + twom32) + two32) - two32;
945*25c28e83SPiotr Jasiukajtis 	yhi = ((y + twom32) + two32) - two32;
946*25c28e83SPiotr Jasiukajtis 	xlo = x - xhi;
947*25c28e83SPiotr Jasiukajtis 	ylo = y - yhi;
948*25c28e83SPiotr Jasiukajtis 	x *= y;
949*25c28e83SPiotr Jasiukajtis 	y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo;
950*25c28e83SPiotr Jasiukajtis 	if (x >= two) {
951*25c28e83SPiotr Jasiukajtis 		x *= half;
952*25c28e83SPiotr Jasiukajtis 		y *= half;
953*25c28e83SPiotr Jasiukajtis 		exy++;
954*25c28e83SPiotr Jasiukajtis 	}
955*25c28e83SPiotr Jasiukajtis 
956*25c28e83SPiotr Jasiukajtis 	/* extract the significands */
957*25c28e83SPiotr Jasiukajtis 	xx.e = x;
958*25c28e83SPiotr Jasiukajtis 	xy0 = xx.i[1];
959*25c28e83SPiotr Jasiukajtis 	xy1 = xx.i[0];
960*25c28e83SPiotr Jasiukajtis 	yy.e = t = y + twom32;
961*25c28e83SPiotr Jasiukajtis 	xy2 = yy.i[0];
962*25c28e83SPiotr Jasiukajtis 	yy.e = (y - (t - twom32)) + twom64;
963*25c28e83SPiotr Jasiukajtis 	xy3 = yy.i[0];
964*25c28e83SPiotr Jasiukajtis 	xy4 = 0;
965*25c28e83SPiotr Jasiukajtis 	z0 = zz.i[1];
966*25c28e83SPiotr Jasiukajtis 	z1 = zz.i[0];
967*25c28e83SPiotr Jasiukajtis 	z2 = z3 = z4 = 0;
968*25c28e83SPiotr Jasiukajtis 
969*25c28e83SPiotr Jasiukajtis 	/*
970*25c28e83SPiotr Jasiukajtis 	 * now x*y is represented by sxy, exy, and xy[0-4], and z is
971*25c28e83SPiotr Jasiukajtis 	 * represented likewise; swap if need be so |xy| <= |z|
972*25c28e83SPiotr Jasiukajtis 	 */
973*25c28e83SPiotr Jasiukajtis 	if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 &&
974*25c28e83SPiotr Jasiukajtis 		(xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) {
975*25c28e83SPiotr Jasiukajtis 		e = sxy; sxy = sz; sz = e;
976*25c28e83SPiotr Jasiukajtis 		e = exy; exy = ez; ez = e;
977*25c28e83SPiotr Jasiukajtis 		e = xy0; xy0 = z0; z0 = e;
978*25c28e83SPiotr Jasiukajtis 		e = xy1; xy1 = z1; z1 = e;
979*25c28e83SPiotr Jasiukajtis 		z2 = xy2; xy2 = 0;
980*25c28e83SPiotr Jasiukajtis 		z3 = xy3; xy3 = 0;
981*25c28e83SPiotr Jasiukajtis 	}
982*25c28e83SPiotr Jasiukajtis 
983*25c28e83SPiotr Jasiukajtis 	/* shift the significand of xy keeping a sticky bit */
984*25c28e83SPiotr Jasiukajtis 	e = ez - exy;
985*25c28e83SPiotr Jasiukajtis 	if (e > 130) {
986*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = 0;
987*25c28e83SPiotr Jasiukajtis 		xy4 = 1;
988*25c28e83SPiotr Jasiukajtis 	} else if (e >= 128) {
989*25c28e83SPiotr Jasiukajtis 		sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e));
990*25c28e83SPiotr Jasiukajtis 		xy4 = xy0 >> (e - 128);
991*25c28e83SPiotr Jasiukajtis 		if (sticky)
992*25c28e83SPiotr Jasiukajtis 			xy4 |= 1;
993*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = xy3 = 0;
994*25c28e83SPiotr Jasiukajtis 	} else if (e >= 96) {
995*25c28e83SPiotr Jasiukajtis 		sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e));
996*25c28e83SPiotr Jasiukajtis 		xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
997*25c28e83SPiotr Jasiukajtis 		if (sticky)
998*25c28e83SPiotr Jasiukajtis 			xy4 |= 1;
999*25c28e83SPiotr Jasiukajtis 		xy3 = xy0 >> (e - 96);
1000*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = xy2 = 0;
1001*25c28e83SPiotr Jasiukajtis 	} else if (e >= 64) {
1002*25c28e83SPiotr Jasiukajtis 		sticky = xy3 | ((xy2 << 1) << (95 - e));
1003*25c28e83SPiotr Jasiukajtis 		xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
1004*25c28e83SPiotr Jasiukajtis 		if (sticky)
1005*25c28e83SPiotr Jasiukajtis 			xy4 |= 1;
1006*25c28e83SPiotr Jasiukajtis 		xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
1007*25c28e83SPiotr Jasiukajtis 		xy2 = xy0 >> (e - 64);
1008*25c28e83SPiotr Jasiukajtis 		xy0 = xy1 = 0;
1009*25c28e83SPiotr Jasiukajtis 	} else if (e >= 32) {
1010*25c28e83SPiotr Jasiukajtis 		sticky = (xy3 << 1) << (63 - e);
1011*25c28e83SPiotr Jasiukajtis 		xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
1012*25c28e83SPiotr Jasiukajtis 		if (sticky)
1013*25c28e83SPiotr Jasiukajtis 			xy4 |= 1;
1014*25c28e83SPiotr Jasiukajtis 		xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
1015*25c28e83SPiotr Jasiukajtis 		xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
1016*25c28e83SPiotr Jasiukajtis 		xy1 = xy0 >> (e - 32);
1017*25c28e83SPiotr Jasiukajtis 		xy0 = 0;
1018*25c28e83SPiotr Jasiukajtis 	} else if (e) {
1019*25c28e83SPiotr Jasiukajtis 		xy4 = (xy3 << 1) << (31 - e);
1020*25c28e83SPiotr Jasiukajtis 		xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
1021*25c28e83SPiotr Jasiukajtis 		xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
1022*25c28e83SPiotr Jasiukajtis 		xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
1023*25c28e83SPiotr Jasiukajtis 		xy0 >>= e;
1024*25c28e83SPiotr Jasiukajtis 	}
1025*25c28e83SPiotr Jasiukajtis 
1026*25c28e83SPiotr Jasiukajtis 	/* if this is a magnitude subtract, negate the significand of xy */
1027*25c28e83SPiotr Jasiukajtis 	if (sxy ^ sz) {
1028*25c28e83SPiotr Jasiukajtis 		xy0 = ~xy0;
1029*25c28e83SPiotr Jasiukajtis 		xy1 = ~xy1;
1030*25c28e83SPiotr Jasiukajtis 		xy2 = ~xy2;
1031*25c28e83SPiotr Jasiukajtis 		xy3 = ~xy3;
1032*25c28e83SPiotr Jasiukajtis 		xy4 = -xy4;
1033*25c28e83SPiotr Jasiukajtis 		if (xy4 == 0)
1034*25c28e83SPiotr Jasiukajtis 			if (++xy3 == 0)
1035*25c28e83SPiotr Jasiukajtis 				if (++xy2 == 0)
1036*25c28e83SPiotr Jasiukajtis 					if (++xy1 == 0)
1037*25c28e83SPiotr Jasiukajtis 						xy0++;
1038*25c28e83SPiotr Jasiukajtis 	}
1039*25c28e83SPiotr Jasiukajtis 
1040*25c28e83SPiotr Jasiukajtis 	/* add, propagating carries */
1041*25c28e83SPiotr Jasiukajtis 	z4 += xy4;
1042*25c28e83SPiotr Jasiukajtis 	carry = (z4 < xy4);
1043*25c28e83SPiotr Jasiukajtis 	z3 += xy3;
1044*25c28e83SPiotr Jasiukajtis 	if (carry) {
1045*25c28e83SPiotr Jasiukajtis 		z3++;
1046*25c28e83SPiotr Jasiukajtis 		carry = (z3 <= xy3);
1047*25c28e83SPiotr Jasiukajtis 	} else
1048*25c28e83SPiotr Jasiukajtis 		carry = (z3 < xy3);
1049*25c28e83SPiotr Jasiukajtis 	z2 += xy2;
1050*25c28e83SPiotr Jasiukajtis 	if (carry) {
1051*25c28e83SPiotr Jasiukajtis 		z2++;
1052*25c28e83SPiotr Jasiukajtis 		carry = (z2 <= xy2);
1053*25c28e83SPiotr Jasiukajtis 	} else
1054*25c28e83SPiotr Jasiukajtis 		carry = (z2 < xy2);
1055*25c28e83SPiotr Jasiukajtis 	z1 += xy1;
1056*25c28e83SPiotr Jasiukajtis 	if (carry) {
1057*25c28e83SPiotr Jasiukajtis 		z1++;
1058*25c28e83SPiotr Jasiukajtis 		carry = (z1 <= xy1);
1059*25c28e83SPiotr Jasiukajtis 	} else
1060*25c28e83SPiotr Jasiukajtis 		carry = (z1 < xy1);
1061*25c28e83SPiotr Jasiukajtis 	z0 += xy0;
1062*25c28e83SPiotr Jasiukajtis 	if (carry) {
1063*25c28e83SPiotr Jasiukajtis 		z0++;
1064*25c28e83SPiotr Jasiukajtis 		carry = (z0 <= xy0);
1065*25c28e83SPiotr Jasiukajtis 	} else
1066*25c28e83SPiotr Jasiukajtis 		carry = (z0 < xy0);
1067*25c28e83SPiotr Jasiukajtis 
1068*25c28e83SPiotr Jasiukajtis 	/* for a magnitude subtract, ignore the last carry out */
1069*25c28e83SPiotr Jasiukajtis 	if (sxy ^ sz)
1070*25c28e83SPiotr Jasiukajtis 		carry = 0;
1071*25c28e83SPiotr Jasiukajtis 
1072*25c28e83SPiotr Jasiukajtis 	/* postnormalize and collect rounding information into z2 */
1073*25c28e83SPiotr Jasiukajtis 	if (ez < 1) {
1074*25c28e83SPiotr Jasiukajtis 		/* result is tiny; shift right until exponent is within range */
1075*25c28e83SPiotr Jasiukajtis 		e = 1 - ez;
1076*25c28e83SPiotr Jasiukajtis 		if (e > 67) {
1077*25c28e83SPiotr Jasiukajtis 			z2 = 1;	/* result can't be exactly zero */
1078*25c28e83SPiotr Jasiukajtis 			z0 = z1 = 0;
1079*25c28e83SPiotr Jasiukajtis 		} else if (e >= 64) {
1080*25c28e83SPiotr Jasiukajtis 			sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e));
1081*25c28e83SPiotr Jasiukajtis 			z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e));
1082*25c28e83SPiotr Jasiukajtis 			if (sticky)
1083*25c28e83SPiotr Jasiukajtis 				z2 |= 1;
1084*25c28e83SPiotr Jasiukajtis 			z1 = carry >> (e - 64);
1085*25c28e83SPiotr Jasiukajtis 			z0 = 0;
1086*25c28e83SPiotr Jasiukajtis 		} else if (e >= 32) {
1087*25c28e83SPiotr Jasiukajtis 			sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e));
1088*25c28e83SPiotr Jasiukajtis 			z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
1089*25c28e83SPiotr Jasiukajtis 			if (sticky)
1090*25c28e83SPiotr Jasiukajtis 				z2 |= 1;
1091*25c28e83SPiotr Jasiukajtis 			z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e));
1092*25c28e83SPiotr Jasiukajtis 			z0 = carry >> (e - 32);
1093*25c28e83SPiotr Jasiukajtis 		} else {
1094*25c28e83SPiotr Jasiukajtis 			sticky = z4 | z3 | (z2 << 1) << (31 - e);
1095*25c28e83SPiotr Jasiukajtis 			z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
1096*25c28e83SPiotr Jasiukajtis 			if (sticky)
1097*25c28e83SPiotr Jasiukajtis 				z2 |= 1;
1098*25c28e83SPiotr Jasiukajtis 			z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
1099*25c28e83SPiotr Jasiukajtis 			z0 = (z0 >> e) | ((carry << 1) << (31 - e));
1100*25c28e83SPiotr Jasiukajtis 		}
1101*25c28e83SPiotr Jasiukajtis 		ez = 1;
1102*25c28e83SPiotr Jasiukajtis 	} else if (carry) {
1103*25c28e83SPiotr Jasiukajtis 		/* carry out; shift right by one */
1104*25c28e83SPiotr Jasiukajtis 		sticky = (z2 & 1) | z3 | z4;
1105*25c28e83SPiotr Jasiukajtis 		z2 = (z2 >> 1) | (z1 << 31);
1106*25c28e83SPiotr Jasiukajtis 		if (sticky)
1107*25c28e83SPiotr Jasiukajtis 			z2 |= 1;
1108*25c28e83SPiotr Jasiukajtis 		z1 = (z1 >> 1) | (z0 << 31);
1109*25c28e83SPiotr Jasiukajtis 		z0 = (z0 >> 1) | 0x80000000;
1110*25c28e83SPiotr Jasiukajtis 		ez++;
1111*25c28e83SPiotr Jasiukajtis 	} else {
1112*25c28e83SPiotr Jasiukajtis 		if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) {
1113*25c28e83SPiotr Jasiukajtis 			/*
1114*25c28e83SPiotr Jasiukajtis 			 * borrow/cancellation; shift left as much as
1115*25c28e83SPiotr Jasiukajtis 			 * exponent allows
1116*25c28e83SPiotr Jasiukajtis 			 */
1117*25c28e83SPiotr Jasiukajtis 			while (!z0 && ez >= 33) {
1118*25c28e83SPiotr Jasiukajtis 				z0 = z1;
1119*25c28e83SPiotr Jasiukajtis 				z1 = z2;
1120*25c28e83SPiotr Jasiukajtis 				z2 = z3;
1121*25c28e83SPiotr Jasiukajtis 				z3 = z4;
1122*25c28e83SPiotr Jasiukajtis 				z4 = 0;
1123*25c28e83SPiotr Jasiukajtis 				ez -= 32;
1124*25c28e83SPiotr Jasiukajtis 			}
1125*25c28e83SPiotr Jasiukajtis 			while (z0 < 0x80000000u && ez > 1) {
1126*25c28e83SPiotr Jasiukajtis 				z0 = (z0 << 1) | (z1 >> 31);
1127*25c28e83SPiotr Jasiukajtis 				z1 = (z1 << 1) | (z2 >> 31);
1128*25c28e83SPiotr Jasiukajtis 				z2 = (z2 << 1) | (z3 >> 31);
1129*25c28e83SPiotr Jasiukajtis 				z3 = (z3 << 1) | (z4 >> 31);
1130*25c28e83SPiotr Jasiukajtis 				z4 <<= 1;
1131*25c28e83SPiotr Jasiukajtis 				ez--;
1132*25c28e83SPiotr Jasiukajtis 			}
1133*25c28e83SPiotr Jasiukajtis 		}
1134*25c28e83SPiotr Jasiukajtis 		if (z3 | z4)
1135*25c28e83SPiotr Jasiukajtis 			z2 |= 1;
1136*25c28e83SPiotr Jasiukajtis 	}
1137*25c28e83SPiotr Jasiukajtis 
1138*25c28e83SPiotr Jasiukajtis 	/* get the rounding mode */
1139*25c28e83SPiotr Jasiukajtis 	rm = oldcwsw & 0x0c000000;
1140*25c28e83SPiotr Jasiukajtis 
1141*25c28e83SPiotr Jasiukajtis 	/* adjust exponent if result is subnormal */
1142*25c28e83SPiotr Jasiukajtis 	tinyafter = 0;
1143*25c28e83SPiotr Jasiukajtis 	if (!(z0 & 0x80000000)) {
1144*25c28e83SPiotr Jasiukajtis 		ez = 0;
1145*25c28e83SPiotr Jasiukajtis 		tinyafter = 1;
1146*25c28e83SPiotr Jasiukajtis 		if (!(z0 | z1 | z2)) { /* exact zero */
1147*25c28e83SPiotr Jasiukajtis 			zz.i[2] = rm == FCW_RM ? 0x8000 : 0;
1148*25c28e83SPiotr Jasiukajtis 			zz.i[1] = zz.i[0] = 0;
1149*25c28e83SPiotr Jasiukajtis 			__fenv_setcwsw(&oldcwsw);
1150*25c28e83SPiotr Jasiukajtis 			return (zz.e);
1151*25c28e83SPiotr Jasiukajtis 		}
1152*25c28e83SPiotr Jasiukajtis 	}
1153*25c28e83SPiotr Jasiukajtis 
1154*25c28e83SPiotr Jasiukajtis 	/*
1155*25c28e83SPiotr Jasiukajtis 	 * flip the sense of directed roundings if the result is negative;
1156*25c28e83SPiotr Jasiukajtis 	 * the logic below applies to a positive result
1157*25c28e83SPiotr Jasiukajtis 	 */
1158*25c28e83SPiotr Jasiukajtis 	if (sz && (rm == FCW_RM || rm == FCW_RP))
1159*25c28e83SPiotr Jasiukajtis 		rm = (FCW_RM + FCW_RP) - rm;
1160*25c28e83SPiotr Jasiukajtis 
1161*25c28e83SPiotr Jasiukajtis 	/* round */
1162*25c28e83SPiotr Jasiukajtis 	if (z2) {
1163*25c28e83SPiotr Jasiukajtis 		if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u ||
1164*25c28e83SPiotr Jasiukajtis 			(z2 == 0x80000000u && (z1 & 1))))) {
1165*25c28e83SPiotr Jasiukajtis 			/* round up and renormalize if necessary */
1166*25c28e83SPiotr Jasiukajtis 			if (++z1 == 0) {
1167*25c28e83SPiotr Jasiukajtis 				if (++z0 == 0) {
1168*25c28e83SPiotr Jasiukajtis 					z0 = 0x80000000;
1169*25c28e83SPiotr Jasiukajtis 					ez++;
1170*25c28e83SPiotr Jasiukajtis 				} else if (z0 == 0x80000000) {
1171*25c28e83SPiotr Jasiukajtis 					/* rounded up to smallest normal */
1172*25c28e83SPiotr Jasiukajtis 					ez = 1;
1173*25c28e83SPiotr Jasiukajtis 					if ((rm == FCW_RP && z2 >
1174*25c28e83SPiotr Jasiukajtis 						0x80000000u) || (rm == FCW_RN &&
1175*25c28e83SPiotr Jasiukajtis 						z2 >= 0xc0000000u))
1176*25c28e83SPiotr Jasiukajtis 						/*
1177*25c28e83SPiotr Jasiukajtis 						 * would have rounded up to
1178*25c28e83SPiotr Jasiukajtis 						 * smallest normal even with
1179*25c28e83SPiotr Jasiukajtis 						 * unbounded range
1180*25c28e83SPiotr Jasiukajtis 						 */
1181*25c28e83SPiotr Jasiukajtis 						tinyafter = 0;
1182*25c28e83SPiotr Jasiukajtis 				}
1183*25c28e83SPiotr Jasiukajtis 			}
1184*25c28e83SPiotr Jasiukajtis 		}
1185*25c28e83SPiotr Jasiukajtis 	}
1186*25c28e83SPiotr Jasiukajtis 
1187*25c28e83SPiotr Jasiukajtis 	/* restore the control and status words, check for over/underflow */
1188*25c28e83SPiotr Jasiukajtis 	__fenv_setcwsw(&oldcwsw);
1189*25c28e83SPiotr Jasiukajtis 	if (ez >= 0x7fff) {
1190*25c28e83SPiotr Jasiukajtis 		if (rm == FCW_RN || rm == FCW_RP) {
1191*25c28e83SPiotr Jasiukajtis 			zz.i[2] = sz | 0x7fff;
1192*25c28e83SPiotr Jasiukajtis 			zz.i[1] = 0x80000000;
1193*25c28e83SPiotr Jasiukajtis 			zz.i[0] = 0;
1194*25c28e83SPiotr Jasiukajtis 		} else {
1195*25c28e83SPiotr Jasiukajtis 			zz.i[2] = sz | 0x7ffe;
1196*25c28e83SPiotr Jasiukajtis 			zz.i[1] = 0xffffffff;
1197*25c28e83SPiotr Jasiukajtis 			zz.i[0] = 0xffffffff;
1198*25c28e83SPiotr Jasiukajtis 		}
1199*25c28e83SPiotr Jasiukajtis 		dummy = huge;
1200*25c28e83SPiotr Jasiukajtis 		dummy *= huge;
1201*25c28e83SPiotr Jasiukajtis 	} else {
1202*25c28e83SPiotr Jasiukajtis 		zz.i[2] = sz | ez;
1203*25c28e83SPiotr Jasiukajtis 		zz.i[1] = z0;
1204*25c28e83SPiotr Jasiukajtis 		zz.i[0] = z1;
1205*25c28e83SPiotr Jasiukajtis 
1206*25c28e83SPiotr Jasiukajtis 		/*
1207*25c28e83SPiotr Jasiukajtis 		 * tinyafter => result rounded w/ unbounded range would be tiny,
1208*25c28e83SPiotr Jasiukajtis 		 * z2 nonzero => result delivered is inexact
1209*25c28e83SPiotr Jasiukajtis 		 */
1210*25c28e83SPiotr Jasiukajtis 		if (tinyafter) {
1211*25c28e83SPiotr Jasiukajtis 			dummy = tiny;
1212*25c28e83SPiotr Jasiukajtis 			if (z2)
1213*25c28e83SPiotr Jasiukajtis 				dummy *= tiny;
1214*25c28e83SPiotr Jasiukajtis 			else
1215*25c28e83SPiotr Jasiukajtis 				dummy -= tiny2;
1216*25c28e83SPiotr Jasiukajtis 		} else if (z2) {
1217*25c28e83SPiotr Jasiukajtis 			dummy = huge;
1218*25c28e83SPiotr Jasiukajtis 			dummy += tiny;
1219*25c28e83SPiotr Jasiukajtis 		}
1220*25c28e83SPiotr Jasiukajtis 	}
1221*25c28e83SPiotr Jasiukajtis 
1222*25c28e83SPiotr Jasiukajtis 	return (zz.e);
1223*25c28e83SPiotr Jasiukajtis }
1224*25c28e83SPiotr Jasiukajtis 
1225*25c28e83SPiotr Jasiukajtis #else
1226*25c28e83SPiotr Jasiukajtis #error Unknown architecture
1227*25c28e83SPiotr Jasiukajtis #endif
1228