math: long double trigonometric cleanup (cosl, sinl, sincosl, tanl)

ld128 support was added to internal kernel functions (__cosl, __sinl,
__tanl, __rem_pio2l) from freebsd (not tested, but should be a good
start for when ld128 arch arrives)

__rem_pio2l had some code cleanup, the freebsd ld128 code seems to
gather the results of a large reduction with precision loss (fixed
the bug but a todo comment was added for later investigation)

the old copyright was removed from the non-kernel wrapper functions
(cosl, sinl, sincosl, tanl) since these are trivial and the interesting
parts and comments had been already rewritten.

src/math/__cosl.c

@@ -1,4 +1,5 @@
 /* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */
+/* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -14,7 +15,8 @@
 
 #include "libm.h"
 
-#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+#if LDBL_MANT_DIG == 64
 /*
  * ld80 version of __cos.c.  See __cos.c for most comments.
  */
@@ -43,7 +45,6 @@
  */
 static const long double
 C1 =  0.0416666666666666666136L;        /*  0xaaaaaaaaaaaaaa9b.0p-68 */
-
 static const double
 C2 = -0.0013888888888888874,            /* -0x16c16c16c16c10.0p-62 */
 C3 =  0.000024801587301571716,          /*  0x1a01a01a018e22.0p-68 */
@@ -51,13 +52,43 @@ C4 = -0.00000027557319215507120,        /* -0x127e4fb7602f22.0p-74 */
 C5 =  0.0000000020876754400407278,      /*  0x11eed8caaeccf1.0p-81 */
 C6 = -1.1470297442401303e-11,           /* -0x19393412bd1529.0p-89 */
 C7 =  4.7383039476436467e-14;           /*  0x1aac9d9af5c43e.0p-97 */
+#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7)))))))
+#elif LDBL_MANT_DIG == 113
+/*
+ * ld128 version of __cos.c.  See __cos.c for most comments.
+ */
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]:
+ * |cos(x) - c(x))| < 2**-122.0
+ *
+ * 113-bit precision requires more care than 64-bit precision, since
+ * simple methods give a minimax polynomial with coefficient for x^2
+ * that is 1 ulp below 0.5, but we want it to be precisely 0.5.  See
+ * above for more details.
+ */
+static const long double
+C1 =  0.04166666666666666666666666666666658424671L,
+C2 = -0.001388888888888888888888888888863490893732L,
+C3 =  0.00002480158730158730158730158600795304914210L,
+C4 = -0.2755731922398589065255474947078934284324e-6L,
+C5 =  0.2087675698786809897659225313136400793948e-8L,
+C6 = -0.1147074559772972315817149986812031204775e-10L,
+C7 =  0.4779477332386808976875457937252120293400e-13L;
+static const double
+C8 = -0.1561920696721507929516718307820958119868e-15,
+C9 =  0.4110317413744594971475941557607804508039e-18,
+C10 = -0.8896592467191938803288521958313920156409e-21,
+C11 =  0.1601061435794535138244346256065192782581e-23;
+#define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+ \
+	z*(C8+z*(C9+z*(C10+z*C11)))))))))))
+#endif
 
 long double __cosl(long double x, long double y)
 {
 	long double hz,z,r,w;
 
 	z  = x*x;
-	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))));
+	r  = POLY(z);
 	hz = 0.5*z;
 	w  = 1.0-hz;
 	return w + (((1.0-w)-hz) + (z*r-x*y));

src/math/__rem_pio2l.c

@@ -13,15 +13,22 @@
  * Optimized by Bruce D. Evans.
  */
 #include "libm.h"
-#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-/* ld80 version of __rem_pio2(x,y)
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+/* ld80 and ld128 version of __rem_pio2(x,y)
  *
  * return the remainder of x rem pi/2 in y[0]+y[1]
  * use __rem_pio2_large() for large x
  */
 
-#define BIAS    (LDBL_MAX_EXP - 1)
-
+#if LDBL_MANT_DIG == 64
+/* u ~< 0x1p25*pi/2 */
+#define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.m>>48) < ((0x3fff + 25)<<16 | 0x921f>>1 | 0x8000))
+#define TOINT 0x1.8p63
+#define QUOBITS(x) ((uint32_t)(int32_t)x & 0x7fffffff)
+#define ROUND1 22
+#define ROUND2 61
+#define NX 3
+#define NY 2
 /*
  * invpio2:  64 bits of 2/pi
  * pio2_1:   first  39 bits of pi/2
@@ -32,60 +39,61 @@
  * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
  */
 static const double
-two24  =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 pio2_1 =  1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
 pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
 pio2_3 =  6.36831716351370313614e-25; /*  0x18a2e037074000.0p-133 */
-
 static const long double
 invpio2 =  6.36619772367581343076e-01L, /*  0xa2f9836e4e44152a.0p-64 */
 pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
 pio2_2t =  6.36831716351095013979e-25L, /*  0xc51701b839a25205.0p-144 */
 pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
+#elif LDBL_MANT_DIG == 113
+/* u ~< 0x1p45*pi/2 */
+#define SMALL(u) (((u.i.se & 0x7fffU)<<16 | u.i.top) < ((0x3fff + 45)<<16 | 0x921f))
+#define TOINT 0x1.8p112
+#define QUOBITS(x) ((uint32_t)(int64_t)x & 0x7fffffff)
+#define ROUND1 51
+#define ROUND2 119
+#define NX 5
+#define NY 3
+static const long double
+invpio2 =  6.3661977236758134307553505349005747e-01L,	/*  0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
+pio2_1  =  1.5707963267948966192292994253909555e+00L,	/*  0x1921fb54442d18469800000000000.0p-112 */
+pio2_1t =  2.0222662487959507323996846200947577e-21L,	/*  0x13198a2e03707344a4093822299f3.0p-181 */
+pio2_2  =  2.0222662487959507323994779168837751e-21L,	/*  0x13198a2e03707344a400000000000.0p-181 */
+pio2_2t =  2.0670321098263988236496903051604844e-43L,	/*  0x127044533e63a0105df531d89cd91.0p-254 */
+pio2_3  =  2.0670321098263988236499468110329591e-43L,	/*  0x127044533e63a0105e00000000000.0p-254 */
+pio2_3t = -2.5650587247459238361625433492959285e-65L;	/* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
+#endif
 
 int __rem_pio2l(long double x, long double *y)
 {
-	union IEEEl2bits u,u1;
+	union ldshape u,uz;
 	long double z,w,t,r,fn;
-	double tx[3],ty[2];
-	int e0,ex,i,j,nx,n;
-	int16_t expsign;
-
-	u.e = x;
-	expsign = u.xbits.expsign;
-	ex = expsign & 0x7fff;
-	if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
-		union IEEEl2bits u2;
-		int ex1;
+	double tx[NX],ty[NY];
+	int ex,ey,n,i;
 
-		/* |x| ~< 2^25*(pi/2), medium size */
-		/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-		fn = x*invpio2 + 0x1.8p63;
-		fn = fn - 0x1.8p63;
-// FIXME
-//#ifdef HAVE_EFFICIENT_IRINT
-//		n = irint(fn);
-//#else
-		n = fn;
-//#endif
+	u.f = x;
+	ex = u.i.se & 0x7fff;
+	if (SMALL(u)) {
+		/* rint(x/(pi/2)), Assume round-to-nearest. */
+		fn = x*invpio2 + TOINT - TOINT;
+		n = QUOBITS(fn);
 		r = x-fn*pio2_1;
-		w = fn*pio2_1t;    /* 1st round good to 102 bit */
-		j = ex;
+		w = fn*pio2_1t;  /* 1st round good to 102/180 bits (ld80/ld128) */
 		y[0] = r-w;
-		u2.e = y[0];
-		ex1 = u2.xbits.expsign & 0x7fff;
-		i = j-ex1;
-		if (i > 22) {  /* 2nd iteration needed, good to 141 */
+		u.f = y[0];
+		ey = u.i.se & 0x7fff;
+		if (ex - ey > ROUND1) {  /* 2nd iteration needed, good to 141/248 (ld80/ld128) */
 			t = r;
 			w = fn*pio2_2;
 			r = t-w;
 			w = fn*pio2_2t-((t-r)-w);
 			y[0] = r-w;
-			u2.e = y[0];
-			ex1 = u2.xbits.expsign & 0x7fff;
-			i = j-ex1;
-			if (i > 61) {  /* 3rd iteration need, 180 bits acc */
-				t = r; /* will cover all possible cases */
+			u.f = y[0];
+			ey = u.i.se & 0x7fff;
+			if (ex - ey > ROUND2) {  /* 3rd iteration, good to 180/316 bits */
+				t = r; /* will cover all possible cases (not verified for ld128) */
 				w = fn*pio2_3;
 				r = t-w;
 				w = fn*pio2_3t-((t-r)-w);
@@ -102,23 +110,26 @@ int __rem_pio2l(long double x, long double *y)
 		y[0] = y[1] = x - x;
 		return 0;
 	}
-	/* set z = scalbn(|x|,ilogb(x)-23) */
-	u1.e = x;
-	e0 = ex - BIAS - 23;            /* e0 = ilogb(|x|)-23; */
-	u1.xbits.expsign = ex - e0;
-	z = u1.e;
-	for (i=0; i<2; i++) {
+	/* set z = scalbn(|x|,-ilogb(x)+23) */
+	uz.f = x;
+	uz.i.se = 0x3fff + 23;
+	z = uz.f;
+	for (i=0; i < NX - 1; i++) {
 		tx[i] = (double)(int32_t)z;
-		z     = (z-tx[i])*two24;
+		z     = (z-tx[i])*0x1p24;
 	}
-	tx[2] = z;
-	nx = 3;
-	while (tx[nx-1] == 0.0)
-		nx--;     /* skip zero term */
-	n = __rem_pio2_large(tx,ty,e0,nx,2);
-	r = (long double)ty[0] + ty[1];
-	w = ty[1] - (r - ty[0]);
-	if (expsign < 0) {
+	tx[i] = z;
+	while (tx[i] == 0)
+		i--;
+	n = __rem_pio2_large(tx, ty, ex-0x3fff-23, i+1, NY);
+	w = ty[1];
+	if (NY == 3)
+		w += ty[2];
+	r = ty[0] + w;
+	/* TODO: for ld128 this does not follow the recommendation of the
+	comments of __rem_pio2_large which seem wrong if |ty[0]| > |ty[1]+ty[2]| */
+	w -= r - ty[0];
+	if (u.i.se >> 15) {
 		y[0] = -r;
 		y[1] = -w;
 		return -n;

src/math/__sinl.c

@@ -1,4 +1,5 @@
 /* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */
+/* origin: FreeBSD /usr/src/lib/msun/ld128/k_sinl.c */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -13,7 +14,8 @@
 
 #include "libm.h"
 
-#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+#if LDBL_MANT_DIG == 64
 /*
  * ld80 version of __sin.c.  See __sin.c for most comments.
  */
@@ -23,10 +25,8 @@
  *
  * See __cosl.c for more details about the polynomial.
  */
-
 static const long double
 S1 = -0.166666666666666666671L;   /* -0xaaaaaaaaaaaaaaab.0p-66 */
-
 static const double
 S2 =  0.0083333333333333332,      /*  0x11111111111111.0p-59 */
 S3 = -0.00019841269841269427,     /* -0x1a01a01a019f81.0p-65 */
@@ -35,6 +35,34 @@ S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
 S6 =  1.6059006598854211e-10,     /*  0x161242b90243b5.0p-85 */
 S7 = -7.6429779983024564e-13,     /* -0x1ae42ebd1b2e00.0p-93 */
 S8 =  2.6174587166648325e-15;     /*  0x179372ea0b3f64.0p-101 */
+#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8))))))
+#elif LDBL_MANT_DIG == 113
+/*
+ * ld128 version of __sin.c.  See __sin.c for most comments.
+ */
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
+ * |sin(x)/x - s(x)| < 2**-122.1
+ *
+ * See __cosl.c for more details about the polynomial.
+ */
+static const long double
+S1 = -0.16666666666666666666666666666666666606732416116558L,
+S2 =  0.0083333333333333333333333333333331135404851288270047L,
+S3 = -0.00019841269841269841269841269839935785325638310428717L,
+S4 =  0.27557319223985890652557316053039946268333231205686e-5L,
+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
+S6 =  0.16059043836821614596571832194524392581082444805729e-9L,
+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
+S8 =  0.28114572543451292625024967174638477283187397621303e-14L;
+static const double
+S9  = -0.82206352458348947812512122163446202498005154296863e-17,
+S10 =  0.19572940011906109418080609928334380560135358385256e-19,
+S11 = -0.38680813379701966970673724299207480965452616911420e-22,
+S12 =  0.64038150078671872796678569586315881020659912139412e-25;
+#define POLY(z) (S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+ \
+	z*(S9+z*(S10+z*(S11+z*S12))))))))))
+#endif
 
 long double __sinl(long double x, long double y, int iy)
 {
@@ -42,7 +70,7 @@ long double __sinl(long double x, long double y, int iy)
 
 	z = x*x;
 	v = z*x;
-	r = S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*S8)))));
+	r = POLY(z);
 	if (iy == 0)
 		return x+v*(S1+z*r);
 	return x-((z*(0.5*y-v*r)-y)-v*S1);

src/math/__tanl.c

@@ -1,4 +1,5 @@
 /* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */
+/* origin: FreeBSD /usr/src/lib/msun/ld128/k_tanl.c */
 /*
  * ====================================================
  * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
@@ -12,7 +13,8 @@
 
 #include "libm.h"
 
-#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+#if LDBL_MANT_DIG == 64
 /*
  * ld80 version of __tan.c.  See __tan.c for most comments.
  */
@@ -22,14 +24,12 @@
  *
  * See __cosl.c for more details about the polynomial.
  */
-
 static const long double
 T3 =  0.333333333333333333180L,         /*  0xaaaaaaaaaaaaaaa5.0p-65 */
 T5 =  0.133333333333333372290L,         /*  0x88888888888893c3.0p-66 */
 T7 =  0.0539682539682504975744L,        /*  0xdd0dd0dd0dc13ba2.0p-68 */
 pio4   =  0.785398163397448309628L,     /*  0xc90fdaa22168c235.0p-64 */
 pio4lo = -1.25413940316708300586e-20L;  /* -0xece675d1fc8f8cbb.0p-130 */
-
 static const double
 T9  =  0.021869488536312216,            /*  0x1664f4882cc1c2.0p-58 */
 T11 =  0.0088632355256619590,           /*  0x1226e355c17612.0p-59 */
@@ -44,6 +44,59 @@ T27 =  0.0000024196006108814377,        /*  0x144c0d80cc6896.0p-71 */
 T29 =  0.0000078293456938132840,        /*  0x106b59141a6cb3.0p-69 */
 T31 = -0.0000032609076735050182,        /* -0x1b5abef3ba4b59.0p-71 */
 T33 =  0.0000023261313142559411;        /*  0x13835436c0c87f.0p-71 */
+#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
+	w * (T25 + w * (T29 + w * T33)))))))
+#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
+	w * (T27 + w * T31))))))
+#elif LDBL_MANT_DIG == 113
+/*
+ * ld128 version of __tan.c.  See __tan.c for most comments.
+ */
+/*
+ * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37]
+ * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37)
+ *
+ * See __cosl.c for more details about the polynomial.
+ */
+static const long double
+T3 = 0x1.5555555555555555555555555553p-2L,
+T5 = 0x1.1111111111111111111111111eb5p-3L,
+T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L,
+T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L,
+T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L,
+T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L,
+T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L,
+T17 = 0x1.355824803674477dfcf726649efep-11L,
+T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L,
+T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L,
+T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L,
+T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L,
+T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L,
+T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L,
+T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L,
+T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L,
+T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L,
+T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L,
+pio4 = 0x1.921fb54442d18469898cc51701b8p-1L,
+pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L;
+static const double
+T39 =  0.000000028443389121318352,	/*  0x1e8a7592977938.0p-78 */
+T41 =  0.000000011981013102001973,	/*  0x19baa1b1223219.0p-79 */
+T43 =  0.0000000038303578044958070,	/*  0x107385dfb24529.0p-80 */
+T45 =  0.0000000034664378216909893,	/*  0x1dc6c702a05262.0p-81 */
+T47 = -0.0000000015090641701997785,	/* -0x19ecef3569ebb6.0p-82 */
+T49 =  0.0000000029449552300483952,	/*  0x194c0668da786a.0p-81 */
+T51 = -0.0000000022006995706097711,	/* -0x12e763b8845268.0p-81 */
+T53 =  0.0000000015468200913196612,	/*  0x1a92fc98c29554.0p-82 */
+T55 = -0.00000000061311613386849674,	/* -0x151106cbc779a9.0p-83 */
+T57 =  1.4912469681508012e-10;		/*  0x147edbdba6f43a.0p-85 */
+#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
+	w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + \
+	w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))))
+#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
+	w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + \
+	w * (T47 + w * (T51 + w * T55))))))))))))
+#endif
 
 long double __tanl(long double x, long double y, int odd) {
 	long double z, r, v, w, s, a, t;
@@ -62,10 +115,8 @@ long double __tanl(long double x, long double y, int odd) {
 	}
 	z = x * x;
 	w = z * z;
-	r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
-	     w * (T25 + w * (T29 + w * T33))))));
-	v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
-	     w * (T27 + w * T31))))));
+	r = RPOLY(w);
+	v = z * VPOLY(w);
 	s = z * x;
 	r = y + z * (s * (r + v) + y) + T3 * s;
 	w = x + r;
@@ -76,7 +127,6 @@ long double __tanl(long double x, long double y, int odd) {
 	}
 	if (!odd)
 		return w;
-
 	/*
 	 * if allow error up to 2 ulp, simply return
 	 * -1.0 / (x+r) here

src/math/cosl.c

@@ -1,34 +1,3 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_cosl.c */
-/*-
- * Copyright (c) 2007 Steven G. Kargl
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice unmodified, this list of conditions, and the following
- *    disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-/*
- * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows
- * an accuracy of <= 0.7412 ULP.
- */
-
 #include "libm.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -38,44 +7,33 @@ long double cosl(long double x) {
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 long double cosl(long double x)
 {
-	union IEEEl2bits z;
+	union ldshape u = {x};
 	unsigned n;
-	long double y[2];
-	long double hi, lo;
-
-	z.e = x;
-	z.bits.sign = 0;
+	long double y[2], hi, lo;
 
-	/* If x = NaN or Inf, then cos(x) = NaN. */
-	if (z.bits.exp == 0x7fff)
-		return (x - x) / (x - x);
-
-	/* |x| < (double)pi/4 */
-	if (z.e < M_PI_4) {
-		/* |x| < 0x1p-64 */
-		if (z.bits.exp < 0x3fff - 64)
+	u.i.se &= 0x7fff;
+	if (u.i.se == 0x7fff)
+		return x - x;
+	x = u.f;
+	if (x < M_PI_4) {
+		if (u.i.se < 0x3fff - LDBL_MANT_DIG)
 			/* raise inexact if x!=0 */
 			return 1.0 + x;
-		return __cosl(z.e, 0);
+		return __cosl(x, 0);
 	}
-
 	n = __rem_pio2l(x, y);
 	hi = y[0];
 	lo = y[1];
 	switch (n & 3) {
 	case 0:
-		hi = __cosl(hi, lo);
-		break;
+		return __cosl(hi, lo);
 	case 1:
-		hi = -__sinl(hi, lo, 1);
-		break;
+		return -__sinl(hi, lo, 1);
 	case 2:
-		hi = -__cosl(hi, lo);
-		break;
+		return -__cosl(hi, lo);
 	case 3:
-		hi = __sinl(hi, lo, 1);
-		break;
+	default:
+		return __sinl(hi, lo, 1);
 	}
-	return hi;
 }
 #endif

src/math/sincosl.c

@@ -9,25 +9,19 @@ void sincosl(long double x, long double *sin, long double *cos)
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 void sincosl(long double x, long double *sin, long double *cos)
 {
-	union IEEEl2bits u;
+	union ldshape u = {x};
 	unsigned n;
 	long double y[2], s, c;
 
-	u.e = x;
-	u.bits.sign = 0;
-
-	/* x = nan or inf */
-	if (u.bits.exp == 0x7fff) {
+	u.i.se &= 0x7fff;
+	if (u.i.se == 0x7fff) {
 		*sin = *cos = x - x;
 		return;
 	}
-
-	/* |x| < (double)pi/4 */
-	if (u.e < M_PI_4) {
-		/* |x| < 0x1p-64 */
-		if (u.bits.exp < 0x3fff - 64) {
+	if (u.f < M_PI_4) {
+		if (u.i.se < 0x3fff - LDBL_MANT_DIG) {
 			/* raise underflow if subnormal */
-			if (u.bits.exp == 0) FORCE_EVAL(x*0x1p-120f);
+			if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f);
 			*sin = x;
 			/* raise inexact if x!=0 */
 			*cos = 1.0 + x;
@@ -37,7 +31,6 @@ void sincosl(long double x, long double *sin, long double *cos)
 		*cos = __cosl(x, 0);
 		return;
 	}
-
 	n = __rem_pio2l(x, y);
 	s = __sinl(y[0], y[1], 1);
 	c = __cosl(y[0], y[1]);

src/math/sinl.c

@@ -1,31 +1,3 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_sinl.c */
-/*-
- * Copyright (c) 2007 Steven G. Kargl
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice unmodified, this list of conditions, and the following
- *    disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-
 #include "libm.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -36,46 +8,34 @@ long double sinl(long double x)
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 long double sinl(long double x)
 {
-	union IEEEl2bits z;
+	union ldshape u = {x};
 	unsigned n;
-	long double y[2];
-	long double hi, lo;
-
-	z.e = x;
-	z.bits.sign = 0;
+	long double y[2], hi, lo;
 
-	/* If x = NaN or Inf, then sin(x) = NaN. */
-	if (z.bits.exp == 0x7fff)
-		return (x - x) / (x - x);
-
-	/* |x| < (double)pi/4 */
-	if (z.e < M_PI_4) {
-		/* |x| < 0x1p-64 */
-		if (z.bits.exp < 0x3fff - 64) {
+	u.i.se &= 0x7fff;
+	if (u.i.se == 0x7fff)
+		return x - x;
+	if (u.f < M_PI_4) {
+		if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
 			/* raise inexact if x!=0 and underflow if subnormal */
-			FORCE_EVAL(z.bits.exp == 0 ? x/0x1p120f : x+0x1p120f);
+			FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
 			return x;
 		}
 		return __sinl(x, 0.0, 0);
 	}
-
 	n = __rem_pio2l(x, y);
 	hi = y[0];
 	lo = y[1];
 	switch (n & 3) {
 	case 0:
-		hi = __sinl(hi, lo, 1);
-		break;
+		return __sinl(hi, lo, 1);
 	case 1:
-		hi = __cosl(hi, lo);
-		break;
+		return __cosl(hi, lo);
 	case 2:
-		hi = -__sinl(hi, lo, 1);
-		break;
+		return -__sinl(hi, lo, 1);
 	case 3:
-		hi = -__cosl(hi, lo);
-		break;
+	default:
+		return -__cosl(hi, lo);
 	}
-	return hi;
 }
 #endif

src/math/tanl.c

@@ -1,35 +1,3 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/s_tanl.c */
-/*-
- * Copyright (c) 2007 Steven G. Kargl
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * 1. Redistributions of source code must retain the above copyright
- *    notice unmodified, this list of conditions, and the following
- *    disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright
- *    notice, this list of conditions and the following disclaimer in the
- *    documentation and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
- * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
- * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
- * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
- * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
- * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
- * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
- * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
- * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-/*
- * Limited testing on pseudorandom numbers drawn within [0:4e8] shows
- * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million
- * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%).
- */
-
 #include "libm.h"
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
@@ -40,28 +8,21 @@ long double tanl(long double x)
 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 long double tanl(long double x)
 {
-	union IEEEl2bits z;
+	union ldshape u = {x};
 	long double y[2];
 	unsigned n;
 
-	z.e = x;
-	z.bits.sign = 0;
-
-	/* If x = NaN or Inf, then tan(x) = NaN. */
-	if (z.bits.exp == 0x7fff)
-		return (x - x) / (x - x);
-
-	/* |x| < (double)pi/4 */
-	if (z.e < M_PI_4) {
-		/* |x| < 0x1p-64 */
-		if (z.bits.exp < 0x3fff - 64) {
+	u.i.se &= 0x7fff;
+	if (u.i.se == 0x7fff)
+		return x - x;
+	if (u.f < M_PI_4) {
+		if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
 			/* raise inexact if x!=0 and underflow if subnormal */
-			FORCE_EVAL(z.bits.exp == 0 ? x/0x1p120f : x+0x1p120f);
+			FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
 			return x;
 		}
 		return __tanl(x, 0, 0);
 	}
-
 	n = __rem_pio2l(x, y);
 	return __tanl(y[0], y[1], n&1);
 }