125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
3025c28e8Piotr Jasiukajtis#pragma weak fmal = __fmal
3125c28e8Piotr Jasiukajtis
3225c28e8Piotr Jasiukajtis#include "libm.h"
3325c28e8Piotr Jasiukajtis#include "fma.h"
3425c28e8Piotr Jasiukajtis#include "fenv_inlines.h"
3525c28e8Piotr Jasiukajtis
3625c28e8Piotr Jasiukajtis#if defined(__sparc)
3725c28e8Piotr Jasiukajtis
3825c28e8Piotr Jasiukajtisstatic const union {
3925c28e8Piotr Jasiukajtis	unsigned i[2];
4025c28e8Piotr Jasiukajtis	double d;
4125c28e8Piotr Jasiukajtis} C[] = {
4225c28e8Piotr Jasiukajtis	{ 0x3fe00000u, 0 },
4325c28e8Piotr Jasiukajtis	{ 0x40000000u, 0 },
4425c28e8Piotr Jasiukajtis	{ 0x3ef00000u, 0 },
4525c28e8Piotr Jasiukajtis	{ 0x3e700000u, 0 },
4625c28e8Piotr Jasiukajtis	{ 0x41300000u, 0 },
4725c28e8Piotr Jasiukajtis	{ 0x3e300000u, 0 },
4825c28e8Piotr Jasiukajtis	{ 0x3b300000u, 0 },
4925c28e8Piotr Jasiukajtis	{ 0x38300000u, 0 },
5025c28e8Piotr Jasiukajtis	{ 0x42300000u, 0 },
5125c28e8Piotr Jasiukajtis	{ 0x3df00000u, 0 },
5225c28e8Piotr Jasiukajtis	{ 0x7fe00000u, 0 },
5325c28e8Piotr Jasiukajtis	{ 0x00100000u, 0 },
5425c28e8Piotr Jasiukajtis	{ 0x00100001u, 0 },
5525c28e8Piotr Jasiukajtis	{ 0, 0 },
5625c28e8Piotr Jasiukajtis	{ 0x7ff00000u, 0 },
5725c28e8Piotr Jasiukajtis	{ 0x7ff00001u, 0 }
5825c28e8Piotr Jasiukajtis};
5925c28e8Piotr Jasiukajtis
6025c28e8Piotr Jasiukajtis#define	half	C[0].d
6125c28e8Piotr Jasiukajtis#define	two	C[1].d
6225c28e8Piotr Jasiukajtis#define	twom16	C[2].d
6325c28e8Piotr Jasiukajtis#define	twom24	C[3].d
6425c28e8Piotr Jasiukajtis#define	two20	C[4].d
6525c28e8Piotr Jasiukajtis#define	twom28	C[5].d
6625c28e8Piotr Jasiukajtis#define	twom76	C[6].d
6725c28e8Piotr Jasiukajtis#define	twom124	C[7].d
6825c28e8Piotr Jasiukajtis#define	two36	C[8].d
6925c28e8Piotr Jasiukajtis#define	twom32	C[9].d
7025c28e8Piotr Jasiukajtis#define	huge	C[10].d
7125c28e8Piotr Jasiukajtis#define	tiny	C[11].d
7225c28e8Piotr Jasiukajtis#define	tiny2	C[12].d
7325c28e8Piotr Jasiukajtis#define	zero	C[13].d
7425c28e8Piotr Jasiukajtis#define	inf	C[14].d
7525c28e8Piotr Jasiukajtis#define	snan	C[15].d
7625c28e8Piotr Jasiukajtis
7725c28e8Piotr Jasiukajtisstatic const unsigned int fsr_rm = 0xc0000000u;
7825c28e8Piotr Jasiukajtis
7925c28e8Piotr Jasiukajtis/*
8025c28e8Piotr Jasiukajtis * fmal for SPARC: 128-bit quad precision, big-endian
8125c28e8Piotr Jasiukajtis */
8225c28e8Piotr Jasiukajtislong double
8325c28e8Piotr Jasiukajtis__fmal(long double x, long double y, long double z) {
8425c28e8Piotr Jasiukajtis	union {
8525c28e8Piotr Jasiukajtis		unsigned int i[4];
8625c28e8Piotr Jasiukajtis		long double q;
8725c28e8Piotr Jasiukajtis	} xx, yy, zz;
8825c28e8Piotr Jasiukajtis	union {
8925c28e8Piotr Jasiukajtis		unsigned int i[2];
9025c28e8Piotr Jasiukajtis		double d;
9125c28e8Piotr Jasiukajtis	} u;
9225c28e8Piotr Jasiukajtis	double dx[5], dy[5], dxy[9], c, s;
9325c28e8Piotr Jasiukajtis	unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
9425c28e8Piotr Jasiukajtis	unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
9525c28e8Piotr Jasiukajtis	unsigned int rm, sticky;
9625c28e8Piotr Jasiukajtis	unsigned int fsr;
9725c28e8Piotr Jasiukajtis	int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit;
9825c28e8Piotr Jasiukajtis	int cx, cy, cz;
9925c28e8Piotr Jasiukajtis	volatile double	dummy;
10025c28e8Piotr Jasiukajtis
10125c28e8Piotr Jasiukajtis	/* extract the high order words of the arguments */
10225c28e8Piotr Jasiukajtis	xx.q = x;
10325c28e8Piotr Jasiukajtis	yy.q = y;
10425c28e8Piotr Jasiukajtis	zz.q = z;
10525c28e8Piotr Jasiukajtis	hx = xx.i[0] & ~0x80000000;
10625c28e8Piotr Jasiukajtis	hy = yy.i[0] & ~0x80000000;
10725c28e8Piotr Jasiukajtis	hz = zz.i[0] & ~0x80000000;
10825c28e8Piotr Jasiukajtis
10925c28e8Piotr Jasiukajtis	/*
11025c28e8Piotr Jasiukajtis	 * distinguish zero, finite nonzero, infinite, and quiet nan
11125c28e8Piotr Jasiukajtis	 * arguments; raise invalid and return for signaling nans
11225c28e8Piotr Jasiukajtis	 */
11325c28e8Piotr Jasiukajtis	if (hx >= 0x7fff0000) {
11425c28e8Piotr Jasiukajtis		if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) {
11525c28e8Piotr Jasiukajtis			if (!(hx & 0x8000)) {
11625c28e8Piotr Jasiukajtis				/* signaling nan, raise invalid */
11725c28e8Piotr Jasiukajtis				dummy = snan;
11825c28e8Piotr Jasiukajtis				dummy += snan;
11925c28e8Piotr Jasiukajtis				xx.i[0] |= 0x8000;
12025c28e8Piotr Jasiukajtis				return (xx.q);
12125c28e8Piotr Jasiukajtis			}
12225c28e8Piotr Jasiukajtis			cx = 3;	/* quiet nan */
12325c28e8Piotr Jasiukajtis		} else
12425c28e8Piotr Jasiukajtis			cx = 2;	/* inf */
12525c28e8Piotr Jasiukajtis	} else if (hx == 0) {
12625c28e8Piotr Jasiukajtis		cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
12725c28e8Piotr Jasiukajtis				/* subnormal or zero */
12825c28e8Piotr Jasiukajtis	} else
12925c28e8Piotr Jasiukajtis		cx = 1;		/* finite nonzero */
13025c28e8Piotr Jasiukajtis
13125c28e8Piotr Jasiukajtis	if (hy >= 0x7fff0000) {
13225c28e8Piotr Jasiukajtis		if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
13325c28e8Piotr Jasiukajtis			if (!(hy & 0x8000)) {
13425c28e8Piotr Jasiukajtis				dummy = snan;
13525c28e8Piotr Jasiukajtis				dummy += snan;
13625c28e8Piotr Jasiukajtis				yy.i[0] |= 0x8000;
13725c28e8Piotr Jasiukajtis				return (yy.q);
13825c28e8Piotr Jasiukajtis			}
13925c28e8Piotr Jasiukajtis			cy = 3;
14025c28e8Piotr Jasiukajtis		} else
14125c28e8Piotr Jasiukajtis			cy = 2;
14225c28e8Piotr Jasiukajtis	} else if (hy == 0) {
14325c28e8Piotr Jasiukajtis		cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
14425c28e8Piotr Jasiukajtis	} else
14525c28e8Piotr Jasiukajtis		cy = 1;
14625c28e8Piotr Jasiukajtis
14725c28e8Piotr Jasiukajtis	if (hz >= 0x7fff0000) {
14825c28e8Piotr Jasiukajtis		if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
14925c28e8Piotr Jasiukajtis			if (!(hz & 0x8000)) {
15025c28e8Piotr Jasiukajtis				dummy = snan;
15125c28e8Piotr Jasiukajtis				dummy += snan;
15225c28e8Piotr Jasiukajtis				zz.i[0] |= 0x8000;
15325c28e8Piotr Jasiukajtis				return (zz.q);
15425c28e8Piotr Jasiukajtis			}
15525c28e8Piotr Jasiukajtis			cz = 3;
15625c28e8Piotr Jasiukajtis		} else
15725c28e8Piotr Jasiukajtis			cz = 2;
15825c28e8Piotr Jasiukajtis	} else if (hz == 0) {
15925c28e8Piotr Jasiukajtis		cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
16025c28e8Piotr Jasiukajtis	} else
16125c28e8Piotr Jasiukajtis		cz = 1;
16225c28e8Piotr Jasiukajtis
16325c28e8Piotr Jasiukajtis	/* get the fsr and clear current exceptions */
16425c28e8Piotr Jasiukajtis	__fenv_getfsr32(&fsr);
16525c28e8Piotr Jasiukajtis	fsr &= ~FSR_CEXC;
16625c28e8Piotr Jasiukajtis
16725c28e8Piotr Jasiukajtis	/* handle all other zero, inf, and nan cases */
16825c28e8Piotr Jasiukajtis	if (cx != 1 || cy != 1 || cz != 1) {
16925c28e8Piotr Jasiukajtis		/* if x or y is a quiet nan, return it */
17025c28e8Piotr Jasiukajtis		if (cx == 3) {
17125c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
17225c28e8Piotr Jasiukajtis			return (x);
17325c28e8Piotr Jasiukajtis		}
17425c28e8Piotr Jasiukajtis		if (cy == 3) {
17525c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
17625c28e8Piotr Jasiukajtis			return (y);
17725c28e8Piotr Jasiukajtis		}
17825c28e8Piotr Jasiukajtis
17925c28e8Piotr Jasiukajtis		/* if x*y is 0*inf, raise invalid and return the default nan */
18025c28e8Piotr Jasiukajtis		if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) {
18125c28e8Piotr Jasiukajtis			dummy = zero;
18225c28e8Piotr Jasiukajtis			dummy *= inf;
18325c28e8Piotr Jasiukajtis			zz.i[0] = 0x7fffffff;
18425c28e8Piotr Jasiukajtis			zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
18525c28e8Piotr Jasiukajtis			return (zz.q);
18625c28e8Piotr Jasiukajtis		}
18725c28e8Piotr Jasiukajtis
18825c28e8Piotr Jasiukajtis		/* if z is a quiet nan, return it */
18925c28e8Piotr Jasiukajtis		if (cz == 3) {
19025c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
19125c28e8Piotr Jasiukajtis			return (z);
19225c28e8Piotr Jasiukajtis		}
19325c28e8Piotr Jasiukajtis
19425c28e8Piotr Jasiukajtis		/*
19525c28e8Piotr Jasiukajtis		 * now none of x, y, or z is nan; handle cases where x or y
19625c28e8Piotr Jasiukajtis		 * is inf
19725c28e8Piotr Jasiukajtis		 */
19825c28e8Piotr Jasiukajtis		if (cx == 2 || cy == 2) {
19925c28e8Piotr Jasiukajtis			/*
20025c28e8Piotr Jasiukajtis			 * if z is also inf, either we have inf-inf or
20125c28e8Piotr Jasiukajtis			 * the result is the same as z depending on signs
20225c28e8Piotr Jasiukajtis			 */
20325c28e8Piotr Jasiukajtis			if (cz == 2) {
20425c28e8Piotr Jasiukajtis				if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
20525c28e8Piotr Jasiukajtis					dummy = inf;
20625c28e8Piotr Jasiukajtis					dummy -= inf;
20725c28e8Piotr Jasiukajtis					zz.i[0] = 0x7fffffff;
20825c28e8Piotr Jasiukajtis					zz.i[1] = zz.i[2] = zz.i[3] =
20925c28e8Piotr Jasiukajtis						0xffffffff;
21025c28e8Piotr Jasiukajtis					return (zz.q);
21125c28e8Piotr Jasiukajtis				}
21225c28e8Piotr Jasiukajtis				__fenv_setfsr32(&fsr);
21325c28e8Piotr Jasiukajtis				return (z);
21425c28e8Piotr Jasiukajtis			}
21525c28e8Piotr Jasiukajtis
21625c28e8Piotr Jasiukajtis			/* otherwise the result is inf with appropriate sign */
21725c28e8Piotr Jasiukajtis			zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) |
21825c28e8Piotr Jasiukajtis				0x7fff0000;
21925c28e8Piotr Jasiukajtis			zz.i[1] = zz.i[2] = zz.i[3] = 0;
22025c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
22125c28e8Piotr Jasiukajtis			return (zz.q);
22225c28e8Piotr Jasiukajtis		}
22325c28e8Piotr Jasiukajtis
22425c28e8Piotr Jasiukajtis		/* if z is inf, return it */
22525c28e8Piotr Jasiukajtis		if (cz == 2) {
22625c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
22725c28e8Piotr Jasiukajtis			return (z);
22825c28e8Piotr Jasiukajtis		}
22925c28e8Piotr Jasiukajtis
23025c28e8Piotr Jasiukajtis		/*
23125c28e8Piotr Jasiukajtis		 * now x, y, and z are all finite; handle cases where x or y
23225c28e8Piotr Jasiukajtis		 * is zero
23325c28e8Piotr Jasiukajtis		 */
23425c28e8Piotr Jasiukajtis		if (cx == 0 || cy == 0) {
23525c28e8Piotr Jasiukajtis			/* either we have 0-0 or the result is the same as z */
23625c28e8Piotr Jasiukajtis			if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
23725c28e8Piotr Jasiukajtis				0) {
23825c28e8Piotr Jasiukajtis				zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
23925c28e8Piotr Jasiukajtis					0;
24025c28e8Piotr Jasiukajtis				__fenv_setfsr32(&fsr);
24125c28e8Piotr Jasiukajtis				return (zz.q);
24225c28e8Piotr Jasiukajtis			}
24325c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
24425c28e8Piotr Jasiukajtis			return (z);
24525c28e8Piotr Jasiukajtis		}
24625c28e8Piotr Jasiukajtis
24725c28e8Piotr Jasiukajtis		/* if we get here, x and y are nonzero finite, z must be zero */
24825c28e8Piotr Jasiukajtis		return (x * y);
24925c28e8Piotr Jasiukajtis	}
25025c28e8Piotr Jasiukajtis
25125c28e8Piotr Jasiukajtis	/*
25225c28e8Piotr Jasiukajtis	 * now x, y, and z are all finite and nonzero; set round-to-
25325c28e8Piotr Jasiukajtis	 * negative-infinity mode
25425c28e8Piotr Jasiukajtis	 */
25525c28e8Piotr Jasiukajtis	__fenv_setfsr32(&fsr_rm);
25625c28e8Piotr Jasiukajtis
25725c28e8Piotr Jasiukajtis	/*
25825c28e8Piotr Jasiukajtis	 * get the signs and exponents and normalize the significands
25925c28e8Piotr Jasiukajtis	 * of x and y
26025c28e8Piotr Jasiukajtis	 */
26125c28e8Piotr Jasiukajtis	sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000;
26225c28e8Piotr Jasiukajtis	ex = hx >> 16;
26325c28e8Piotr Jasiukajtis	hx &= 0xffff;
26425c28e8Piotr Jasiukajtis	if (!ex) {
26525c28e8Piotr Jasiukajtis		if (hx | (xx.i[1] & 0xfffe0000)) {
26625c28e8Piotr Jasiukajtis			ex = 1;
26725c28e8Piotr Jasiukajtis		} else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) {
26825c28e8Piotr Jasiukajtis			hx = xx.i[1];
26925c28e8Piotr Jasiukajtis			xx.i[1] = xx.i[2];
27025c28e8Piotr Jasiukajtis			xx.i[2] = xx.i[3];
27125c28e8Piotr Jasiukajtis			xx.i[3] = 0;
27225c28e8Piotr Jasiukajtis			ex = -31;
27325c28e8Piotr Jasiukajtis		} else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) {
27425c28e8Piotr Jasiukajtis			hx = xx.i[2];
27525c28e8Piotr Jasiukajtis			xx.i[1] = xx.i[3];
27625c28e8Piotr Jasiukajtis			xx.i[2] = xx.i[3] = 0;
27725c28e8Piotr Jasiukajtis			ex = -63;
27825c28e8Piotr Jasiukajtis		} else {
27925c28e8Piotr Jasiukajtis			hx = xx.i[3];
28025c28e8Piotr Jasiukajtis			xx.i[1] = xx.i[2] = xx.i[3] = 0;
28125c28e8Piotr Jasiukajtis			ex = -95;
28225c28e8Piotr Jasiukajtis		}
28325c28e8Piotr Jasiukajtis		while ((hx & 0x10000) == 0) {
28425c28e8Piotr Jasiukajtis			hx = (hx << 1) | (xx.i[1] >> 31);
28525c28e8Piotr Jasiukajtis			xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
28625c28e8Piotr Jasiukajtis			xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
28725c28e8Piotr Jasiukajtis			xx.i[3] <<= 1;
28825c28e8Piotr Jasiukajtis			ex--;
28925c28e8Piotr Jasiukajtis		}
29025c28e8Piotr Jasiukajtis	} else
29125c28e8Piotr Jasiukajtis		hx |= 0x10000;
29225c28e8Piotr Jasiukajtis	ey = hy >> 16;
29325c28e8Piotr Jasiukajtis	hy &= 0xffff;
29425c28e8Piotr Jasiukajtis	if (!ey) {
29525c28e8Piotr Jasiukajtis		if (hy | (yy.i[1] & 0xfffe0000)) {
29625c28e8Piotr Jasiukajtis			ey = 1;
29725c28e8Piotr Jasiukajtis		} else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
29825c28e8Piotr Jasiukajtis			hy = yy.i[1];
29925c28e8Piotr Jasiukajtis			yy.i[1] = yy.i[2];
30025c28e8Piotr Jasiukajtis			yy.i[2] = yy.i[3];
30125c28e8Piotr Jasiukajtis			yy.i[3] = 0;
30225c28e8Piotr Jasiukajtis			ey = -31;
30325c28e8Piotr Jasiukajtis		} else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) {
30425c28e8Piotr Jasiukajtis			hy = yy.i[2];
30525c28e8Piotr Jasiukajtis			yy.i[1] = yy.i[3];
30625c28e8Piotr Jasiukajtis			yy.i[2] = yy.i[3] = 0;
30725c28e8Piotr Jasiukajtis			ey = -63;
30825c28e8Piotr Jasiukajtis		} else {
30925c28e8Piotr Jasiukajtis			hy = yy.i[3];
31025c28e8Piotr Jasiukajtis			yy.i[1] = yy.i[2] = yy.i[3] = 0;
31125c28e8Piotr Jasiukajtis			ey = -95;
31225c28e8Piotr Jasiukajtis		}
31325c28e8Piotr Jasiukajtis		while ((hy & 0x10000) == 0) {
31425c28e8Piotr Jasiukajtis			hy = (hy << 1) | (yy.i[1] >> 31);
31525c28e8Piotr Jasiukajtis			yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
31625c28e8Piotr Jasiukajtis			yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
31725c28e8Piotr Jasiukajtis			yy.i[3] <<= 1;
31825c28e8Piotr Jasiukajtis			ey--;
31925c28e8Piotr Jasiukajtis		}
32025c28e8Piotr Jasiukajtis	} else
32125c28e8Piotr Jasiukajtis		hy |= 0x10000;
32225c28e8Piotr Jasiukajtis	exy = ex + ey - 0x3fff;
32325c28e8Piotr Jasiukajtis
32425c28e8Piotr Jasiukajtis	/* convert the significands of x and y to doubles */
32525c28e8Piotr Jasiukajtis	c = twom16;
32625c28e8Piotr Jasiukajtis	dx[0] = (double) ((int) hx) * c;
32725c28e8Piotr Jasiukajtis	dy[0] = (double) ((int) hy) * c;
32825c28e8Piotr Jasiukajtis
32925c28e8Piotr Jasiukajtis	c *= twom24;
33025c28e8Piotr Jasiukajtis	dx[1] = (double) ((int) (xx.i[1] >> 8)) * c;
33125c28e8Piotr Jasiukajtis	dy[1] = (double) ((int) (yy.i[1] >> 8)) * c;
33225c28e8Piotr Jasiukajtis
33325c28e8Piotr Jasiukajtis	c *= twom24;
33425c28e8Piotr Jasiukajtis	dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) &
33525c28e8Piotr Jasiukajtis	    0xffffff)) * c;
33625c28e8Piotr Jasiukajtis	dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) &
33725c28e8Piotr Jasiukajtis	    0xffffff)) * c;
33825c28e8Piotr Jasiukajtis
33925c28e8Piotr Jasiukajtis	c *= twom24;
34025c28e8Piotr Jasiukajtis	dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) &
34125c28e8Piotr Jasiukajtis	    0xffffff)) * c;
34225c28e8Piotr Jasiukajtis	dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) &
34325c28e8Piotr Jasiukajtis	    0xffffff)) * c;
34425c28e8Piotr Jasiukajtis
34525c28e8Piotr Jasiukajtis	c *= twom24;
34625c28e8Piotr Jasiukajtis	dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c;
34725c28e8Piotr Jasiukajtis	dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c;
34825c28e8Piotr Jasiukajtis
34925c28e8Piotr Jasiukajtis	/* form the "digits" of the product */
35025c28e8Piotr Jasiukajtis	dxy[0] = dx[0] * dy[0];
35125c28e8Piotr Jasiukajtis	dxy[1] = dx[0] * dy[1] + dx[1] * dy[0];
35225c28e8Piotr Jasiukajtis	dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
35325c28e8Piotr Jasiukajtis	dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] +
35425c28e8Piotr Jasiukajtis	    dx[3] * dy[0];
35525c28e8Piotr Jasiukajtis	dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] +
35625c28e8Piotr Jasiukajtis	    dx[3] * dy[1] + dx[4] * dy[0];
35725c28e8Piotr Jasiukajtis	dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] +
35825c28e8Piotr Jasiukajtis	    dx[4] * dy[1];
35925c28e8Piotr Jasiukajtis	dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
36025c28e8Piotr Jasiukajtis	dxy[7] = dx[3] * dy[4] + dx[4] * dy[3];
36125c28e8Piotr Jasiukajtis	dxy[8] = dx[4] * dy[4];
36225c28e8Piotr Jasiukajtis
36325c28e8Piotr Jasiukajtis	/* split odd-numbered terms and combine into even-numbered terms */
36425c28e8Piotr Jasiukajtis	c = (dxy[1] + two20) - two20;
36525c28e8Piotr Jasiukajtis	dxy[0] += c;
36625c28e8Piotr Jasiukajtis	dxy[1] -= c;
36725c28e8Piotr Jasiukajtis	c = (dxy[3] + twom28) - twom28;
36825c28e8Piotr Jasiukajtis	dxy[2] += c + dxy[1];
36925c28e8Piotr Jasiukajtis	dxy[3] -= c;
37025c28e8Piotr Jasiukajtis	c = (dxy[5] + twom76) - twom76;
37125c28e8Piotr Jasiukajtis	dxy[4] += c + dxy[3];
37225c28e8Piotr Jasiukajtis	dxy[5] -= c;
37325c28e8Piotr Jasiukajtis	c = (dxy[7] + twom124) - twom124;
37425c28e8Piotr Jasiukajtis	dxy[6] += c + dxy[5];
37525c28e8Piotr Jasiukajtis	dxy[8] += (dxy[7] - c);
37625c28e8Piotr Jasiukajtis
37725c28e8Piotr Jasiukajtis	/* propagate carries, adjusting the exponent if need be */
37825c28e8Piotr Jasiukajtis	dxy[7] = dxy[6] + dxy[8];
37925c28e8Piotr Jasiukajtis	dxy[5] = dxy[4] + dxy[7];
38025c28e8Piotr Jasiukajtis	dxy[3] = dxy[2] + dxy[5];
38125c28e8Piotr Jasiukajtis	dxy[1] = dxy[0] + dxy[3];
38225c28e8Piotr Jasiukajtis	if (dxy[1] >= two) {
38325c28e8Piotr Jasiukajtis		dxy[0] *= half;
38425c28e8Piotr Jasiukajtis		dxy[1] *= half;
38525c28e8Piotr Jasiukajtis		dxy[2] *= half;
38625c28e8Piotr Jasiukajtis		dxy[3] *= half;
38725c28e8Piotr Jasiukajtis		dxy[4] *= half;
38825c28e8Piotr Jasiukajtis		dxy[5] *= half;
38925c28e8Piotr Jasiukajtis		dxy[6] *= half;
39025c28e8Piotr Jasiukajtis		dxy[7] *= half;
39125c28e8Piotr Jasiukajtis		dxy[8] *= half;
39225c28e8Piotr Jasiukajtis		exy++;
39325c28e8Piotr Jasiukajtis	}
39425c28e8Piotr Jasiukajtis
39525c28e8Piotr Jasiukajtis	/* extract the significand of x*y */
39625c28e8Piotr Jasiukajtis	s = two36;
39725c28e8Piotr Jasiukajtis	u.d = c = dxy[1] + s;
39825c28e8Piotr Jasiukajtis	xy0 = u.i[1];
39925c28e8Piotr Jasiukajtis	c -= s;
40025c28e8Piotr Jasiukajtis	dxy[1] -= c;
40125c28e8Piotr Jasiukajtis	dxy[0] -= c;
40225c28e8Piotr Jasiukajtis
40325c28e8Piotr Jasiukajtis	s *= twom32;
40425c28e8Piotr Jasiukajtis	u.d = c = dxy[1] + s;
40525c28e8Piotr Jasiukajtis	xy1 = u.i[1];
40625c28e8Piotr Jasiukajtis	c -= s;
40725c28e8Piotr Jasiukajtis	dxy[2] += (dxy[0] - c);
40825c28e8Piotr Jasiukajtis	dxy[3] = dxy[2] + dxy[5];
40925c28e8Piotr Jasiukajtis
41025c28e8Piotr Jasiukajtis	s *= twom32;
41125c28e8Piotr Jasiukajtis	u.d = c = dxy[3] + s;
41225c28e8Piotr Jasiukajtis	xy2 = u.i[1];
41325c28e8Piotr Jasiukajtis	c -= s;
41425c28e8Piotr Jasiukajtis	dxy[4] += (dxy[2] - c);
41525c28e8Piotr Jasiukajtis	dxy[5] = dxy[4] + dxy[7];
41625c28e8Piotr Jasiukajtis
41725c28e8Piotr Jasiukajtis	s *= twom32;
41825c28e8Piotr Jasiukajtis	u.d = c = dxy[5] + s;
41925c28e8Piotr Jasiukajtis	xy3 = u.i[1];
42025c28e8Piotr Jasiukajtis	c -= s;
42125c28e8Piotr Jasiukajtis	dxy[4] -= c;
42225c28e8Piotr Jasiukajtis	dxy[5] = dxy[4] + dxy[7];
42325c28e8Piotr Jasiukajtis
42425c28e8Piotr Jasiukajtis	s *= twom32;
42525c28e8Piotr Jasiukajtis	u.d = c = dxy[5] + s;
42625c28e8Piotr Jasiukajtis	xy4 = u.i[1];
42725c28e8Piotr Jasiukajtis	c -= s;
42825c28e8Piotr Jasiukajtis	dxy[6] += (dxy[4] - c);
42925c28e8Piotr Jasiukajtis	dxy[7] = dxy[6] + dxy[8];
43025c28e8Piotr Jasiukajtis
43125c28e8Piotr Jasiukajtis	s *= twom32;
43225c28e8Piotr Jasiukajtis	u.d = c = dxy[7] + s;
43325c28e8Piotr Jasiukajtis	xy5 = u.i[1];
43425c28e8Piotr Jasiukajtis	c -= s;
43525c28e8Piotr Jasiukajtis	dxy[8] += (dxy[6] - c);
43625c28e8Piotr Jasiukajtis
43725c28e8Piotr Jasiukajtis	s *= twom32;
43825c28e8Piotr Jasiukajtis	u.d = c = dxy[8] + s;
43925c28e8Piotr Jasiukajtis	xy6 = u.i[1];
44025c28e8Piotr Jasiukajtis	c -= s;
44125c28e8Piotr Jasiukajtis	dxy[8] -= c;
44225c28e8Piotr Jasiukajtis
44325c28e8Piotr Jasiukajtis	s *= twom32;
44425c28e8Piotr Jasiukajtis	u.d = c = dxy[8] + s;
44525c28e8Piotr Jasiukajtis	xy7 = u.i[1];
44625c28e8Piotr Jasiukajtis
44725c28e8Piotr Jasiukajtis	/* extract the sign, exponent, and significand of z */
44825c28e8Piotr Jasiukajtis	sz = zz.i[0] & 0x80000000;
44925c28e8Piotr Jasiukajtis	ez = hz >> 16;
45025c28e8Piotr Jasiukajtis	z0 = hz & 0xffff;
45125c28e8Piotr Jasiukajtis	if (!ez) {
45225c28e8Piotr Jasiukajtis		if (z0 | (zz.i[1] & 0xfffe0000)) {
45325c28e8Piotr Jasiukajtis			z1 = zz.i[1];
45425c28e8Piotr Jasiukajtis			z2 = zz.i[2];
45525c28e8Piotr Jasiukajtis			z3 = zz.i[3];
45625c28e8Piotr Jasiukajtis			ez = 1;
45725c28e8Piotr Jasiukajtis		} else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) {
45825c28e8Piotr Jasiukajtis			z0 = zz.i[1];
45925c28e8Piotr Jasiukajtis			z1 = zz.i[2];
46025c28e8Piotr Jasiukajtis			z2 = zz.i[3];
46125c28e8Piotr Jasiukajtis			z3 = 0;
46225c28e8Piotr Jasiukajtis			ez = -31;
46325c28e8Piotr Jasiukajtis		} else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) {
46425c28e8Piotr Jasiukajtis			z0 = zz.i[2];
46525c28e8Piotr Jasiukajtis			z1 = zz.i[3];
46625c28e8Piotr Jasiukajtis			z2 = z3 = 0;
46725c28e8Piotr Jasiukajtis			ez = -63;
46825c28e8Piotr Jasiukajtis		} else {
46925c28e8Piotr Jasiukajtis			z0 = zz.i[3];
47025c28e8Piotr Jasiukajtis			z1 = z2 = z3 = 0;
47125c28e8Piotr Jasiukajtis			ez = -95;
47225c28e8Piotr Jasiukajtis		}
47325c28e8Piotr Jasiukajtis		while ((z0 & 0x10000) == 0) {
47425c28e8Piotr Jasiukajtis			z0 = (z0 << 1) | (z1 >> 31);
47525c28e8Piotr Jasiukajtis			z1 = (z1 << 1) | (z2 >> 31);
47625c28e8Piotr Jasiukajtis			z2 = (z2 << 1) | (z3 >> 31);
47725c28e8Piotr Jasiukajtis			z3 <<= 1;
47825c28e8Piotr Jasiukajtis			ez--;
47925c28e8Piotr Jasiukajtis		}
48025c28e8Piotr Jasiukajtis	} else {
48125c28e8Piotr Jasiukajtis		z0 |= 0x10000;
48225c28e8Piotr Jasiukajtis		z1 = zz.i[1];
48325c28e8Piotr Jasiukajtis		z2 = zz.i[2];
48425c28e8Piotr Jasiukajtis		z3 = zz.i[3];
48525c28e8Piotr Jasiukajtis	}
48625c28e8Piotr Jasiukajtis	z4 = z5 = z6 = z7 = 0;
48725c28e8Piotr Jasiukajtis
48825c28e8Piotr Jasiukajtis	/*
48925c28e8Piotr Jasiukajtis	 * now x*y is represented by sxy, exy, and xy[0-7], and z is
49025c28e8Piotr Jasiukajtis	 * represented likewise; swap if need be so |xy| <= |z|
49125c28e8Piotr Jasiukajtis	 */
49225c28e8Piotr Jasiukajtis	if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
49325c28e8Piotr Jasiukajtis		(xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 ||
49425c28e8Piotr Jasiukajtis		(xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
49525c28e8Piotr Jasiukajtis		e = sxy; sxy = sz; sz = e;
49625c28e8Piotr Jasiukajtis		e = exy; exy = ez; ez = e;
49725c28e8Piotr Jasiukajtis		e = xy0; xy0 = z0; z0 = e;
49825c28e8Piotr Jasiukajtis		e = xy1; xy1 = z1; z1 = e;
49925c28e8Piotr Jasiukajtis		e = xy2; xy2 = z2; z2 = e;
50025c28e8Piotr Jasiukajtis		e = xy3; xy3 = z3; z3 = e;
50125c28e8Piotr Jasiukajtis		z4 = xy4; xy4 = 0;
50225c28e8Piotr Jasiukajtis		z5 = xy5; xy5 = 0;
50325c28e8Piotr Jasiukajtis		z6 = xy6; xy6 = 0;
50425c28e8Piotr Jasiukajtis		z7 = xy7; xy7 = 0;
50525c28e8Piotr Jasiukajtis	}
50625c28e8Piotr Jasiukajtis
50725c28e8Piotr Jasiukajtis	/* shift the significand of xy keeping a sticky bit */
50825c28e8Piotr Jasiukajtis	e = ez - exy;
50925c28e8Piotr Jasiukajtis	if (e > 236) {
51025c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
51125c28e8Piotr Jasiukajtis		xy7 = 1;
51225c28e8Piotr Jasiukajtis	} else if (e >= 224) {
51325c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 |
51425c28e8Piotr Jasiukajtis			((xy0 << 1) << (255 - e));
51525c28e8Piotr Jasiukajtis		xy7 = xy0 >> (e - 224);
51625c28e8Piotr Jasiukajtis		if (sticky)
51725c28e8Piotr Jasiukajtis			xy7 |= 1;
51825c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
51925c28e8Piotr Jasiukajtis	} else if (e >= 192) {
52025c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 |
52125c28e8Piotr Jasiukajtis			((xy1 << 1) << (223 - e));
52225c28e8Piotr Jasiukajtis		xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
52325c28e8Piotr Jasiukajtis		if (sticky)
52425c28e8Piotr Jasiukajtis			xy7 |= 1;
52525c28e8Piotr Jasiukajtis		xy6 = xy0 >> (e - 192);
52625c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0;
52725c28e8Piotr Jasiukajtis	} else if (e >= 160) {
52825c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | xy5 | xy4 | xy3 |
52925c28e8Piotr Jasiukajtis			((xy2 << 1) << (191 - e));
53025c28e8Piotr Jasiukajtis		xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
53125c28e8Piotr Jasiukajtis		if (sticky)
53225c28e8Piotr Jasiukajtis			xy7 |= 1;
53325c28e8Piotr Jasiukajtis		xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
53425c28e8Piotr Jasiukajtis		xy5 = xy0 >> (e - 160);
53525c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = xy3 = xy4 = 0;
53625c28e8Piotr Jasiukajtis	} else if (e >= 128) {
53725c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
53825c28e8Piotr Jasiukajtis		xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
53925c28e8Piotr Jasiukajtis		if (sticky)
54025c28e8Piotr Jasiukajtis			xy7 |= 1;
54125c28e8Piotr Jasiukajtis		xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
54225c28e8Piotr Jasiukajtis		xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
54325c28e8Piotr Jasiukajtis		xy4 = xy0 >> (e - 128);
54425c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = xy3 = 0;
54525c28e8Piotr Jasiukajtis	} else if (e >= 96) {
54625c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
54725c28e8Piotr Jasiukajtis		xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
54825c28e8Piotr Jasiukajtis		if (sticky)
54925c28e8Piotr Jasiukajtis			xy7 |= 1;
55025c28e8Piotr Jasiukajtis		xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
55125c28e8Piotr Jasiukajtis		xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
55225c28e8Piotr Jasiukajtis		xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
55325c28e8Piotr Jasiukajtis		xy3 = xy0 >> (e - 96);
55425c28e8Piotr Jasiukajtis		xy0 = xy1 = xy2 = 0;
55525c28e8Piotr Jasiukajtis	} else if (e >= 64) {
55625c28e8Piotr Jasiukajtis		sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
55725c28e8Piotr Jasiukajtis		xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
55825c28e8Piotr Jasiukajtis		if (sticky)
55925c28e8Piotr Jasiukajtis			xy7 |= 1;
56025c28e8Piotr Jasiukajtis		xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
56125c28e8Piotr Jasiukajtis		xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
56225c28e8Piotr Jasiukajtis		xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
56325c28e8Piotr Jasiukajtis		xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
56425c28e8Piotr Jasiukajtis		xy2 = xy0 >> (e - 64);
56525c28e8Piotr Jasiukajtis		xy0 = xy1 = 0;
56625c28e8Piotr Jasiukajtis	} else if (e >= 32) {
56725c28e8Piotr Jasiukajtis		sticky = xy7 | ((xy6 << 1) << (63 - e));
56825c28e8Piotr Jasiukajtis		xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
56925c28e8Piotr Jasiukajtis		if (sticky)
57025c28e8Piotr Jasiukajtis			xy7 |= 1;
57125c28e8Piotr Jasiukajtis		xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
57225c28e8Piotr Jasiukajtis		xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
57325c28e8Piotr Jasiukajtis		xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
57425c28e8Piotr Jasiukajtis		xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
57525c28e8Piotr Jasiukajtis		xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
57625c28e8Piotr Jasiukajtis		xy1 = xy0 >> (e - 32);
57725c28e8Piotr Jasiukajtis		xy0 = 0;
57825c28e8Piotr Jasiukajtis	} else if (e) {
57925c28e8Piotr Jasiukajtis		sticky = (xy7 << 1) << (31 - e);
58025c28e8Piotr Jasiukajtis		xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
58125c28e8Piotr Jasiukajtis		if (sticky)
58225c28e8Piotr Jasiukajtis			xy7 |= 1;
58325c28e8Piotr Jasiukajtis		xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
58425c28e8Piotr Jasiukajtis		xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
58525c28e8Piotr Jasiukajtis		xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
58625c28e8Piotr Jasiukajtis		xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
58725c28e8Piotr Jasiukajtis		xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
58825c28e8Piotr Jasiukajtis		xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
58925c28e8Piotr Jasiukajtis		xy0 >>= e;
59025c28e8Piotr Jasiukajtis	}
59125c28e8Piotr Jasiukajtis
59225c28e8Piotr Jasiukajtis	/* if this is a magnitude subtract, negate the significand of xy */
59325c28e8Piotr Jasiukajtis	if (sxy ^ sz) {
59425c28e8Piotr Jasiukajtis		xy0 = ~xy0;
59525c28e8Piotr Jasiukajtis		xy1 = ~xy1;
59625c28e8Piotr Jasiukajtis		xy2 = ~xy2;
59725c28e8Piotr Jasiukajtis		xy3 = ~xy3;
59825c28e8Piotr Jasiukajtis		xy4 = ~xy4;
59925c28e8Piotr Jasiukajtis		xy5 = ~xy5;
60025c28e8Piotr Jasiukajtis		xy6 = ~xy6;
60125c28e8Piotr Jasiukajtis		xy7 = -xy7;
60225c28e8Piotr Jasiukajtis		if (xy7 == 0)
60325c28e8Piotr Jasiukajtis			if (++xy6 == 0)
60425c28e8Piotr Jasiukajtis				if (++xy5 == 0)
60525c28e8Piotr Jasiukajtis					if (++xy4 == 0)
60625c28e8Piotr Jasiukajtis						if (++xy3 == 0)
60725c28e8Piotr Jasiukajtis							if (++xy2 == 0)
60825c28e8Piotr Jasiukajtis								if (++xy1 == 0)
60925c28e8Piotr Jasiukajtis									xy0++;
61025c28e8Piotr Jasiukajtis	}
61125c28e8Piotr Jasiukajtis
61225c28e8Piotr Jasiukajtis	/* add, propagating carries */
61325c28e8Piotr Jasiukajtis	z7 += xy7;
61425c28e8Piotr Jasiukajtis	e = (z7 < xy7);
61525c28e8Piotr Jasiukajtis	z6 += xy6;
61625c28e8Piotr Jasiukajtis	if (e) {
61725c28e8Piotr Jasiukajtis		z6++;
61825c28e8Piotr Jasiukajtis		e = (z6 <= xy6);
61925c28e8Piotr Jasiukajtis	} else
62025c28e8Piotr Jasiukajtis		e = (z6 < xy6);
62125c28e8Piotr Jasiukajtis	z5 += xy5;
62225c28e8Piotr Jasiukajtis	if (e) {
62325c28e8Piotr Jasiukajtis		z5++;
62425c28e8Piotr Jasiukajtis		e = (z5 <= xy5);
62525c28e8Piotr Jasiukajtis	} else
62625c28e8Piotr Jasiukajtis		e = (z5 < xy5);
62725c28e8Piotr Jasiukajtis	z4 += xy4;
62825c28e8Piotr Jasiukajtis	if (e) {
62925c28e8Piotr Jasiukajtis		z4++;
63025c28e8Piotr Jasiukajtis		e = (z4 <= xy4);
63125c28e8Piotr Jasiukajtis	} else
63225c28e8Piotr Jasiukajtis		e = (z4 < xy4);
63325c28e8Piotr Jasiukajtis	z3 += xy3;
63425c28e8Piotr Jasiukajtis	if (e) {
63525c28e8Piotr Jasiukajtis		z3++;
63625c28e8Piotr Jasiukajtis		e = (z3 <= xy3);
63725c28e8Piotr Jasiukajtis	} else
63825c28e8Piotr Jasiukajtis		e = (z3 < xy3);
63925c28e8Piotr Jasiukajtis	z2 += xy2;
64025c28e8Piotr Jasiukajtis	if (e) {
64125c28e8Piotr Jasiukajtis		z2++;
64225c28e8Piotr Jasiukajtis		e = (z2 <= xy2);
64325c28e8Piotr Jasiukajtis	} else
64425c28e8Piotr Jasiukajtis		e = (z2 < xy2);
64525c28e8Piotr Jasiukajtis	z1 += xy1;
64625c28e8Piotr Jasiukajtis	if (e) {
64725c28e8Piotr Jasiukajtis		z1++;
64825c28e8Piotr Jasiukajtis		e = (z1 <= xy1);
64925c28e8Piotr Jasiukajtis	} else
65025c28e8Piotr Jasiukajtis		e = (z1 < xy1);
65125c28e8Piotr Jasiukajtis	z0 += xy0;
65225c28e8Piotr Jasiukajtis	if (e)
65325c28e8Piotr Jasiukajtis		z0++;
65425c28e8Piotr Jasiukajtis
65525c28e8Piotr Jasiukajtis	/* postnormalize and collect rounding information into z4 */
65625c28e8Piotr Jasiukajtis	if (ez < 1) {
65725c28e8Piotr Jasiukajtis		/* result is tiny; shift right until exponent is within range */
65825c28e8Piotr Jasiukajtis		e = 1 - ez;
65925c28e8Piotr Jasiukajtis		if (e > 116) {
66025c28e8Piotr Jasiukajtis			z4 = 1; /* result can't be exactly zero */
66125c28e8Piotr Jasiukajtis			z0 = z1 = z2 = z3 = 0;
66225c28e8Piotr Jasiukajtis		} else if (e >= 96) {
66325c28e8Piotr Jasiukajtis			sticky = z7 | z6 | z5 | z4 | z3 | z2 |
66425c28e8Piotr Jasiukajtis				((z1 << 1) << (127 - e));
66525c28e8Piotr Jasiukajtis			z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
66625c28e8Piotr Jasiukajtis			if (sticky)
66725c28e8Piotr Jasiukajtis				z4 |= 1;
66825c28e8Piotr Jasiukajtis			z3 = z0 >> (e - 96);
66925c28e8Piotr Jasiukajtis			z0 = z1 = z2 = 0;
67025c28e8Piotr Jasiukajtis		} else if (e >= 64) {
67125c28e8Piotr Jasiukajtis			sticky = z7 | z6 | z5 | z4 | z3 |
67225c28e8Piotr Jasiukajtis				((z2 << 1) << (95 - e));
67325c28e8Piotr Jasiukajtis			z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
67425c28e8Piotr Jasiukajtis			if (sticky)
67525c28e8Piotr Jasiukajtis				z4 |= 1;
67625c28e8Piotr Jasiukajtis			z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
67725c28e8Piotr Jasiukajtis			z2 = z0 >> (e - 64);
67825c28e8Piotr Jasiukajtis			z0 = z1 = 0;
67925c28e8Piotr Jasiukajtis		} else if (e >= 32) {
68025c28e8Piotr Jasiukajtis			sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
68125c28e8Piotr Jasiukajtis			z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
68225c28e8Piotr Jasiukajtis			if (sticky)
68325c28e8Piotr Jasiukajtis				z4 |= 1;
68425c28e8Piotr Jasiukajtis			z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
68525c28e8Piotr Jasiukajtis			z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
68625c28e8Piotr Jasiukajtis			z1 = z0 >> (e - 32);
68725c28e8Piotr Jasiukajtis			z0 = 0;
68825c28e8Piotr Jasiukajtis		} else {
68925c28e8Piotr Jasiukajtis			sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
69025c28e8Piotr Jasiukajtis			z4 = (z4 >> e) | ((z3 << 1) << (31 - e));
69125c28e8Piotr Jasiukajtis			if (sticky)
69225c28e8Piotr Jasiukajtis				z4 |= 1;
69325c28e8Piotr Jasiukajtis			z3 = (z3 >> e) | ((z2 << 1) << (31 - e));
69425c28e8Piotr Jasiukajtis			z2 = (z2 >> e) | ((z1 << 1) << (31 - e));
69525c28e8Piotr Jasiukajtis			z1 = (z1 >> e) | ((z0 << 1) << (31 - e));
69625c28e8Piotr Jasiukajtis			z0 >>= e;
69725c28e8Piotr Jasiukajtis		}
69825c28e8Piotr Jasiukajtis		ez = 1;
69925c28e8Piotr Jasiukajtis	} else if (z0 >= 0x20000) {
70025c28e8Piotr Jasiukajtis		/* carry out; shift right by one */
70125c28e8Piotr Jasiukajtis		sticky = (z4 & 1) | z5 | z6 | z7;
70225c28e8Piotr Jasiukajtis		z4 = (z4 >> 1) | (z3 << 31);
70325c28e8Piotr Jasiukajtis		if (sticky)
70425c28e8Piotr Jasiukajtis			z4 |= 1;
70525c28e8Piotr Jasiukajtis		z3 = (z3 >> 1) | (z2 << 31);
70625c28e8Piotr Jasiukajtis		z2 = (z2 >> 1) | (z1 << 31);
70725c28e8Piotr Jasiukajtis		z1 = (z1 >> 1) | (z0 << 31);
70825c28e8Piotr Jasiukajtis		z0 >>= 1;
70925c28e8Piotr Jasiukajtis		ez++;
71025c28e8Piotr Jasiukajtis	} else {
71125c28e8Piotr Jasiukajtis		if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7)
71225c28e8Piotr Jasiukajtis			!= 0) {
71325c28e8Piotr Jasiukajtis			/*
71425c28e8Piotr Jasiukajtis			 * borrow/cancellation; shift left as much as
71525c28e8Piotr Jasiukajtis			 * exponent allows
71625c28e8Piotr Jasiukajtis			 */
71725c28e8Piotr Jasiukajtis			while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
71825c28e8Piotr Jasiukajtis				z0 = z1;
71925c28e8Piotr Jasiukajtis				z1 = z2;
72025c28e8Piotr Jasiukajtis				z2 = z3;
72125c28e8Piotr Jasiukajtis				z3 = z4;
72225c28e8Piotr Jasiukajtis				z4 = z5;
72325c28e8Piotr Jasiukajtis				z5 = z6;
72425c28e8Piotr Jasiukajtis				z6 = z7;
72525c28e8Piotr Jasiukajtis				z7 = 0;
72625c28e8Piotr Jasiukajtis				ez -= 32;
72725c28e8Piotr Jasiukajtis			}
72825c28e8Piotr Jasiukajtis			while (z0 < 0x10000 && ez > 1) {
72925c28e8Piotr Jasiukajtis				z0 = (z0 << 1) | (z1 >> 31);
73025c28e8Piotr Jasiukajtis				z1 = (z1 << 1) | (z2 >> 31);
73125c28e8Piotr Jasiukajtis				z2 = (z2 << 1) | (z3 >> 31);
73225c28e8Piotr Jasiukajtis				z3 = (z3 << 1) | (z4 >> 31);
73325c28e8Piotr Jasiukajtis				z4 = (z4 << 1) | (z5 >> 31);
73425c28e8Piotr Jasiukajtis				z5 = (z5 << 1) | (z6 >> 31);
73525c28e8Piotr Jasiukajtis				z6 = (z6 << 1) | (z7 >> 31);
73625c28e8Piotr Jasiukajtis				z7 <<= 1;
73725c28e8Piotr Jasiukajtis				ez--;
73825c28e8Piotr Jasiukajtis			}
73925c28e8Piotr Jasiukajtis		}
74025c28e8Piotr Jasiukajtis		if (z5 | z6 | z7)
74125c28e8Piotr Jasiukajtis			z4 |= 1;
74225c28e8Piotr Jasiukajtis	}
74325c28e8Piotr Jasiukajtis
74425c28e8Piotr Jasiukajtis	/* get the rounding mode */
74525c28e8Piotr Jasiukajtis	rm = fsr >> 30;
74625c28e8Piotr Jasiukajtis
74725c28e8Piotr Jasiukajtis	/* strip off the integer bit, if there is one */
74825c28e8Piotr Jasiukajtis	ibit = z0 & 0x10000;
74925c28e8Piotr Jasiukajtis	if (ibit)
75025c28e8Piotr Jasiukajtis		z0 -= 0x10000;
75125c28e8Piotr Jasiukajtis	else {
75225c28e8Piotr Jasiukajtis		ez = 0;
75325c28e8Piotr Jasiukajtis		if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
75425c28e8Piotr Jasiukajtis			zz.i[0] = rm == FSR_RM ? 0x80000000 : 0;
75525c28e8Piotr Jasiukajtis			zz.i[1] = zz.i[2] = zz.i[3] = 0;
75625c28e8Piotr Jasiukajtis			__fenv_setfsr32(&fsr);
75725c28e8Piotr Jasiukajtis			return (zz.q);
75825c28e8Piotr Jasiukajtis		}
75925c28e8Piotr Jasiukajtis	}
76025c28e8Piotr Jasiukajtis
76125c28e8Piotr Jasiukajtis	/*
76225c28e8Piotr Jasiukajtis	 * flip the sense of directed roundings if the result is negative;
76325c28e8Piotr Jasiukajtis	 * the logic below applies to a positive result
76425c28e8Piotr Jasiukajtis	 */
76525c28e8Piotr Jasiukajtis	if (sz)
76625c28e8Piotr Jasiukajtis		rm ^= rm >> 1;
76725c28e8Piotr Jasiukajtis
76825c28e8Piotr Jasiukajtis	/* round and raise exceptions */
76925c28e8Piotr Jasiukajtis	if (z4) {
77025c28e8Piotr Jasiukajtis		fsr |= FSR_NXC;
77125c28e8Piotr Jasiukajtis
77225c28e8Piotr Jasiukajtis		/* decide whether to round the fraction up */
77325c28e8Piotr Jasiukajtis		if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u ||
77425c28e8Piotr Jasiukajtis			(z4 == 0x80000000u && (z3 & 1))))) {
77525c28e8Piotr Jasiukajtis			/* round up and renormalize if necessary */
77625c28e8Piotr Jasiukajtis			if (++z3 == 0)
77725c28e8Piotr Jasiukajtis				if (++z2 == 0)
77825c28e8Piotr Jasiukajtis					if (++z1 == 0)
77925c28e8Piotr Jasiukajtis						if (++z0 == 0x10000) {
78025c28e8Piotr Jasiukajtis							z0 = 0;
78125c28e8Piotr Jasiukajtis							ez++;
78225c28e8Piotr Jasiukajtis						}
78325c28e8Piotr Jasiukajtis		}
78425c28e8Piotr Jasiukajtis	}
78525c28e8Piotr Jasiukajtis
78625c28e8Piotr Jasiukajtis	/* check for under/overflow */
78725c28e8Piotr Jasiukajtis	if (ez >= 0x7fff) {
78825c28e8Piotr Jasiukajtis		if (rm == FSR_RN || rm == FSR_RP) {
78925c28e8Piotr Jasiukajtis			zz.i[0] = sz | 0x7fff0000;
79025c28e8Piotr Jasiukajtis			zz.i[1] = zz.i[2] = zz.i[3] = 0;
79125c28e8Piotr Jasiukajtis		} else {
79225c28e8Piotr Jasiukajtis			zz.i[0] = sz | 0x7ffeffff;
79325c28e8Piotr Jasiukajtis			zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff;
79425c28e8Piotr Jasiukajtis		}
79525c28e8Piotr Jasiukajtis		fsr |= FSR_OFC | FSR_NXC;
79625c28e8Piotr Jasiukajtis	} else {
79725c28e8Piotr Jasiukajtis		zz.i[0] = sz | (ez << 16) | z0;
79825c28e8Piotr Jasiukajtis		zz.i[1] = z1;
79925c28e8Piotr Jasiukajtis		zz.i[2] = z2;
80025c28e8Piotr Jasiukajtis		zz.i[3] = z3;
80125c28e8Piotr Jasiukajtis
80225c28e8Piotr Jasiukajtis		/*
80325c28e8Piotr Jasiukajtis		 * !ibit => exact result was tiny before rounding,
80425c28e8Piotr Jasiukajtis		 * z4 nonzero => result delivered is inexact
80525c28e8Piotr Jasiukajtis		 */
80625c28e8Piotr Jasiukajtis		if (!ibit) {
80725c28e8Piotr Jasiukajtis			if (z4)
80825c28e8Piotr Jasiukajtis				fsr |= FSR_UFC | FSR_NXC;
80925c28e8Piotr Jasiukajtis			else if (fsr & FSR_UFM)
81025c28e8Piotr Jasiukajtis				fsr |= FSR_UFC;
81125c28e8Piotr Jasiukajtis		}
81225c28e8Piotr Jasiukajtis	}
81325c28e8Piotr Jasiukajtis
81425c28e8Piotr Jasiukajtis	/* restore the fsr and emulate exceptions as needed */
81525c28e8Piotr Jasiukajtis	if ((fsr & FSR_CEXC) & (fsr >> 23)) {
81625c28e8Piotr Jasiukajtis		__fenv_setfsr32(&fsr);
81725c28e8Piotr Jasiukajtis		if (fsr & FSR_OFC) {
81825c28e8Piotr Jasiukajtis			dummy = huge;
81925c28e8Piotr Jasiukajtis			dummy *= huge;
82025c28e8Piotr Jasiukajtis		} else if (fsr & FSR_UFC) {
82125c28e8Piotr Jasiukajtis			dummy = tiny;
82225c28e8Piotr Jasiukajtis			if (fsr & FSR_NXC)
82325c28e8Piotr Jasiukajtis				dummy *= tiny;
82425c28e8Piotr Jasiukajtis			else
82525c28e8Piotr Jasiukajtis				dummy -= tiny2;
82625c28e8Piotr Jasiukajtis		} else {
82725c28e8Piotr Jasiukajtis			dummy = huge;
82825c28e8Piotr Jasiukajtis			dummy += tiny;
82925c28e8Piotr Jasiukajtis		}
83025c28e8Piotr Jasiukajtis	} else {
83125c28e8Piotr Jasiukajtis		fsr |= (fsr & 0x1f) << 5;
83225c28e8Piotr Jasiukajtis		__fenv_setfsr32(&fsr);
83325c28e8Piotr Jasiukajtis	}
83425c28e8Piotr Jasiukajtis	return (zz.q);
83525c28e8Piotr Jasiukajtis}
83625c28e8Piotr Jasiukajtis
83725c28e8Piotr Jasiukajtis#elif defined(__x86)
83825c28e8Piotr Jasiukajtis
83925c28e8Piotr Jasiukajtisstatic const union {
84025c28e8Piotr Jasiukajtis	unsigned i[2];
84125c28e8Piotr Jasiukajtis	double d;
84225c28e8Piotr Jasiukajtis} C[] = {
84325c28e8Piotr Jasiukajtis	{ 0, 0x3fe00000u },
84425c28e8Piotr Jasiukajtis	{ 0, 0x40000000u },
84525c28e8Piotr Jasiukajtis	{ 0, 0x3df00000u },
84625c28e8Piotr Jasiukajtis	{ 0, 0x3bf00000u },
84725c28e8Piotr Jasiukajtis	{ 0, 0x41f00000u },
84825c28e8Piotr Jasiukajtis	{ 0, 0x43e00000u },
84925c28e8Piotr Jasiukajtis	{ 0, 0x7fe00000u },
85025c28e8Piotr Jasiukajtis	{ 0, 0x00100000u },
85125c28e8Piotr Jasiukajtis	{ 0, 0x00100001u }
85225c28e8Piotr Jasiukajtis};
85325c28e8Piotr Jasiukajtis
85425c28e8Piotr Jasiukajtis#define	half	C[0].d
85525c28e8Piotr Jasiukajtis#define	two	C[1].d
85625c28e8Piotr Jasiukajtis#define	twom32	C[2].d
85725c28e8Piotr Jasiukajtis#define	twom64	C[3].d
85825c28e8Piotr Jasiukajtis#define	two32	C[4].d
85925c28e8Piotr Jasiukajtis#define	two63	C[5].d
86025c28e8Piotr Jasiukajtis#define	huge	C[6].d
86125c28e8Piotr Jasiukajtis#define	tiny	C[7].d
86225c28e8Piotr Jasiukajtis#define	tiny2	C[8].d
86325c28e8Piotr Jasiukajtis
86425c28e8Piotr Jasiukajtis#if defined(__amd64)
86525c28e8Piotr Jasiukajtis#define	NI	4
86625c28e8Piotr Jasiukajtis#else
86725c28e8Piotr Jasiukajtis#define	NI	3
86825c28e8Piotr Jasiukajtis#endif
86925c28e8Piotr Jasiukajtis
87025c28e8Piotr Jasiukajtis/*
87125c28e8Piotr Jasiukajtis * fmal for x86: 80-bit extended double precision, little-endian
87225c28e8Piotr Jasiukajtis */
87325c28e8Piotr Jasiukajtislong double
87425c28e8Piotr Jasiukajtis__fmal(long double x, long double y, long double z) {
87525c28e8Piotr Jasiukajtis	union {
87625c28e8Piotr Jasiukajtis		unsigned i[NI];
87725c28e8Piotr Jasiukajtis		long double e;
87825c28e8Piotr Jasiukajtis	} xx, yy, zz;
87925c28e8Piotr Jasiukajtis	long double xhi, yhi, xlo, ylo, t;
88025c28e8Piotr Jasiukajtis	unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
88125c28e8Piotr Jasiukajtis	unsigned oldcwsw, cwsw, rm, sticky, carry;
88225c28e8Piotr Jasiukajtis	int ex, ey, ez, exy, sxy, sz, e, tinyafter;
88325c28e8Piotr Jasiukajtis	volatile double	dummy;
88425c28e8Piotr Jasiukajtis
88525c28e8Piotr Jasiukajtis	/* extract the exponents of the arguments */
88625c28e8Piotr Jasiukajtis	xx.e = x;
88725c28e8Piotr Jasiukajtis	yy.e = y;
88825c28e8Piotr Jasiukajtis	zz.e = z;
88925c28e8Piotr Jasiukajtis	ex = xx.i[2] & 0x7fff;
89025c28e8Piotr Jasiukajtis	ey = yy.i[2] & 0x7fff;
89125c28e8Piotr Jasiukajtis	ez = zz.i[2] & 0x7fff;
89225c28e8Piotr Jasiukajtis
89325c28e8Piotr Jasiukajtis	/* dispense with inf, nan, and zero cases */
89425c28e8Piotr Jasiukajtis	if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 ||
89525c28e8Piotr Jasiukajtis		(ey | yy.i[1] | yy.i[0]) == 0)	/* x or y is inf, nan, or 0 */
89625c28e8Piotr Jasiukajtis		return (x * y + z);
89725c28e8Piotr Jasiukajtis
89825c28e8Piotr Jasiukajtis	if (ez == 0x7fff)			/* z is inf or nan */
89925c28e8Piotr Jasiukajtis		return (x + z);	/* avoid spurious under/overflow in x * y */
90025c28e8Piotr Jasiukajtis
90125c28e8Piotr Jasiukajtis	if ((ez | zz.i[1] | zz.i[0]) == 0)	/* z is zero */
90225c28e8Piotr Jasiukajtis		/*
90325c28e8Piotr Jasiukajtis		 * x * y isn't zero but could underflow to zero,
90425c28e8Piotr Jasiukajtis		 * so don't add z, lest we perturb the sign
90525c28e8Piotr Jasiukajtis		 */
90625c28e8Piotr Jasiukajtis		return (x * y);
90725c28e8Piotr Jasiukajtis
90825c28e8Piotr Jasiukajtis	/*
90925c28e8Piotr Jasiukajtis	 * now x, y, and z are all finite and nonzero; extract signs and
91025c28e8Piotr Jasiukajtis	 * normalize the significands (this will raise the denormal operand
91125c28e8Piotr Jasiukajtis	 * exception if need be)
91225c28e8Piotr Jasiukajtis	 */
91325c28e8Piotr Jasiukajtis	sxy = (xx.i[2] ^ yy.i[2]) & 0x8000;
91425c28e8Piotr Jasiukajtis	sz = zz.i[2] & 0x8000;
91525c28e8Piotr Jasiukajtis	if (!ex) {
91625c28e8Piotr Jasiukajtis		xx.e = x * two63;
91725c28e8Piotr Jasiukajtis		ex = (xx.i[2] & 0x7fff) - 63;
91825c28e8Piotr Jasiukajtis	}
91925c28e8Piotr Jasiukajtis	if (!ey) {
92025c28e8Piotr Jasiukajtis		yy.e = y * two63;
92125c28e8Piotr Jasiukajtis		ey = (yy.i[2] & 0x7fff) - 63;
92225c28e8Piotr Jasiukajtis	}
92325c28e8Piotr Jasiukajtis	if (!ez) {
92425c28e8Piotr Jasiukajtis		zz.e = z * two63;
925