xref: /illumos-gate/usr/src/lib/libm/common/Q/sqrtl.c (revision ddc0e0b5)
125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis  * CDDL HEADER START
325c28e83SPiotr Jasiukajtis  *
425c28e83SPiotr Jasiukajtis  * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis  * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis  * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis  *
825c28e83SPiotr Jasiukajtis  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis  * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis  * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis  * and limitations under the License.
1225c28e83SPiotr Jasiukajtis  *
1325c28e83SPiotr Jasiukajtis  * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis  * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis  * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis  * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis  *
1925c28e83SPiotr Jasiukajtis  * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis  */
2125c28e83SPiotr Jasiukajtis 
2225c28e83SPiotr Jasiukajtis /*
2325c28e83SPiotr Jasiukajtis  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e83SPiotr Jasiukajtis  */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e83SPiotr Jasiukajtis  * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis  */
2925c28e83SPiotr Jasiukajtis 
30*ddc0e0b5SRichard Lowe #pragma weak __sqrtl = sqrtl
3125c28e83SPiotr Jasiukajtis 
3225c28e83SPiotr Jasiukajtis #include "libm.h"
3325c28e83SPiotr Jasiukajtis #include "longdouble.h"
3425c28e83SPiotr Jasiukajtis 
3525c28e83SPiotr Jasiukajtis extern int __swapTE(int);
3625c28e83SPiotr Jasiukajtis extern int __swapEX(int);
3725c28e83SPiotr Jasiukajtis extern enum fp_direction_type __swapRD(enum fp_direction_type);
3825c28e83SPiotr Jasiukajtis 
3925c28e83SPiotr Jasiukajtis /*
4025c28e83SPiotr Jasiukajtis  * in struct longdouble, msw consists of
4125c28e83SPiotr Jasiukajtis  *	unsigned short	sgn:1;
4225c28e83SPiotr Jasiukajtis  *	unsigned short	exp:15;
4325c28e83SPiotr Jasiukajtis  *	unsigned short	frac1:16;
4425c28e83SPiotr Jasiukajtis  */
4525c28e83SPiotr Jasiukajtis 
4625c28e83SPiotr Jasiukajtis #ifdef __LITTLE_ENDIAN
4725c28e83SPiotr Jasiukajtis 
4825c28e83SPiotr Jasiukajtis /* array indices used to access words within a double */
4925c28e83SPiotr Jasiukajtis #define	HIWORD	1
5025c28e83SPiotr Jasiukajtis #define	LOWORD	0
5125c28e83SPiotr Jasiukajtis 
5225c28e83SPiotr Jasiukajtis /* structure used to access words within a quad */
5325c28e83SPiotr Jasiukajtis union longdouble {
5425c28e83SPiotr Jasiukajtis 	struct {
5525c28e83SPiotr Jasiukajtis 		unsigned int	frac4;
5625c28e83SPiotr Jasiukajtis 		unsigned int	frac3;
5725c28e83SPiotr Jasiukajtis 		unsigned int	frac2;
5825c28e83SPiotr Jasiukajtis 		unsigned int	msw;
5925c28e83SPiotr Jasiukajtis 	} l;
6025c28e83SPiotr Jasiukajtis 	long double	d;
6125c28e83SPiotr Jasiukajtis };
6225c28e83SPiotr Jasiukajtis 
6325c28e83SPiotr Jasiukajtis /* default NaN returned for sqrt(neg) */
6425c28e83SPiotr Jasiukajtis static const union longdouble
6525c28e83SPiotr Jasiukajtis 	qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
6625c28e83SPiotr Jasiukajtis 
6725c28e83SPiotr Jasiukajtis /* signalling NaN used to raise invalid */
6825c28e83SPiotr Jasiukajtis static const union {
6925c28e83SPiotr Jasiukajtis 	unsigned u[2];
7025c28e83SPiotr Jasiukajtis 	double d;
7125c28e83SPiotr Jasiukajtis } snan = { 0, 0x7ff00001 };
7225c28e83SPiotr Jasiukajtis 
7325c28e83SPiotr Jasiukajtis #else
7425c28e83SPiotr Jasiukajtis 
7525c28e83SPiotr Jasiukajtis /* array indices used to access words within a double */
7625c28e83SPiotr Jasiukajtis #define	HIWORD	0
7725c28e83SPiotr Jasiukajtis #define	LOWORD	1
7825c28e83SPiotr Jasiukajtis 
7925c28e83SPiotr Jasiukajtis /* structure used to access words within a quad */
8025c28e83SPiotr Jasiukajtis union longdouble {
8125c28e83SPiotr Jasiukajtis 	struct {
8225c28e83SPiotr Jasiukajtis 		unsigned int	msw;
8325c28e83SPiotr Jasiukajtis 		unsigned int	frac2;
8425c28e83SPiotr Jasiukajtis 		unsigned int	frac3;
8525c28e83SPiotr Jasiukajtis 		unsigned int	frac4;
8625c28e83SPiotr Jasiukajtis 	} l;
8725c28e83SPiotr Jasiukajtis 	long double	d;
8825c28e83SPiotr Jasiukajtis };
8925c28e83SPiotr Jasiukajtis 
9025c28e83SPiotr Jasiukajtis /* default NaN returned for sqrt(neg) */
9125c28e83SPiotr Jasiukajtis static const union longdouble
9225c28e83SPiotr Jasiukajtis 	qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
9325c28e83SPiotr Jasiukajtis 
9425c28e83SPiotr Jasiukajtis /* signalling NaN used to raise invalid */
9525c28e83SPiotr Jasiukajtis static const union {
9625c28e83SPiotr Jasiukajtis 	unsigned u[2];
9725c28e83SPiotr Jasiukajtis 	double d;
9825c28e83SPiotr Jasiukajtis } snan = { 0x7ff00001, 0 };
9925c28e83SPiotr Jasiukajtis 
10025c28e83SPiotr Jasiukajtis #endif /* __LITTLE_ENDIAN */
10125c28e83SPiotr Jasiukajtis 
10225c28e83SPiotr Jasiukajtis 
10325c28e83SPiotr Jasiukajtis static const double
10425c28e83SPiotr Jasiukajtis 	zero = 0.0,
10525c28e83SPiotr Jasiukajtis 	half = 0.5,
10625c28e83SPiotr Jasiukajtis 	one = 1.0,
10725c28e83SPiotr Jasiukajtis 	huge = 1.0e300,
10825c28e83SPiotr Jasiukajtis 	tiny = 1.0e-300,
10925c28e83SPiotr Jasiukajtis 	two36 = 6.87194767360000000000e+10,
11025c28e83SPiotr Jasiukajtis 	two30 = 1.07374182400000000000e+09,
11125c28e83SPiotr Jasiukajtis 	two6 = 6.40000000000000000000e+01,
11225c28e83SPiotr Jasiukajtis 	two4 = 1.60000000000000000000e+01,
11325c28e83SPiotr Jasiukajtis 	twom18 = 3.81469726562500000000e-06,
11425c28e83SPiotr Jasiukajtis 	twom28 = 3.72529029846191406250e-09,
11525c28e83SPiotr Jasiukajtis 	twom42 = 2.27373675443232059479e-13,
11625c28e83SPiotr Jasiukajtis 	twom60 = 8.67361737988403547206e-19,
11725c28e83SPiotr Jasiukajtis 	twom62 = 2.16840434497100886801e-19,
11825c28e83SPiotr Jasiukajtis 	twom66 = 1.35525271560688054251e-20,
11925c28e83SPiotr Jasiukajtis 	twom90 = 8.07793566946316088742e-28,
12025c28e83SPiotr Jasiukajtis 	twom113 = 9.62964972193617926528e-35,
12125c28e83SPiotr Jasiukajtis 	twom124 = 4.70197740328915003187e-38;
12225c28e83SPiotr Jasiukajtis 
12325c28e83SPiotr Jasiukajtis 
12425c28e83SPiotr Jasiukajtis /*
12525c28e83SPiotr Jasiukajtis *	Extract the exponent and normalized significand (represented as
12625c28e83SPiotr Jasiukajtis *	an array of five doubles) from a finite, nonzero quad.
12725c28e83SPiotr Jasiukajtis */
12825c28e83SPiotr Jasiukajtis static int
__q_unpack(const union longdouble * x,double * s)12925c28e83SPiotr Jasiukajtis __q_unpack(const union longdouble *x, double *s)
13025c28e83SPiotr Jasiukajtis {
13125c28e83SPiotr Jasiukajtis 	union {
13225c28e83SPiotr Jasiukajtis 		double			d;
13325c28e83SPiotr Jasiukajtis 		unsigned int	l[2];
13425c28e83SPiotr Jasiukajtis 	} u;
13525c28e83SPiotr Jasiukajtis 	double			b;
13625c28e83SPiotr Jasiukajtis 	unsigned int	lx, w[3];
13725c28e83SPiotr Jasiukajtis 	int				ex;
13825c28e83SPiotr Jasiukajtis 
13925c28e83SPiotr Jasiukajtis 	/* get the normalized significand and exponent */
14025c28e83SPiotr Jasiukajtis 	ex = (int) ((x->l.msw & 0x7fffffff) >> 16);
14125c28e83SPiotr Jasiukajtis 	lx = x->l.msw & 0xffff;
14225c28e83SPiotr Jasiukajtis 	if (ex)
14325c28e83SPiotr Jasiukajtis 	{
14425c28e83SPiotr Jasiukajtis 		lx |= 0x10000;
14525c28e83SPiotr Jasiukajtis 		w[0] = x->l.frac2;
14625c28e83SPiotr Jasiukajtis 		w[1] = x->l.frac3;
14725c28e83SPiotr Jasiukajtis 		w[2] = x->l.frac4;
14825c28e83SPiotr Jasiukajtis 	}
14925c28e83SPiotr Jasiukajtis 	else
15025c28e83SPiotr Jasiukajtis 	{
15125c28e83SPiotr Jasiukajtis 		if (lx | (x->l.frac2 & 0xfffe0000))
15225c28e83SPiotr Jasiukajtis 		{
15325c28e83SPiotr Jasiukajtis 			w[0] = x->l.frac2;
15425c28e83SPiotr Jasiukajtis 			w[1] = x->l.frac3;
15525c28e83SPiotr Jasiukajtis 			w[2] = x->l.frac4;
15625c28e83SPiotr Jasiukajtis 			ex = 1;
15725c28e83SPiotr Jasiukajtis 		}
15825c28e83SPiotr Jasiukajtis 		else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000))
15925c28e83SPiotr Jasiukajtis 		{
16025c28e83SPiotr Jasiukajtis 			lx = x->l.frac2;
16125c28e83SPiotr Jasiukajtis 			w[0] = x->l.frac3;
16225c28e83SPiotr Jasiukajtis 			w[1] = x->l.frac4;
16325c28e83SPiotr Jasiukajtis 			w[2] = 0;
16425c28e83SPiotr Jasiukajtis 			ex = -31;
16525c28e83SPiotr Jasiukajtis 		}
16625c28e83SPiotr Jasiukajtis 		else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000))
16725c28e83SPiotr Jasiukajtis 		{
16825c28e83SPiotr Jasiukajtis 			lx = x->l.frac3;
16925c28e83SPiotr Jasiukajtis 			w[0] = x->l.frac4;
17025c28e83SPiotr Jasiukajtis 			w[1] = w[2] = 0;
17125c28e83SPiotr Jasiukajtis 			ex = -63;
17225c28e83SPiotr Jasiukajtis 		}
17325c28e83SPiotr Jasiukajtis 		else
17425c28e83SPiotr Jasiukajtis 		{
17525c28e83SPiotr Jasiukajtis 			lx = x->l.frac4;
17625c28e83SPiotr Jasiukajtis 			w[0] = w[1] = w[2] = 0;
17725c28e83SPiotr Jasiukajtis 			ex = -95;
17825c28e83SPiotr Jasiukajtis 		}
17925c28e83SPiotr Jasiukajtis 		while ((lx & 0x10000) == 0)
18025c28e83SPiotr Jasiukajtis 		{
18125c28e83SPiotr Jasiukajtis 			lx = (lx << 1) | (w[0] >> 31);
18225c28e83SPiotr Jasiukajtis 			w[0] = (w[0] << 1) | (w[1] >> 31);
18325c28e83SPiotr Jasiukajtis 			w[1] = (w[1] << 1) | (w[2] >> 31);
18425c28e83SPiotr Jasiukajtis 			w[2] <<= 1;
18525c28e83SPiotr Jasiukajtis 			ex--;
18625c28e83SPiotr Jasiukajtis 		}
18725c28e83SPiotr Jasiukajtis 	}
18825c28e83SPiotr Jasiukajtis 
18925c28e83SPiotr Jasiukajtis 	/* extract the significand into five doubles */
19025c28e83SPiotr Jasiukajtis 	u.l[HIWORD] = 0x42300000;
19125c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = 0;
19225c28e83SPiotr Jasiukajtis 	b = u.d;
19325c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = lx;
19425c28e83SPiotr Jasiukajtis 	s[0] = u.d - b;
19525c28e83SPiotr Jasiukajtis 
19625c28e83SPiotr Jasiukajtis 	u.l[HIWORD] = 0x40300000;
19725c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = 0;
19825c28e83SPiotr Jasiukajtis 	b = u.d;
19925c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = w[0] & 0xffffff00;
20025c28e83SPiotr Jasiukajtis 	s[1] = u.d - b;
20125c28e83SPiotr Jasiukajtis 
20225c28e83SPiotr Jasiukajtis 	u.l[HIWORD] = 0x3e300000;
20325c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = 0;
20425c28e83SPiotr Jasiukajtis 	b = u.d;
20525c28e83SPiotr Jasiukajtis 	u.l[HIWORD] |= w[0] & 0xff;
20625c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = w[1] & 0xffff0000;
20725c28e83SPiotr Jasiukajtis 	s[2] = u.d - b;
20825c28e83SPiotr Jasiukajtis 
20925c28e83SPiotr Jasiukajtis 	u.l[HIWORD] = 0x3c300000;
21025c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = 0;
21125c28e83SPiotr Jasiukajtis 	b = u.d;
21225c28e83SPiotr Jasiukajtis 	u.l[HIWORD] |= w[1] & 0xffff;
21325c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = w[2] & 0xff000000;
21425c28e83SPiotr Jasiukajtis 	s[3] = u.d - b;
21525c28e83SPiotr Jasiukajtis 
21625c28e83SPiotr Jasiukajtis 	u.l[HIWORD] = 0x3c300000;
21725c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = 0;
21825c28e83SPiotr Jasiukajtis 	b = u.d;
21925c28e83SPiotr Jasiukajtis 	u.l[LOWORD] = w[2] & 0xffffff;
22025c28e83SPiotr Jasiukajtis 	s[4] = u.d - b;
22125c28e83SPiotr Jasiukajtis 
22225c28e83SPiotr Jasiukajtis 	return ex - 0x3fff;
22325c28e83SPiotr Jasiukajtis }
22425c28e83SPiotr Jasiukajtis 
22525c28e83SPiotr Jasiukajtis 
22625c28e83SPiotr Jasiukajtis /*
22725c28e83SPiotr Jasiukajtis *	Pack an exponent and array of three doubles representing a finite,
22825c28e83SPiotr Jasiukajtis *	nonzero number into a quad.  Assume the sign is already there and
22925c28e83SPiotr Jasiukajtis *	the rounding mode has been fudged accordingly.
23025c28e83SPiotr Jasiukajtis */
23125c28e83SPiotr Jasiukajtis static void
__q_pack(const double * z,int exp,enum fp_direction_type rm,union longdouble * x,int * inexact)23225c28e83SPiotr Jasiukajtis __q_pack(const double *z, int exp, enum fp_direction_type rm,
23325c28e83SPiotr Jasiukajtis 	union longdouble *x, int *inexact)
23425c28e83SPiotr Jasiukajtis {
23525c28e83SPiotr Jasiukajtis 	union {
23625c28e83SPiotr Jasiukajtis 		double			d;
23725c28e83SPiotr Jasiukajtis 		unsigned int	l[2];
23825c28e83SPiotr Jasiukajtis 	} u;
23925c28e83SPiotr Jasiukajtis 	double			s[3], t, t2;
24025c28e83SPiotr Jasiukajtis 	unsigned int	msw, frac2, frac3, frac4;
24125c28e83SPiotr Jasiukajtis 
24225c28e83SPiotr Jasiukajtis 	/* bias exponent and strip off integer bit */
24325c28e83SPiotr Jasiukajtis 	exp += 0x3fff;
24425c28e83SPiotr Jasiukajtis 	s[0] = z[0] - one;
24525c28e83SPiotr Jasiukajtis 	s[1] = z[1];
24625c28e83SPiotr Jasiukajtis 	s[2] = z[2];
24725c28e83SPiotr Jasiukajtis 
24825c28e83SPiotr Jasiukajtis 	/*
24925c28e83SPiotr Jasiukajtis 	 * chop the significand to obtain the fraction;
25025c28e83SPiotr Jasiukajtis 	 * use round-to-minus-infinity to ensure chopping
25125c28e83SPiotr Jasiukajtis 	 */
25225c28e83SPiotr Jasiukajtis 	(void) __swapRD(fp_negative);
25325c28e83SPiotr Jasiukajtis 
25425c28e83SPiotr Jasiukajtis 	/* extract the first eighty bits of fraction */
25525c28e83SPiotr Jasiukajtis 	t = s[1] + s[2];
25625c28e83SPiotr Jasiukajtis 	u.d = two36 + (s[0] + t);
25725c28e83SPiotr Jasiukajtis 	msw = u.l[LOWORD];
25825c28e83SPiotr Jasiukajtis 	s[0] -= (u.d - two36);
25925c28e83SPiotr Jasiukajtis 
26025c28e83SPiotr Jasiukajtis 	u.d = two4 + (s[0] + t);
26125c28e83SPiotr Jasiukajtis 	frac2 = u.l[LOWORD];
26225c28e83SPiotr Jasiukajtis 	s[0] -= (u.d - two4);
26325c28e83SPiotr Jasiukajtis 
26425c28e83SPiotr Jasiukajtis 	u.d = twom28 + (s[0] + t);
26525c28e83SPiotr Jasiukajtis 	frac3 = u.l[LOWORD];
26625c28e83SPiotr Jasiukajtis 	s[0] -= (u.d - twom28);
26725c28e83SPiotr Jasiukajtis 
26825c28e83SPiotr Jasiukajtis 	/* condense the remaining fraction; errors here won't matter */
26925c28e83SPiotr Jasiukajtis 	t = s[0] + s[1];
27025c28e83SPiotr Jasiukajtis 	s[1] = ((s[0] - t) + s[1]) + s[2];
27125c28e83SPiotr Jasiukajtis 	s[0] = t;
27225c28e83SPiotr Jasiukajtis 
27325c28e83SPiotr Jasiukajtis 	/* get the last word of fraction */
27425c28e83SPiotr Jasiukajtis 	u.d = twom60 + (s[0] + s[1]);
27525c28e83SPiotr Jasiukajtis 	frac4 = u.l[LOWORD];
27625c28e83SPiotr Jasiukajtis 	s[0] -= (u.d - twom60);
27725c28e83SPiotr Jasiukajtis 
27825c28e83SPiotr Jasiukajtis 	/*
27925c28e83SPiotr Jasiukajtis 	 * keep track of what's left for rounding; note that
28025c28e83SPiotr Jasiukajtis 	 * t2 will be non-negative due to rounding mode
28125c28e83SPiotr Jasiukajtis 	 */
28225c28e83SPiotr Jasiukajtis 	t = s[0] + s[1];
28325c28e83SPiotr Jasiukajtis 	t2 = (s[0] - t) + s[1];
28425c28e83SPiotr Jasiukajtis 
28525c28e83SPiotr Jasiukajtis 	if (t != zero)
28625c28e83SPiotr Jasiukajtis 	{
28725c28e83SPiotr Jasiukajtis 		*inexact = 1;
28825c28e83SPiotr Jasiukajtis 
28925c28e83SPiotr Jasiukajtis 		/* decide whether to round the fraction up */
29025c28e83SPiotr Jasiukajtis 		if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
29125c28e83SPiotr Jasiukajtis 			(t == twom113 && (t2 != zero || frac4 & 1)))))
29225c28e83SPiotr Jasiukajtis 		{
29325c28e83SPiotr Jasiukajtis 			/* round up and renormalize if necessary */
29425c28e83SPiotr Jasiukajtis 			if (++frac4 == 0)
29525c28e83SPiotr Jasiukajtis 				if (++frac3 == 0)
29625c28e83SPiotr Jasiukajtis 					if (++frac2 == 0)
29725c28e83SPiotr Jasiukajtis 						if (++msw == 0x10000)
29825c28e83SPiotr Jasiukajtis 						{
29925c28e83SPiotr Jasiukajtis 							msw = 0;
30025c28e83SPiotr Jasiukajtis 							exp++;
30125c28e83SPiotr Jasiukajtis 						}
30225c28e83SPiotr Jasiukajtis 		}
30325c28e83SPiotr Jasiukajtis 	}
30425c28e83SPiotr Jasiukajtis 
30525c28e83SPiotr Jasiukajtis 	/* assemble the result */
30625c28e83SPiotr Jasiukajtis 	x->l.msw |= msw | (exp << 16);
30725c28e83SPiotr Jasiukajtis 	x->l.frac2 = frac2;
30825c28e83SPiotr Jasiukajtis 	x->l.frac3 = frac3;
30925c28e83SPiotr Jasiukajtis 	x->l.frac4 = frac4;
31025c28e83SPiotr Jasiukajtis }
31125c28e83SPiotr Jasiukajtis 
31225c28e83SPiotr Jasiukajtis 
31325c28e83SPiotr Jasiukajtis /*
31425c28e83SPiotr Jasiukajtis *	Compute the square root of x and place the TP result in s.
31525c28e83SPiotr Jasiukajtis */
31625c28e83SPiotr Jasiukajtis static void
__q_tp_sqrt(const double * x,double * s)31725c28e83SPiotr Jasiukajtis __q_tp_sqrt(const double *x, double *s)
31825c28e83SPiotr Jasiukajtis {
31925c28e83SPiotr Jasiukajtis 	double	c, rr, r[3], tt[3], t[5];
32025c28e83SPiotr Jasiukajtis 
32125c28e83SPiotr Jasiukajtis 	/* approximate the divisor for the Newton iteration */
32225c28e83SPiotr Jasiukajtis 	c = sqrt((x[0] + x[1]) + x[2]);
32325c28e83SPiotr Jasiukajtis 	rr = half / c;
32425c28e83SPiotr Jasiukajtis 
32525c28e83SPiotr Jasiukajtis 	/* compute the first five "digits" of the square root */
32625c28e83SPiotr Jasiukajtis 	t[0] = (c + two30) - two30;
32725c28e83SPiotr Jasiukajtis 	tt[0] = t[0] + t[0];
32825c28e83SPiotr Jasiukajtis 	r[0] = ((x[0] - t[0] * t[0]) + x[1]) + x[2];
32925c28e83SPiotr Jasiukajtis 
33025c28e83SPiotr Jasiukajtis 	t[1] = (rr * (r[0] + x[3]) + two6) - two6;
33125c28e83SPiotr Jasiukajtis 	tt[1] = t[1] + t[1];
33225c28e83SPiotr Jasiukajtis 	r[0] -= tt[0] * t[1];
33325c28e83SPiotr Jasiukajtis 	r[1] = x[3] - t[1] * t[1];
33425c28e83SPiotr Jasiukajtis 	c = (r[1] + twom18) - twom18;
33525c28e83SPiotr Jasiukajtis 	r[0] += c;
33625c28e83SPiotr Jasiukajtis 	r[1] = (r[1] - c) + x[4];
33725c28e83SPiotr Jasiukajtis 
33825c28e83SPiotr Jasiukajtis 	t[2] = (rr * (r[0] + r[1]) + twom18) - twom18;
33925c28e83SPiotr Jasiukajtis 	tt[2] = t[2] + t[2];
34025c28e83SPiotr Jasiukajtis 	r[0] -= tt[0] * t[2];
34125c28e83SPiotr Jasiukajtis 	r[1] -= tt[1] * t[2];
34225c28e83SPiotr Jasiukajtis 	c = (r[1] + twom42) - twom42;
34325c28e83SPiotr Jasiukajtis 	r[0] += c;
34425c28e83SPiotr Jasiukajtis 	r[1] = (r[1] - c) - t[2] * t[2];
34525c28e83SPiotr Jasiukajtis 
34625c28e83SPiotr Jasiukajtis 	t[3] = (rr * (r[0] + r[1]) + twom42) - twom42;
34725c28e83SPiotr Jasiukajtis 	r[0] = ((r[0] - tt[0] * t[3]) + r[1]) - tt[1] * t[3];
34825c28e83SPiotr Jasiukajtis 	r[1] = -tt[2] * t[3];
34925c28e83SPiotr Jasiukajtis 	c = (r[1] + twom90) - twom90;
35025c28e83SPiotr Jasiukajtis 	r[0] += c;
35125c28e83SPiotr Jasiukajtis 	r[1] = (r[1] - c) - t[3] * t[3];
35225c28e83SPiotr Jasiukajtis 
35325c28e83SPiotr Jasiukajtis 	t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
35425c28e83SPiotr Jasiukajtis 
35525c28e83SPiotr Jasiukajtis 	/* here we just need to get the sign of the remainder */
35625c28e83SPiotr Jasiukajtis 	c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1])
35725c28e83SPiotr Jasiukajtis 		- tt[2] * t[4]) - (t[3] + t[3]) * t[4]) - t[4] * t[4];
35825c28e83SPiotr Jasiukajtis 
35925c28e83SPiotr Jasiukajtis 	/* reduce to three doubles */
36025c28e83SPiotr Jasiukajtis 	t[0] += t[1];
36125c28e83SPiotr Jasiukajtis 	t[1] = t[2] + t[3];
36225c28e83SPiotr Jasiukajtis 	t[2] = t[4];
36325c28e83SPiotr Jasiukajtis 
36425c28e83SPiotr Jasiukajtis 	/* if the third term might lie on a rounding boundary, perturb it */
36525c28e83SPiotr Jasiukajtis 	if (c != zero && t[2] == (twom62 + t[2]) - twom62)
36625c28e83SPiotr Jasiukajtis 	{
36725c28e83SPiotr Jasiukajtis 		if (c < zero)
36825c28e83SPiotr Jasiukajtis 			t[2] -= twom124;
36925c28e83SPiotr Jasiukajtis 		else
37025c28e83SPiotr Jasiukajtis 			t[2] += twom124;
37125c28e83SPiotr Jasiukajtis 	}
37225c28e83SPiotr Jasiukajtis 
37325c28e83SPiotr Jasiukajtis 	/* condense the square root */
37425c28e83SPiotr Jasiukajtis 	c = t[1] + t[2];
37525c28e83SPiotr Jasiukajtis 	t[2] += (t[1] - c);
37625c28e83SPiotr Jasiukajtis 	t[1] = c;
37725c28e83SPiotr Jasiukajtis 	c = t[0] + t[1];
37825c28e83SPiotr Jasiukajtis 	s[1] = t[1] + (t[0] - c);
37925c28e83SPiotr Jasiukajtis 	s[0] = c;
38025c28e83SPiotr Jasiukajtis 	if (s[1] == zero)
38125c28e83SPiotr Jasiukajtis 	{
38225c28e83SPiotr Jasiukajtis 		c = s[0] + t[2];
38325c28e83SPiotr Jasiukajtis 		s[1] = t[2] + (s[0] - c);
38425c28e83SPiotr Jasiukajtis 		s[0] = c;
38525c28e83SPiotr Jasiukajtis 		s[2] = zero;
38625c28e83SPiotr Jasiukajtis 	}
38725c28e83SPiotr Jasiukajtis 	else
38825c28e83SPiotr Jasiukajtis 	{
38925c28e83SPiotr Jasiukajtis 		c = s[1] + t[2];
39025c28e83SPiotr Jasiukajtis 		s[2] = t[2] + (s[1] - c);
39125c28e83SPiotr Jasiukajtis 		s[1] = c;
39225c28e83SPiotr Jasiukajtis 	}
39325c28e83SPiotr Jasiukajtis }
39425c28e83SPiotr Jasiukajtis 
39525c28e83SPiotr Jasiukajtis 
39625c28e83SPiotr Jasiukajtis long double
sqrtl(long double ldx)39725c28e83SPiotr Jasiukajtis sqrtl(long double ldx)
39825c28e83SPiotr Jasiukajtis {
39925c28e83SPiotr Jasiukajtis 	union	longdouble		x;
40025c28e83SPiotr Jasiukajtis 	volatile double			t;
40125c28e83SPiotr Jasiukajtis 	double					xx[5], zz[3];
40225c28e83SPiotr Jasiukajtis 	enum fp_direction_type	rm;
40325c28e83SPiotr Jasiukajtis 	int				ex, inexact, exc, traps;
40425c28e83SPiotr Jasiukajtis 
40525c28e83SPiotr Jasiukajtis 	/* clear cexc */
40625c28e83SPiotr Jasiukajtis 	t = zero;
40725c28e83SPiotr Jasiukajtis 	t -= zero;
40825c28e83SPiotr Jasiukajtis 
40925c28e83SPiotr Jasiukajtis 	/* check for zero operand */
41025c28e83SPiotr Jasiukajtis 	x.d = ldx;
41125c28e83SPiotr Jasiukajtis 	if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
41225c28e83SPiotr Jasiukajtis 		return ldx;
41325c28e83SPiotr Jasiukajtis 
41425c28e83SPiotr Jasiukajtis 	/* handle nan and inf cases */
41525c28e83SPiotr Jasiukajtis 	if ((x.l.msw & 0x7fffffff) >= 0x7fff0000)
41625c28e83SPiotr Jasiukajtis 	{
41725c28e83SPiotr Jasiukajtis 		if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4)
41825c28e83SPiotr Jasiukajtis 		{
41925c28e83SPiotr Jasiukajtis 			if (!(x.l.msw & 0x8000))
42025c28e83SPiotr Jasiukajtis 			{
42125c28e83SPiotr Jasiukajtis 				/* snan, signal invalid */
42225c28e83SPiotr Jasiukajtis 				t += snan.d;
42325c28e83SPiotr Jasiukajtis 			}
42425c28e83SPiotr Jasiukajtis 			x.l.msw |= 0x8000;
42525c28e83SPiotr Jasiukajtis 			return x.d;
42625c28e83SPiotr Jasiukajtis 		}
42725c28e83SPiotr Jasiukajtis 		if (x.l.msw & 0x80000000)
42825c28e83SPiotr Jasiukajtis 		{
42925c28e83SPiotr Jasiukajtis 			/* sqrt(-inf), signal invalid */
43025c28e83SPiotr Jasiukajtis 			t = -one;
43125c28e83SPiotr Jasiukajtis 			t = sqrt(t);
43225c28e83SPiotr Jasiukajtis 			return qnan.d;
43325c28e83SPiotr Jasiukajtis 		}
43425c28e83SPiotr Jasiukajtis 		/* sqrt(inf), return inf */
43525c28e83SPiotr Jasiukajtis 		return x.d;
43625c28e83SPiotr Jasiukajtis 	}
43725c28e83SPiotr Jasiukajtis 
43825c28e83SPiotr Jasiukajtis 	/* handle negative numbers */
43925c28e83SPiotr Jasiukajtis 	if (x.l.msw & 0x80000000)
44025c28e83SPiotr Jasiukajtis 	{
44125c28e83SPiotr Jasiukajtis 		t = -one;
44225c28e83SPiotr Jasiukajtis 		t = sqrt(t);
44325c28e83SPiotr Jasiukajtis 		return qnan.d;
44425c28e83SPiotr Jasiukajtis 	}
44525c28e83SPiotr Jasiukajtis 
44625c28e83SPiotr Jasiukajtis 	/* now x is finite, positive */
44725c28e83SPiotr Jasiukajtis 
44825c28e83SPiotr Jasiukajtis 	traps = __swapTE(0);
44925c28e83SPiotr Jasiukajtis 	exc = __swapEX(0);
45025c28e83SPiotr Jasiukajtis 	rm = __swapRD(fp_nearest);
45125c28e83SPiotr Jasiukajtis 
45225c28e83SPiotr Jasiukajtis 	ex = __q_unpack(&x, xx);
45325c28e83SPiotr Jasiukajtis 	if (ex & 1)
45425c28e83SPiotr Jasiukajtis 	{
45525c28e83SPiotr Jasiukajtis 		/* make exponent even */
45625c28e83SPiotr Jasiukajtis 		xx[0] += xx[0];
45725c28e83SPiotr Jasiukajtis 		xx[1] += xx[1];
45825c28e83SPiotr Jasiukajtis 		xx[2] += xx[2];
45925c28e83SPiotr Jasiukajtis 		xx[3] += xx[3];
46025c28e83SPiotr Jasiukajtis 		xx[4] += xx[4];
46125c28e83SPiotr Jasiukajtis 		ex--;
46225c28e83SPiotr Jasiukajtis 	}
46325c28e83SPiotr Jasiukajtis 	__q_tp_sqrt(xx, zz);
46425c28e83SPiotr Jasiukajtis 
46525c28e83SPiotr Jasiukajtis 	/* put everything together */
46625c28e83SPiotr Jasiukajtis 	x.l.msw = 0;
46725c28e83SPiotr Jasiukajtis 	inexact = 0;
46825c28e83SPiotr Jasiukajtis 	__q_pack(zz, ex >> 1, rm, &x, &inexact);
46925c28e83SPiotr Jasiukajtis 
47025c28e83SPiotr Jasiukajtis 	(void) __swapRD(rm);
47125c28e83SPiotr Jasiukajtis 	(void) __swapEX(exc);
47225c28e83SPiotr Jasiukajtis 	(void) __swapTE(traps);
47325c28e83SPiotr Jasiukajtis 	if (inexact)
47425c28e83SPiotr Jasiukajtis 	{
47525c28e83SPiotr Jasiukajtis 		t = huge;
47625c28e83SPiotr Jasiukajtis 		t += tiny;
47725c28e83SPiotr Jasiukajtis 	}
47825c28e83SPiotr Jasiukajtis 	return x.d;
47925c28e83SPiotr Jasiukajtis }
480