Add implementations of sinl(), cosl(), and tanl().

Submitted by:	Steve Kargl <sgk@apl.washington.edu>

7 files changed, 367 insertions, 0 deletions

diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 41545eb..94bd919 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -405,7 +405,9 @@ long double	ceill(long double);
 long double	copysignl(long double, long double) __pure2;
 #if 0
 long double	coshl(long double);
+#endif
 long double	cosl(long double);
+#if 0
 long double	erfcl(long double);
 long double	erfl(long double);
 #endif
@@ -463,10 +465,14 @@ long double	scalblnl(long double, long);
 long double	scalbnl(long double, int);
 #if 0
 long double	sinhl(long double);
+#endif
 long double	sinl(long double);
+#if 0
 long double	sqrtl(long double);
 long double	tanhl(long double);
+#endif
 long double	tanl(long double);
+#if 0
 long double	tgammal(long double);
 #endif
 long double	truncl(long double);
diff --git a/lib/msun/src/s_cos.c b/lib/msun/src/s_cos.c
index 17afa34..b2fae38 100644
--- a/lib/msun/src/s_cos.c
+++ b/lib/msun/src/s_cos.c
@@ -45,6 +45,8 @@ static char rcsid[] = "$FreeBSD$";
  *	TRIG(x) returns trig(x) nearly rounded
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -80,3 +82,7 @@ cos(double x)
 	    }
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cos, cosl);
+#endif
diff --git a/lib/msun/src/s_cosl.c b/lib/msun/src/s_cosl.c
new file mode 100644
index 0000000..d184aa9
--- /dev/null
+++ b/lib/msun/src/s_cosl.c
@@ -0,0 +1,114 @@
+/*-
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * Compute cos(x) for x where x is reduced to y = x - k * pi / 2.
+ * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows
+ * an accuracy of <= 0.7412 ULP.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+cosl(long double x)
+{
+	union IEEEl2bits z;
+	int i, e0;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	z.bits.sign = 0;
+
+	/* If x = +-0 or x is a subnormal number, then cos(x) = 1 */
+	if (z.bits.exp == 0)
+		return (1.0);
+
+	/* If x = NaN or Inf, then cos(x) = NaN. */
+	if (z.bits.exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4)
+		return (__kernel_cosl(z.e, 0));
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+	
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	    hi = __kernel_cosl(hi, lo);
+	    break;
+	case 1:
+	    hi = - __kernel_sinl(hi, lo, 1);
+	    break;
+	case 2:
+	    hi = - __kernel_cosl(hi, lo);
+	    break;
+	case 3:
+	    hi = __kernel_sinl(hi, lo, 1);
+	    break;
+	}
+	
+	return (hi);
+}
diff --git a/lib/msun/src/s_sin.c b/lib/msun/src/s_sin.c
index af9ee62..7cd3877 100644
--- a/lib/msun/src/s_sin.c
+++ b/lib/msun/src/s_sin.c
@@ -45,6 +45,8 @@ static char rcsid[] = "$FreeBSD$";
  *	TRIG(x) returns trig(x) nearly rounded
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -80,3 +82,7 @@ sin(double x)
 	    }
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sin, sinl);
+#endif
diff --git a/lib/msun/src/s_sinl.c b/lib/msun/src/s_sinl.c
new file mode 100644
index 0000000..65d7026
--- /dev/null
+++ b/lib/msun/src/s_sinl.c
@@ -0,0 +1,115 @@
+/*-
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * Compute sin(x) for x where x is reduced to y = x - k * pi / 2.
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+sinl(long double x)
+{
+	union IEEEl2bits z;
+	int i, e0, s;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	s = z.bits.sign;
+	z.bits.sign = 0;
+
+	/* If x = +-0 or x is a subnormal number, then sin(x) = x */
+	if (z.bits.exp == 0)
+		return (x);
+
+	/* If x = NaN or Inf, then sin(x) = NaN. */
+	if (z.bits.exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4) {
+		hi = __kernel_sinl(z.e, 0, 0);
+		return  (s ? -hi : hi);
+	}
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	    hi = __kernel_sinl(hi, lo, 1);
+	    break;
+	case 1:
+	    hi = __kernel_cosl(hi, lo);
+	    break;
+	case 2:
+	    hi = - __kernel_sinl(hi, lo, 1);
+	    break;
+	case 3:
+	    hi = - __kernel_cosl(hi, lo);
+	    break;
+	}
+	
+	return (s ? -hi : hi);
+}
diff --git a/lib/msun/src/s_tan.c b/lib/msun/src/s_tan.c
index d707ecd..d9d37b6 100644
--- a/lib/msun/src/s_tan.c
+++ b/lib/msun/src/s_tan.c
@@ -44,6 +44,8 @@ static char rcsid[] = "$FreeBSD$";
  *	TRIG(x) returns trig(x) nearly rounded
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -74,3 +76,7 @@ tan(double x)
 							-1 -- n odd */
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(tan, tanl);
+#endif
diff --git a/lib/msun/src/s_tanl.c b/lib/msun/src/s_tanl.c
new file mode 100644
index 0000000..2592dcd
--- /dev/null
+++ b/lib/msun/src/s_tanl.c
@@ -0,0 +1,114 @@
+/*-
+ * 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 <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * Compute tan(x) for x where x is reduced to y = x - k * pi / 2.
+ * 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 <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+tanl(long double x)
+{
+	union IEEEl2bits z;
+	int i, e0, s;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	s = z.bits.sign;
+	z.bits.sign = 0;
+
+	/* If x = +-0 or x is subnormal, then tan(x) = x. */
+	if (z.bits.exp == 0)
+		return (x);
+
+	/* If x = NaN or Inf, then tan(x) = NaN. */
+	if (z.bits.exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4) {
+		hi = __kernel_tanl(z.e, 0, 0);
+		return (s ? -hi : hi);
+	}
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+	
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	case 2:
+	    hi = __kernel_tanl(hi, lo, 0);
+	    break;
+	case 1:
+	case 3:
+	    hi = __kernel_tanl(hi, lo, 1);
+	    break;
+	}
+
+	return (s ? -hi : hi);
+}