diff options
author | das <das@FreeBSD.org> | 2008-02-17 07:33:12 +0000 |
---|---|---|
committer | das <das@FreeBSD.org> | 2008-02-17 07:33:12 +0000 |
commit | 42e85f679f3473a2fdd4855f9fc6d109919b4ee6 (patch) | |
tree | 0107204fcf6b50c64cef34e3fb9b5079daec83ce /lib | |
parent | 61222ca5ae1148ade49b4c5249c16846498e35f7 (diff) | |
download | FreeBSD-src-42e85f679f3473a2fdd4855f9fc6d109919b4ee6.zip FreeBSD-src-42e85f679f3473a2fdd4855f9fc6d109919b4ee6.tar.gz |
Add implementations of sinl(), cosl(), and tanl().
Submitted by: Steve Kargl <sgk@apl.washington.edu>
Diffstat (limited to 'lib')
-rw-r--r-- | lib/msun/src/math.h | 6 | ||||
-rw-r--r-- | lib/msun/src/s_cos.c | 6 | ||||
-rw-r--r-- | lib/msun/src/s_cosl.c | 114 | ||||
-rw-r--r-- | lib/msun/src/s_sin.c | 6 | ||||
-rw-r--r-- | lib/msun/src/s_sinl.c | 115 | ||||
-rw-r--r-- | lib/msun/src/s_tan.c | 6 | ||||
-rw-r--r-- | lib/msun/src/s_tanl.c | 114 |
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); +} |