d/freebsd-nq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

[LIBM] implementations of sinpi[fl], cospi[fl], and tanpi[fl]

Browse Source 
Both IEEE-754 2008 and ISO/IEC TS 18661-4 define the half-cycle
trignometric functions cospi, sinpi, and tanpi.  The attached
patch implements cospi[fl], sinpi[fl], and tanpi[fl].  Limited
testing on the cospi and sinpi reveal a max ULP less than 0.89;
while tanpi is more problematic with a max ULP less than 2.01
in the interval [0,0.5].  The algorithms used in these functions
are documented in {ks}_cospi.c, {ks}_sinpi.c, and s_tanpi.c.

Note.  I no longer have access to a system with ld128 and
adequate support to compile and test the ld128 implementations
of these functions.  Given the almost complete lack of input from
others on improvements to libm, I doubt that anyone cares.  If
someone does care, the ld128 files contain a number of FIXME comments,
and in particular, while the polynomial coefficients are given
I did not update the polynomial algorithms to properly use the
coefficients.

PR:	218514
MFC after:	2 weeks
This commit is contained in:
Steve Kargl 2021-10-25 16:13:52 +03:00 committed by Konstantin Belousov
parent 179219ea04
commit dce5f3abed
24 changed files with 2203 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
lib/msun/Makefile
View File

@ -126,6 +126,12 @@ COMMON_SRCS+= catrigl.c \
# See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=130067
.if ${COMPILER_TYPE} == "gcc"
CFLAGS.e_powl.c+= -Wno-error=overflow
# IEEE-754 2008 and ISO/IEC TS 18661-4 half-cycle trignometric functions
COMMON_SRCS+= s_cospi.c s_cospif.c s_cospil.c \
s_sinpi.c s_sinpif.c s_sinpil.c \
s_tanpi.c s_tanpif.c s_tanpil.c
.endif
.endif
@ -154,7 +160,8 @@ INCS+= fenv.h math.h
MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
cimag.3 clog.3 copysign.3 cos.3 cosh.3 cpow.3 csqrt.3 erf.3 \
cimag.3 clog.3 copysign.3 cos.3 cosh.3 cospi.3 \
cpow.3 csqrt.3 erf.3 \
exp.3 fabs.3 fdim.3 \
feclearexcept.3 feenableexcept.3 fegetenv.3 \
fegetround.3 fenv.3 floor.3 \
@ -162,7 +169,7 @@ MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
nextafter.3 remainder.3 rint.3 \
round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
sinh.3 sinpi.3 sqrt.3 tan.3 tanh.3 tanpi.3 trunc.3 \
complex.3
MLINKS+=acos.3 acosf.3 acos.3 acosl.3
@ -192,6 +199,7 @@ MLINKS+=clog.3 clogf.3 clog.3 clogl.3
MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
MLINKS+=cos.3 cosf.3 cos.3 cosl.3
MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
MLINKS+=cospi.3 cospif.3 cospi.3 cospil.3
MLINKS+=cpow.3 cpowf.3 cpow.3 cpowl.3
MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
@ -244,10 +252,12 @@ MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
MLINKS+=sin.3 sinf.3 sin.3 sinl.3
MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
MLINKS+=sinpi.3 sinpif.3 sinpi.3 sinpil.3
MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
sqrt.3 sqrtl.3
MLINKS+=tan.3 tanf.3 tan.3 tanl.3
MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3
MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
.include <src.opts.mk>

13
lib/msun/Symbol.map
View File

@ -304,3 +304,16 @@ FBSD_1.5 {
sincosf;
sincosl;
};
/* First added in 14.0-CURRENT */
FBSD_1.7 {
cospi;
cospif;
cospil;
sinpi;
sinpif;
sinpil;
tanpi;
tanpif;
tanpil;
};

42
lib/msun/ld128/k_cospil.h Normal file
View File

@ -0,0 +1,42 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/k_cospi.c for implementation details.
*/
static inline long double
__kernel_cospil(long double x)
{
long double hi, lo;
hi = (double)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_cosl(hi, lo));
}

42
lib/msun/ld128/k_sinpil.h Normal file
View File

@ -0,0 +1,42 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/k_sinpi.c for implementation details.
*/
static inline long double
__kernel_sinpil(long double x)
{
long double hi, lo;
hi = (double)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_sinl(hi, lo, 1));
}

109
lib/msun/ld128/s_cospil.c Normal file
View File

@ -0,0 +1,109 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_cospi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
#include "math_private.h"
/*
* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
*/
static const long double
pi_hi = 3.14159265358979322702026593105983920e+00L,
pi_lo = 1.14423774522196636802434264184180742e-17L;
#include "k_cospil.h"
#include "k_sinpil.h"
volatile static const double vzero = 0;
long double
cospil(long double x)
{
long double ax, c, xf;
uint32_t ix;
ax = fabsl(x);
if (ax < 1) {
if (ax < 0.25) {
if (ax < 0x1p-60) {
if ((int)x == 0)
return (1);
}
return (__kernel_cospil(ax));
}
if (ax < 0.5)
c = __kernel_sinpil(0.5 - ax);
else if (ax < 0.75) {
if (ax == 0.5)
return (0);
c = -__kernel_sinpil(ax - 0.5);
} else
c = -__kernel_cospil(1 - ax);
return (c);
}
if (ax < 0x1p112) {
xf = floorl(ax);
ax -= xf;
if (x < 0.5) {
if (x < 0.25)
c = ax == 0 ? 1 : __kernel_cospil(ax);
else
c = __kernel_sinpil(0.5 - ax);
} else {
if (x < 0.75) {
if (ax == 0.5)
return (0);
c = -__kernel_sinpil(ax - 0.5);
} else
c = -__kernel_cospil(1 - ax);
}
if (xf > 0x1p50)
xf -= 0x1p50;
if (xf > 0x1p30)
xf -= 0x1p30;
ix = (uint32_t)xf;
return (ix & 1 ? -c : c);
}
if (isinf(x) || isnan(x))
return (vzero / vzero);
/*
* |x| >= 0x1p112 is always an even integer, so return 1.
*/
return (1);
}

118
lib/msun/ld128/s_sinpil.c Normal file
View File

@ -0,0 +1,118 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_sinpi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
#include "math_private.h"
/*
* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
*/
static const long double
pi_hi = 3.14159265358979322702026593105983920e+00L,
pi_lo = 1.14423774522196636802434264184180742e-17L;
#include "k_cospil.h"
#include "k_sinpil.h"
volatile static const double vzero = 0;
long double
sinpil(long double x)
{
long double ax, hi, lo, s, xf, xhi, xlo;
uint32_t ix;
ax = fabsl(x);
if (ax < 1) {
if (ax < 0.25) {
if (ax < 0x1p-60) {
if (x == 0)
return (x);
hi = (double)x;
hi *= 0x1p113L;
lo = x * 0x1p113L - hi;
s = (pi_lo + pi_hi) * lo + pi_lo * lo +
pi_hi * hi;
return (s * 0x1p-113L);
}
s = __kernel_sinpil(ax);
return (copysignl(s, x));
}
if (ax < 0.5)
s = __kernel_cospil(0.5 - ax);
else if (ax < 0.75)
s = __kernel_cospil(ax - 0.5);
else
s = __kernel_sinpil(1 - ax);
return (copysignl(s, x));
}
if (ax < 0x1p112) {
xf = floorl(ax);
ax -= xf;
if (ax == 0) {
s = 0;
} else {
if (ax < 0.5) {
if (ax <= 0.25)
s = __kernel_sinpil(ax);
else
s = __kernel_cospil(0.5 - ax);
} else {
if (ax < 0.75)
s = __kernel_cospil(ax - 0.5);
else
s = __kernel_sinpil(1 - ax);
}
if (xf > 0x1p50)
xf -= 0x1p50;
if (xf > 0x1p30)
xf -= 0x1p30;
ix = (uint32_t)xf;
if (ix & 1) s = -s;
}
return (copysignl(s, x));
}
if (isinf(x) || isnan(x))
return (vzero / vzero);
/*
* |x| >= 0x1p112 is always an integer, so return +-0.
*/
return (copysignl(0, x));
}

120
lib/msun/ld128/s_tanpil.c Normal file
View File

@ -0,0 +1,120 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_tanpi.c for implementation details.
*
* FIXME: This has not been compiled nor has it been tested for accuracy.
* FIXME: This should use bit twiddling.
*/
#include "math.h"
#include "math_private.h"
/*
* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
*/
static const long double
pi_hi = 3.14159265358979322702026593105983920e+00L,
pi_lo = 1.14423774522196636802434264184180742e-17L;
static inline long double
__kernel_tanpi(long double x)
{
long double hi, lo, t;
if (x < 0.25) {
hi = (double)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = __kernel_tanl(hi, lo, -1);
} else if (x > 0.25) {
x = 0.5 - x;
hi = (double)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = - __kernel_tanl(hi, lo, 1);
} else
t = 1;
return (t);
}
volatile static const double vzero = 0;
long double
tanpil(long double x)
{
long double ax, hi, lo, xf;
uint32_t ix;
ax = fabsl(ax);
if (ax < 1) {
if (ax < 0.5) {
if (ax < 0x1p-60) {
if (x == 0)
return (x);
hi = (double)x;
hi *= 0x1p113L
lo = x * 0x1p113L - hi;
t = (pi_lo + pi_hi) * lo + pi_lo * lo +
pi_hi * hi;
return (t * 0x1p-113L);
}
t = __kernel_tanpil(ax);
} else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
else
t = -__kernel_tanpil(1 - ax);
return (copysignl(t, x));
}
if (ix < 0x1p112) {
xf = floorl(ax);
ax -= xf;
if (ax < 0.5)
t = ax == 0 ? 0 : __kernel_tanpil(ax);
else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
else
t = -__kernel_tanpil(1 - ax);
return (copysignl(t, x));
}
/* x = +-inf or nan. */
if (isinf(x) || isnan(x))
return (vzero / vzero);
/*
* |x| >= 0x1p53 is always an integer, so return +-0.
*/
return (copysignl(0, x));
}

42
lib/msun/ld80/k_cospil.h Normal file
View File

@ -0,0 +1,42 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/k_cospi.c for implementation details.
*/
static inline long double
__kernel_cospil(long double x)
{
long double hi, lo;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_cosl(hi, lo));
}

42
lib/msun/ld80/k_sinpil.h Normal file
View File

@ -0,0 +1,42 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/k_sinpi.c for implementation details.
*/
static inline long double
__kernel_sinpil(long double x)
{
long double hi, lo;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_sinl(hi, lo, 1));
}

129
lib/msun/ld80/s_cospil.c Normal file
View File

@ -0,0 +1,129 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_cospi.c for implementation details.
*/
#ifdef __i386__
#include <ieeefp.h>
#endif
#include <stdint.h>
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
static const double
pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
#include "k_cospil.h"
#include "k_sinpil.h"
volatile static const double vzero = 0;
long double
cospil(long double x)
{
long double ax, c;
uint64_t lx, m;
uint32_t j0;
uint16_t hx, ix;
EXTRACT_LDBL80_WORDS(hx, lx, x);
ix = hx & 0x7fff;
INSERT_LDBL80_WORDS(ax, ix, lx);
ENTERI();
if (ix < 0x3fff) { /* |x| < 1 */
if (ix < 0x3ffd) { /* |x| < 0.25 */
if (ix < 0x3fdd) { /* |x| < 0x1p-34 */
if ((int)x == 0)
RETURNI(1);
}
RETURNI(__kernel_cospil(ax));
}
if (ix < 0x3ffe) /* |x| < 0.5 */
c = __kernel_sinpil(0.5 - ax);
else if (lx < 0xc000000000000000ull) { /* |x| < 0.75 */
if (ax == 0.5)
RETURNI(0);
c = -__kernel_sinpil(ax - 0.5);
} else
c = -__kernel_cospil(1 - ax);
RETURNI(c);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);
if (ix < 0x3ffe) { /* |x| < 0.5 */
if (ix < 0x3ffd) /* |x| < 0.25 */
c = ix == 0 ? 1 : __kernel_cospil(ax);
else
c = __kernel_sinpil(0.5 - ax);
} else {
if (lx < 0xc000000000000000ull) { /* |x| < 0.75 */
if (ax == 0.5)
RETURNI(0);
c = -__kernel_sinpil(ax - 0.5);
} else
c = -__kernel_cospil(1 - ax);
}
if (j0 > 40)
x -= 0x1p40;
if (j0 > 30)
x -= 0x1p30;
j0 = (uint32_t)x;
RETURNI(j0 & 1 ? -c : c);
}
if (ix >= 0x7fff)
RETURNI(vzero / vzero);
/*
* |x| >= 0x1p63 is always an even integer, so return 1.
*/
RETURNI(1);
}

140
lib/msun/ld80/s_sinpil.c Normal file
View File

@ -0,0 +1,140 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_sinpi.c for implementation details.
*/
#ifdef __i386__
#include <ieeefp.h>
#endif
#include <stdint.h>
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
static const union IEEEl2bits
pi_hi_u = LD80C(0xc90fdaa200000000, 1, 3.14159265346825122833e+00L),
pi_lo_u = LD80C(0x85a308d313198a2e, -33, 1.21542010130123852029e-10L);
#define pi_hi (pi_hi_u.e)
#define pi_lo (pi_lo_u.e)
#include "k_cospil.h"
#include "k_sinpil.h"
volatile static const double vzero = 0;
long double
sinpil(long double x)
{
long double ax, hi, lo, s;
uint64_t lx, m;
uint32_t j0;
uint16_t hx, ix;
EXTRACT_LDBL80_WORDS(hx, lx, x);
ix = hx & 0x7fff;
INSERT_LDBL80_WORDS(ax, ix, lx);
ENTERI();
if (ix < 0x3fff) { /* |x| < 1 */
if (ix < 0x3ffd) { /* |x| < 0.25 */
if (ix < 0x3fdd) { /* |x| < 0x1p-34 */
if (x == 0)
RETURNI(x);
INSERT_LDBL80_WORDS(hi, hx,
lx & 0xffffffff00000000ull);
hi *= 0x1p63L;
lo = x * 0x1p63L - hi;
s = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
RETURNI(s * 0x1p-63L);
}
s = __kernel_sinpil(ax);
RETURNI((hx & 0x8000) ? -s : s);
}
if (ix < 0x3ffe) /* |x| < 0.5 */
s = __kernel_cospil(0.5 - ax);
else if (lx < 0xc000000000000000ull) /* |x| < 0.75 */
s = __kernel_cospil(ax - 0.5);
else
s = __kernel_sinpil(1 - ax);
RETURNI((hx & 0x8000) ? -s : s);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);
if (ix == 0) {
s = 0;
} else {
if (ix < 0x3ffe) { /* |x| < 0.5 */
if (ix < 0x3ffd) /* |x| < 0.25 */
s = __kernel_sinpil(ax);
else
s = __kernel_cospil(0.5 - ax);
} else {
/* |x| < 0.75 */
if (lx < 0xc000000000000000ull)
s = __kernel_cospil(ax - 0.5);
else
s = __kernel_sinpil(1 - ax);
}
if (j0 > 40)
x -= 0x1p40;
if (j0 > 30)
x -= 0x1p30;
j0 = (uint32_t)x;
if (j0 & 1) s = -s;
}
RETURNI((hx & 0x8000) ? -s : s);
}
/* x = +-inf or nan. */
if (ix >= 0x7fff)
RETURNI(vzero / vzero);
/*
* |x| >= 0x1p63 is always an integer, so return +-0.
*/
RETURNI(copysignl(0, x));
}

139
lib/msun/ld80/s_tanpil.c Normal file
View File

@ -0,0 +1,139 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_tanpi.c for implementation details.
*/
#ifdef __i386__
#include <ieeefp.h>
#endif
#include <stdint.h>
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
static const double
pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
pi_lo = -2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
static inline long double
__kernel_tanpil(long double x)
{
long double hi, lo, t;
if (x < 0.25) {
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = __kernel_tanl(hi, lo, -1);
} else if (x > 0.25) {
x = 0.5 - x;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = - __kernel_tanl(hi, lo, 1);
} else
t = 1;
return (t);
}
volatile static const double vzero = 0;
long double
tanpil(long double x)
{
long double ax, hi, lo, t;
uint64_t lx, m;
uint32_t j0;
uint16_t hx, ix;
EXTRACT_LDBL80_WORDS(hx, lx, x);
ix = hx & 0x7fff;
INSERT_LDBL80_WORDS(ax, ix, lx);
ENTERI();
if (ix < 0x3fff) { /* |x| < 1 */
if (ix < 0x3ffe) { /* |x| < 0.5 */
if (ix < 0x3fdd) { /* |x| < 0x1p-34 */
if (x == 0)
RETURNI(x);
INSERT_LDBL80_WORDS(hi, hx,
lx & 0xffffffff00000000ull);
hi *= 0x1p63L;
lo = x * 0x1p63L - hi;
t = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
RETURNI(t * 0x1p-63L);
}
t = __kernel_tanpil(ax);
} else if (ax == 0.5)
RETURNI((ax - ax) / (ax - ax));
else
t = -__kernel_tanpil(1 - ax);
RETURNI((hx & 0x8000) ? -t : t);
}
if (ix < 0x403e) { /* 1 <= |x| < 0x1p63 */
/* Determine integer part of ax. */
j0 = ix - 0x3fff + 1;
if (j0 < 32) {
lx = (lx >> 32) << 32;
lx &= ~(((lx << 32)-1) >> j0);
} else {
m = (uint64_t)-1 >> (j0 + 1);
if (lx & m) lx &= ~m;
}
INSERT_LDBL80_WORDS(x, ix, lx);
ax -= x;
EXTRACT_LDBL80_WORDS(ix, lx, ax);
if (ix < 0x3ffe) /* |x| < 0.5 */
t = ax == 0 ? 0 : __kernel_tanpil(ax);
else if (ax == 0.5)
RETURNI((ax - ax) / (ax - ax));
else
t = -__kernel_tanpil(1 - ax);
RETURNI((hx & 0x8000) ? -t : t);
}
/* x = +-inf or nan. */
if (ix >= 0x7fff)
RETURNI(vzero / vzero);
/*
* |x| >= 0x1p63 is always an integer, so return +-0.
*/
RETURNI(copysignl(0, x));
}

111
lib/msun/man/cospi.3 Normal file
View File

@ -0,0 +1,111 @@
.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" 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, 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 REGENTS AND CONTRIBUTORS ``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 REGENTS OR CONTRIBUTORS 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.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt COSPI 3
.Os
.Sh NAME
.Nm cospi ,
.Nm cospif ,
.Nm cospil
.Nd half\(encycle cosine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn cospi "double x"
.Ft float
.Fn cospif "float x"
.Ft long double
.Fn cospil "long double x"
.Sh DESCRIPTION
The
.Fn cospi ,
.Fn cospif ,
and
.Fn cospil
functions compute the cosine of
.Fa "\(*p \(mu x" .
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn cospi ,
.Fn cospif ,
and
.Fn cospil
functions returns
.Fn cos "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of
.Ar x ,
then the returned value is 1 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn cospi \*(Pm0
returns 1.
.It
.Fn cospi \*(Pmn/2
returns 0 for positive integers
.Ar n .
.It
.Fn cospi n
returns 1 for even integers
.Ar n .
.It
.Fn cospi n
returns \-1 for odd integers
.Ar n .
.It
.Fn cospi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn cospi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr sinpi 3 ,
.Xr tan 3 ,
.Xr tanpi 3
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 ,
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.

102
lib/msun/man/sinpi.3 Normal file
View File

@ -0,0 +1,102 @@
.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" 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, 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 REGENTS AND CONTRIBUTORS ``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 REGENTS OR CONTRIBUTORS 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.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt SINPI 3
.Os
.Sh NAME
.Nm sinpi ,
.Nm sinpif ,
.Nm sinpil
.Nd half\(encycle sine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn sinpi "double x"
.Ft float
.Fn sinpif "float x"
.Ft long double
.Fn sinpil "long double x"
.Sh DESCRIPTION
The
.Fn sinpi ,
.Fn sinpif ,
and
.Fn sinpil
functions compute the sine of
.Fa "\(*p \(mu x" .
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn sinpi ,
.Fn sinpif ,
and
.Fn sinpil
functions returns
.Fn sin "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of
.Ar x ,
then the returned value is \*(Pm0 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn sinpi \*(Pm0
returns \*(Pm0.
.It
.Fn sinpi \*(Pmn
returns \*(Pm0 for positive integers
.Ar n .
.It
.Fn sinpi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn sinpi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr cospi 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr tan 3 ,
.Xr tanpi 3
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 ,
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.

106
lib/msun/man/tanpi.3 Normal file
View File

@ -0,0 +1,106 @@
.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" 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, 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 REGENTS AND CONTRIBUTORS ``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 REGENTS OR CONTRIBUTORS 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.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt TANPI 3
.Os
.Sh NAME
.Nm tanpi ,
.Nm tanpif ,
.Nm tanpil
.Nd half\(encycle tangent functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn tanpi "double x"
.Ft float
.Fn tanpif "float x"
.Ft long double
.Fn tanpil "long double x"
.Sh DESCRIPTION
The
.Fn tanpi ,
.Fn tanpif ,
and
.Fn tanpil
functions compute the tangent of
.Fa "\(*p \(mu x"
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn tanpi ,
.Fn tanpif ,
and
.Fn tanpil
functions returns
.Fn tan "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of
.Ar x ,
then the returned value is \*(Pm0 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn tanpi \*(Pm0
returns \*(Pm0.
.It
.Fn tanpi \*(Pmn
returns \*(Pm0 for positive integers
.Ar n .
.It
.Fn tanpi \*(Pmn/2
returns \*(Na for n > 0 and raises an FE_INVALID exception.
.It
.Fn tanpi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn tanpi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr cospi 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr sinpi 3 ,
.Xr tan 3 ,
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 ,
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.

44
lib/msun/src/k_cospi.h Normal file
View File

@ -0,0 +1,44 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* The basic kernel for x in [0,0.25]. To use the kernel for cos(x), the
* argument to __kernel_cospi() must be multiplied by pi.
*/
static inline double
__kernel_cospi(double x)
{
double_t hi, lo;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_cos(hi, lo));
}

43
lib/msun/src/k_sinpi.h Normal file
View File

@ -0,0 +1,43 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* The basic kernel for x in [0,0.25]. To use the kernel for sin(x), the
* argument to __kernel_sinpi() must be multiplied by pi.
*/
static inline double
__kernel_sinpi(double x)
{
double_t hi, lo;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
return (__kernel_sin(hi, lo, 1));
}

9
lib/msun/src/math.h
View File

@ -508,6 +508,15 @@ long double lgammal_r(long double, int *);
void sincos(double, double *, double *);
void sincosf(float, float *, float *);
void sincosl(long double, long double *, long double *);
double cospi(double);
float cospif(float);
long double cospil(long double);
double sinpi(double);
float sinpif(float);
long double sinpil(long double);
double tanpi(double);
float tanpif(float);
long double tanpil(long double);
#endif
__END_DECLS

151
lib/msun/src/s_cospi.c Normal file
View File

@ -0,0 +1,151 @@
/*-
* Copyright (c) 2017 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.
*/
/**
* cospi(x) computes cos(pi*x) without multiplication by pi (almost). First,
* note that cospi(-x) = cospi(x), so the algorithm considers only |x|. The
* method used depends on the magnitude of x.
*
* 1. For small |x|, cospi(x) = 1 with FE_INEXACT raised where a sloppy
* threshold is used. The threshold is |x| < 0x1pN with N = -(P/2+M).
* P is the precision of the floating-point type and M = 2 to 4.
*
* 2. For |x| < 1, argument reduction is not required and sinpi(x) is
* computed by calling a kernel that leverages the kernels for sin(x)
* ans cos(x). See k_sinpi.c and k_cospi.c for details.
*
* 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
* |x| = j0 + r with j0 an integer and the remainder r satisfies
* 0 <= r < 1. With the given domain, a simplified inline floor(x)
* is used. Also, note the following identity
*
* cospi(x) = cos(pi*(j0+r))
* = cos(pi*j0) * cos(pi*r) - sin(pi*j0) * sin(pi*r)
* = cos(pi*j0) * cos(pi*r)
* = +-cospi(r)
*
* If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
* cospi(r) is then computed via an appropriate kernel.
*
* 4. For |x| >= 0x1p(P-1), |x| is integral and cospi(x) = 1.
*
* 5. Special cases:
*
* cospi(+-0) = 1.
* cospi(n.5) = 0 for n an integer.
* cospi(+-inf) = nan. Raises the "invalid" floating-point exception.
* cospi(nan) = nan. Raises the "invalid" floating-point exception.
*/
#include "math.h"
#include "math_private.h"
static const double
pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
#include "k_cospi.h"
#include "k_sinpi.h"
volatile static const double vzero = 0;
double
cospi(double x)
{
double ax, c;
uint32_t hx, ix, j0, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
INSERT_WORDS(ax, ix, lx);
if (ix < 0x3ff00000) { /* |x| < 1 */
if (ix < 0x3fd00000) { /* |x| < 0.25 */
if (ix < 0x3e200000) { /* |x| < 0x1p-29 */
if ((int)ax == 0)
return (1);
}
return (__kernel_cospi(ax));
}
if (ix < 0x3fe00000) /* |x| < 0.5 */
c = __kernel_sinpi(0.5 - ax);
else if (ix < 0x3fe80000){ /* |x| < 0.75 */
if (ax == 0.5)
return (0);
c = -__kernel_sinpi(ax - 0.5);
} else
c = -__kernel_cospi(1 - ax);
return (c);
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix < 0x3fd00000) /* |x| < 0.25 */
c = ix == 0 ? 1 : __kernel_cospi(ax);
else
c = __kernel_sinpi(0.5 - ax);
} else {
if (ix < 0x3fe80000) { /* |x| < 0.75 */
if (ax == 0.5)
return (0);
c = -__kernel_sinpi(ax - 0.5);
} else
c = -__kernel_cospi(1 - ax);
}
if (j0 > 30)
x -= 0x1p30;
j0 = (uint32_t)x;
return (j0 & 1 ? -c : c);
}
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an even integer, so return 1.
*/
return (1);
}
#if LDBL_MANT_DIG == 53
__weak_reference(cospi, cospil);
#endif

109
lib/msun/src/s_cospif.c Normal file
View File

@ -0,0 +1,109 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_cospi.c for implementation details.
*/
#define INLINE_KERNEL_SINDF
#define INLINE_KERNEL_COSDF
#include "math.h"
#include "math_private.h"
#include "k_cosf.c"
#include "k_sinf.c"
#define __kernel_cospif(x) (__kernel_cosdf(M_PI * (x)))
#define __kernel_sinpif(x) (__kernel_sindf(M_PI * (x)))
volatile static const float vzero = 0;
float
cospif(float x)
{
float ax, c;
uint32_t ix, j0;
GET_FLOAT_WORD(ix, x);
ix = ix & 0x7fffffff;
SET_FLOAT_WORD(ax, ix);
if (ix < 0x3f800000) { /* |x| < 1 */
if (ix < 0x3e800000) { /* |x| < 0.25 */
if (ix < 0x38800000) { /* |x| < 0x1p-14 */
/* Raise inexact iff != 0. */
if ((int)ax == 0)
return (1);
}
return (__kernel_cospif(ax));
}
if (ix < 0x3f000000) /* |x| < 0.5 */
c = __kernel_sinpif(0.5F - ax);
else if (ix < 0x3f400000) { /* |x| < 0.75 */
if (ix == 0x3f000000)
return (0);
c = -__kernel_sinpif(ax - 0.5F);
} else
c = -__kernel_cospif(1 - ax);
return (c);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
ax -= x;
GET_FLOAT_WORD(ix, ax);
if (ix < 0x3f000000) { /* |x| < 0.5 */
if (ix < 0x3e800000) /* |x| < 0.25 */
c = ix == 0 ? 1 : __kernel_cospif(ax);
else
c = __kernel_sinpif(0.5F - ax);
} else {
if (ix < 0x3f400000) { /* |x| < 0.75 */
if (ix == 0x3f000000)
return (0);
c = -__kernel_sinpif(ax - 0.5F);
} else
c = -__kernel_cospif(1 - ax);
}
j0 = (uint32_t)x;
return (j0 & 1 ? -c : c);
}
/* x = +-inf or nan. */
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p23 is always an even integer, so return 1.
*/
return (1);
}

168
lib/msun/src/s_sinpi.c Normal file
View File

@ -0,0 +1,168 @@
/*-
* Copyright (c) 2017 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.
*/
/**
* sinpi(x) computes sin(pi*x) without multiplication by pi (almost). First,
* note that sinpi(-x) = -sinpi(x), so the algorithm considers only |x| and
* includes reflection symmetry by considering the sign of x on output. The
* method used depends on the magnitude of x.
*
* 1. For small |x|, sinpi(x) = pi * x where a sloppy threshold is used. The
* threshold is |x| < 0x1pN with N = -(P/2+M). P is the precision of the
* floating-point type and M = 2 to 4. To achieve high accuracy, pi is
* decomposed into high and low parts with the high part containing a
* number of trailing zero bits. x is also split into high and low parts.
*
* 2. For |x| < 1, argument reduction is not required and sinpi(x) is
* computed by calling a kernel that leverages the kernels for sin(x)
* ans cos(x). See k_sinpi.c and k_cospi.c for details.
*
* 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
* |x| = j0 + r with j0 an integer and the remainder r satisfies
* 0 <= r < 1. With the given domain, a simplified inline floor(x)
* is used. Also, note the following identity
*
* sinpi(x) = sin(pi*(j0+r))
* = sin(pi*j0) * cos(pi*r) + cos(pi*j0) * sin(pi*r)
* = cos(pi*j0) * sin(pi*r)
* = +-sinpi(r)
*
* If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
* sinpi(r) is then computed via an appropriate kernel.
*
* 4. For |x| >= 0x1p(P-1), |x| is integral and sinpi(x) = copysign(0,x).
*
* 5. Special cases:
*
* sinpi(+-0) = +-0
* sinpi(+-n) = +-0, for positive integers n.
* sinpi(+-inf) = nan. Raises the "invalid" floating-point exception.
* sinpi(nan) = nan. Raises the "invalid" floating-point exception.
*/
#include "math.h"
#include "math_private.h"
static const double
pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
#include "k_cospi.h"
#include "k_sinpi.h"
volatile static const double vzero = 0;
double
sinpi(double x)
{
double ax, hi, lo, s;
uint32_t hx, ix, j0, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
INSERT_WORDS(ax, ix, lx);
if (ix < 0x3ff00000) { /* |x| < 1 */
if (ix < 0x3fd00000) { /* |x| < 0.25 */
if (ix < 0x3e200000) { /* |x| < 0x1p-29 */
if (x == 0)
return (x);
/*
* To avoid issues with subnormal values,
* scale the computation and rescale on
* return.
*/
INSERT_WORDS(hi, hx, 0);
hi *= 0x1p53;
lo = x * 0x1p53 - hi;
s = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
return (s * 0x1p-53);
}
s = __kernel_sinpi(ax);
return ((hx & 0x80000000) ? -s : s);
}
if (ix < 0x3fe00000) /* |x| < 0.5 */
s = __kernel_cospi(0.5 - ax);
else if (ix < 0x3fe80000) /* |x| < 0.75 */
s = __kernel_cospi(ax - 0.5);
else
s = __kernel_sinpi(1 - ax);
return ((hx & 0x80000000) ? -s : s);
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
}
INSERT_WORDS(x, ix, lx);
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix == 0)
s = 0;
else {
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix < 0x3fd00000) /* |x| < 0.25 */
s = __kernel_sinpi(ax);
else
s = __kernel_cospi(0.5 - ax);
} else {
if (ix < 0x3fe80000) /* |x| < 0.75 */
s = __kernel_cospi(ax - 0.5);
else
s = __kernel_sinpi(1 - ax);
}
if (j0 > 30)
x -= 0x1p30;
j0 = (uint32_t)x;
if (j0 & 1) s = -s;
}
return ((hx & 0x80000000) ? -s : s);
}
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an integer, so return +-0.
*/
return (copysign(0, x));
}
#if LDBL_MANT_DIG == 53
__weak_reference(sinpi, sinpil);
#endif

122
lib/msun/src/s_sinpif.c Normal file
View File

@ -0,0 +1,122 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_sinpi.c for implementation details.
*/
#define INLINE_KERNEL_SINDF
#define INLINE_KERNEL_COSDF
#include "math.h"
#include "math_private.h"
#include "k_cosf.c"
#include "k_sinf.c"
#define __kernel_cospif(x) (__kernel_cosdf(M_PI * (x)))
#define __kernel_sinpif(x) (__kernel_sindf(M_PI * (x)))
static const float
pi_hi = 3.14160156e+00F, /* 0x40491000 */
pi_lo = -8.90890988e-06F; /* 0xb715777a */
volatile static const float vzero = 0;
float
sinpif(float x)
{
float ax, hi, lo, s;
uint32_t hx, ix, j0;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
SET_FLOAT_WORD(ax, ix);
if (ix < 0x3f800000) { /* |x| < 1 */
if (ix < 0x3e800000) { /* |x| < 0.25 */
if (ix < 0x38800000) { /* |x| < 0x1p-14 */
if (x == 0)
return (x);
SET_FLOAT_WORD(hi, hx & 0xffff0000);
hi *= 0x1p23F;
lo = x * 0x1p23F - hi;
s = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
return (s * 0x1p-23F);
}
s = __kernel_sinpif(ax);
return ((hx & 0x80000000) ? -s : s);
}
if (ix < 0x3f000000) /* |x| < 0.5 */
s = __kernel_cospif(0.5F - ax);
else if (ix < 0x3f400000) /* |x| < 0.75 */
s = __kernel_cospif(ax - 0.5F);
else
s = __kernel_sinpif(1 - ax);
return ((hx & 0x80000000) ? -s : s);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
ax -= x;
GET_FLOAT_WORD(ix, ax);
if (ix == 0)
s = 0;
else {
if (ix < 0x3f000000) { /* |x| < 0.5 */
if (ix < 0x3e800000) /* |x| < 0.25 */
s = __kernel_sinpif(ax);
else
s = __kernel_cospif(0.5F - ax);
} else {
if (ix < 0x3f400000) /* |x| < 0.75 */
s = __kernel_cospif(ax - 0.5F);
else
s = __kernel_sinpif(1 - ax);
}
j0 = (uint32_t)x;
s = (j0 & 1) ? -s : s;
}
return ((hx & 0x80000000) ? -s : s);
}
/* x = +-inf or nan. */
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p23 is always an integer, so return +-0.
*/
return (copysignf(0, x));
}

176
lib/msun/src/s_tanpi.c Normal file
View File

@ -0,0 +1,176 @@
/*-
* Copyright (c) 2017 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.
*/
/**
* tanpi(x) computes tan(pi*x) without multiplication by pi (almost). First,
* note that tanpi(-x) = -tanpi(x), so the algorithm considers only |x| and
* includes reflection symmetry by considering the sign of x on output. The
* method used depends on the magnitude of x.
*
* 1. For small |x|, tanpi(x) = pi * x where a sloppy threshold is used. The
* threshold is |x| < 0x1pN with N = -(P/2+M). P is the precision of the
* floating-point type and M = 2 to 4. To achieve high accuracy, pi is
* decomposed into high and low parts with the high part containing a
* number of trailing zero bits. x is also split into high and low parts.
*
* 2. For |x| < 1, argument reduction is not required and tanpi(x) is
* computed by a direct call to a kernel, which uses the kernel for
* tan(x). See below.
*
* 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
* |x| = j0 + r with j0 an integer and the remainder r satisfies
* 0 <= r < 1. With the given domain, a simplified inline floor(x)
* is used. Also, note the following identity
*
* tan(pi*j0) + tan(pi*r)
* tanpi(x) = tan(pi*(j0+r)) = ---------------------------- = tanpi(r)
* 1 - tan(pi*j0) * tan(pi*r)
*
* So, after argument reduction, the kernel is again invoked.
*
* 4. For |x| >= 0x1p(P-1), |x| is integral and tanpi(x) = copysign(0,x).
*
* 5. Special cases:
*
* tanpi(+-0) = +-0
* tanpi(+-n) = +-0, for positive integers n.
* tanpi(+-n+1/4) = +-1, for positive integers n.
* tanpi(+-n+1/2) = NaN, for positive integers n.
* tanpi(+-inf) = NaN. Raises the "invalid" floating-point exception.
* tanpi(nan) = NaN. Raises the "invalid" floating-point exception.
*/
#include "math.h"
#include "math_private.h"
static const double
pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
pi_lo = -2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
/*
* The kernel for tanpi(x) multiplies x by an 80-bit approximation of
* pi, where the hi and lo parts are used with with kernel for tan(x).
*/
static inline double
__kernel_tanpi(double x)
{
double_t hi, lo, t;
if (x < 0.25) {
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = __kernel_tan(hi, lo, 1);
} else if (x > 0.25) {
x = 0.5 - x;
hi = (float)x;
lo = x - hi;
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
hi *= pi_hi;
_2sumF(hi, lo);
t = - __kernel_tan(hi, lo, -1);
} else
t = 1;
return (t);
}
volatile static const double vzero = 0;
double
tanpi(double x)
{
double ax, hi, lo, t;
uint32_t hx, ix, j0, lx;
EXTRACT_WORDS(hx, lx, x);
ix = hx & 0x7fffffff;
INSERT_WORDS(ax, ix, lx);
if (ix < 0x3ff00000) { /* |x| < 1 */
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix < 0x3e200000) { /* |x| < 0x1p-29 */
if (x == 0)
return (x);
/*
* To avoid issues with subnormal values,
* scale the computation and rescale on
* return.
*/
INSERT_WORDS(hi, hx, 0);
hi *= 0x1p53;
lo = x * 0x1p53 - hi;
t = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
return (t * 0x1p-53);
}
t = __kernel_tanpi(ax);
} else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
else
t = - __kernel_tanpi(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
if (ix < 0x43300000) { /* 1 <= |x| < 0x1p52 */
/* Determine integer part of ax. */
j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
ix &= ~(0x000fffff >> j0);
lx = 0;
} else {
lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
}
INSERT_WORDS(x,ix,lx);
ax -= x;
EXTRACT_WORDS(ix, lx, ax);
if (ix < 0x3fe00000) /* |x| < 0.5 */
t = ax == 0 ? 0 : __kernel_tanpi(ax);
else if (ax == 0.5)
return ((ax - ax) / (ax - ax));
else
t = - __kernel_tanpi(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
/* x = +-inf or nan. */
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p52 is always an integer, so return +-0.
*/
return (copysign(0, x));
}
#if LDBL_MANT_DIG == 53
__weak_reference(tanpi, tanpil);
#endif

114
lib/msun/src/s_tanpif.c Normal file
View File

@ -0,0 +1,114 @@
/*-
* Copyright (c) 2017 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.
*/
/*
* See ../src/s_tanpi.c for implementation details.
*/
#define INLINE_KERNEL_TANDF
#include "math.h"
#include "math_private.h"
#include "k_tanf.c"
static const float
pi_hi = 3.14160156e+00F, /* 0x40491000 */
pi_lo = -8.90890988e-06F; /* 0xb715777a */
static inline float
__kernel_tanpif(float x)
{
float t;
if (x < 0.25F)
t = __kernel_tandf(M_PI * x, 1);
else if (x > 0.25F)
t = -__kernel_tandf(M_PI * (0.5 - x), -1);
else
t = 1;
return (t);
}
volatile static const float vzero = 0;
float
tanpif(float x)
{
float ax, hi, lo, t;
uint32_t hx, ix, j0;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
SET_FLOAT_WORD(ax, ix);
if (ix < 0x3f800000) { /* |x| < 1 */
if (ix < 0x3f000000) { /* |x| < 0.5 */
if (ix < 0x38800000) { /* |x| < 0x1p-14 */
if (ix == 0)
return (x);
SET_FLOAT_WORD(hi, hx & 0xffff0000);
hi *= 0x1p23F;
lo = x * 0x1p23F - hi;
t = (pi_lo + pi_hi) * lo + pi_lo * hi +
pi_hi * hi;
return (t * 0x1p-23F);
}
t = __kernel_tanpif(ax);
} else if (ix == 0x3f000000)
return ((ax - ax) / (ax - ax));
else
t = - __kernel_tanpif(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
if (ix < 0x4b000000) { /* 1 <= |x| < 0x1p23 */
/* Determine integer part of ax. */
j0 = ((ix >> 23) & 0xff) - 0x7f;
ix &= ~(0x007fffff >> j0);
SET_FLOAT_WORD(x, ix);
ax -= x;
GET_FLOAT_WORD(ix, ax);
if (ix < 0x3f000000) /* |x| < 0.5 */
t = ix == 0 ? 0 : __kernel_tanpif(ax);
else if (ix == 0x3f000000)
return ((ax - ax) / (ax - ax));
else
t = - __kernel_tanpif(1 - ax);
return ((hx & 0x80000000) ? -t : t);
}
/* x = +-inf or nan. */
if (ix >= 0x7f800000)
return (vzero / vzero);
/*
* |x| >= 0x1p23 is always an integer, so return +-0.
*/
return (copysignf(0, x));
}