From bcd797a5ba4631c031919dad832d670e564212e9 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Tue, 3 Sep 2013 15:02:10 +0000
Subject: math: long double inverse trigonometric cleanup (acosl, asinl, atanl,
 atan2l)

* added ld128 support from freebsd fdlibm (untested)
* using new ldshape union instead of IEEEl2bits
* inexact status flag is not supported
---
 src/math/asinl.c | 44 ++++++++++++++++++++++----------------------
 1 file changed, 22 insertions(+), 22 deletions(-)

(limited to 'src/math/asinl.c')

diff --git a/src/math/asinl.c b/src/math/asinl.c
index 8799341d..347c5356 100644
--- a/src/math/asinl.c
+++ b/src/math/asinl.c
@@ -23,27 +23,29 @@ long double asinl(long double x)
 }
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 #include "__invtrigl.h"
-/* 0.95 */
-#define THRESH  ((0xe666666666666666ULL>>(64-(LDBL_MANH_SIZE-1)))|LDBL_NBIT)
+#if LDBL_MANT_DIG == 64
+#define CLOSETO1(u) (u.i.m>>56 >= 0xf7)
+#define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32)
+#elif LDBL_MANT_DIG == 113
+#define CLOSETO1(u) (u.i.top >= 0xee00)
+#define CLEARBOTTOM(u) (u.i.lo = 0)
+#endif
 
 long double asinl(long double x)
 {
-	union IEEEl2bits u;
-	long double z,r,s;
-	uint16_t expsign, expt;
+	union ldshape u = {x};
+	long double z, r, s;
+	uint16_t e = u.i.se & 0x7fff;
+	int sign = u.i.se >> 15;
 
-	u.e = x;
-	expsign = u.xbits.expsign;
-	expt = expsign & 0x7fff;
-	if (expt >= 0x3fff) {   /* |x| >= 1 or nan */
-		if (expt == 0x3fff &&
-		    ((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
-			/* asin(+-1)=+-pi/2 with inexact */
+	if (e >= 0x3fff) {   /* |x| >= 1 or nan */
+		/* asin(+-1)=+-pi/2 with inexact */
+		if (x == 1 || x == -1)
 			return x*pio2_hi + 0x1p-120f;
 		return 0/(x-x);
 	}
-	if (expt < 0x3fff - 1) {  /* |x| < 0.5 */
-		if (expt < 0x3fff - 32) {  /* |x|<0x1p-32, asinl(x)=x */
+	if (e < 0x3fff - 1) {  /* |x| < 0.5 */
+		if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) {
 			/* return x with inexact if x!=0 */
 			FORCE_EVAL(x + 0x1p120f);
 			return x;
@@ -54,18 +56,16 @@ long double asinl(long double x)
 	z = (1.0 - fabsl(x))*0.5;
 	s = sqrtl(z);
 	r = __invtrigl_R(z);
-	if (u.bits.manh >= THRESH) { /* if |x| is close to 1 */
+	if (CLOSETO1(u)) {
 		x = pio2_hi - (2*(s+s*r)-pio2_lo);
 	} else {
 		long double f, c;
-		u.e = s;
-		u.bits.manl = 0;
-		f = u.e;
-		c = (z-f*f)/(s+f);
+		u.f = s;
+		CLEARBOTTOM(u);
+		f = u.f;
+		c = (z - f*f)/(s + f);
 		x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
 	}
-	if (expsign>>15)
-		return -x;
-	return x;
+	return sign ? -x : x;
 }
 #endif
-- 
cgit v1.2.1