Attachment #229037 for bug #218514

Lines 126-131 COMMON_SRCS+= catrigl.c \ Link Here

(-)b/lib/msun/Makefile (-2 / +12 lines)
126	# See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=130067	126	# See also: https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=130067
127	.if ${COMPILER_TYPE} == "gcc"	127	.if ${COMPILER_TYPE} == "gcc"
128	CFLAGS.e_powl.c+= -Wno-error=overflow	128	CFLAGS.e_powl.c+= -Wno-error=overflow
		129
		130	# IEEE-754 2008 and ISO/IEC TS 18661-4 half-cycle trignometric functions
		131	COMMON_SRCS+= s_cospi.c s_cospif.c s_cospil.c \
		132	s_sinpi.c s_sinpif.c s_sinpil.c \
		133	s_tanpi.c s_tanpif.c s_tanpil.c
		134
129	.endif	135	.endif
130	.endif	136	.endif
131		137
Lines 154-160 INCS+= fenv.h math.h Link Here
154		160
155	MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \	161	MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
156	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \	162	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
157	cimag.3 clog.3 copysign.3 cos.3 cosh.3 cpow.3 csqrt.3 erf.3 \	163	cimag.3 clog.3 copysign.3 cos.3 cosh.3 cospi.3 \
		164	cpow.3 csqrt.3 erf.3 \
158	exp.3 fabs.3 fdim.3 \	165	exp.3 fabs.3 fdim.3 \
159	feclearexcept.3 feenableexcept.3 fegetenv.3 \	166	feclearexcept.3 feenableexcept.3 fegetenv.3 \
160	fegetround.3 fenv.3 floor.3 \	167	fegetround.3 fenv.3 floor.3 \
Lines 162-168 MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \ Link Here
162	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \	169	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
163	nextafter.3 remainder.3 rint.3 \	170	nextafter.3 remainder.3 rint.3 \
164	round.3 scalbn.3 signbit.3 sin.3 sincos.3 \	171	round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
165	sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \	172	sinh.3 sinpi.3 sqrt.3 tan.3 tanh.3 tanpi.3 trunc.3 \
166	complex.3	173	complex.3
167		174
168	MLINKS+=acos.3 acosf.3 acos.3 acosl.3	175	MLINKS+=acos.3 acosf.3 acos.3 acosl.3
Lines 192-197 MLINKS+=clog.3 clogf.3 clog.3 clogl.3 Link Here
192	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3	199	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
193	MLINKS+=cos.3 cosf.3 cos.3 cosl.3	200	MLINKS+=cos.3 cosf.3 cos.3 cosl.3
194	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3	201	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
		202	MLINKS+=cospi.3 cospif.3 cospi.3 cospil.3
195	MLINKS+=cpow.3 cpowf.3 cpow.3 cpowl.3	203	MLINKS+=cpow.3 cpowf.3 cpow.3 cpowl.3
196	MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3	204	MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
197	MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3	205	MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
Lines 244-253 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3 Link Here
244	MLINKS+=sin.3 sinf.3 sin.3 sinl.3	252	MLINKS+=sin.3 sinf.3 sin.3 sinl.3
245	MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3	253	MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
246	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3	254	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
		255	MLINKS+=sinpi.3 sinpif.3 sinpi.3 sinpil.3
247	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \	256	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
248	sqrt.3 sqrtl.3	257	sqrt.3 sqrtl.3
249	MLINKS+=tan.3 tanf.3 tan.3 tanl.3	258	MLINKS+=tan.3 tanf.3 tan.3 tanl.3
250	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3	259	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
		260	MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3
251	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3	261	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
252		262
253	.include <src.opts.mk>	263	.include <src.opts.mk>




	sincosf;
	sincosl;
};

/* First added in 14.0-CURRENT */
FBSD_1.7 {
	cospi;
	cospif;
	cospil;
	sinpi;
	sinpif;
	sinpil;
	tanpi;
	tanpif;
	tanpil;
};




/*-
 * 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));
}




/*-
 * 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));
}




/*-
 * 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);
}




/*-
 * 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));
}




/*-
 * 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));
}




/*-
 * 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));
}




/*-
 * 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));
}




/*-
 * 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);
}




/*-
 * 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));
}




/*-
 * 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));
}




.\" 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.






.\" 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.





.\" 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.






/*-
 * 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));
}





/*-
 * 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));
}

Lines 508-513 long double lgammal_r(long double, int *); Link Here

(-)b/lib/msun/src/math.h (+9 lines)
508	void sincos(double, double , double );	508	void sincos(double, double , double );
509	void sincosf(float, float , float );	509	void sincosf(float, float , float );
510	void sincosl(long double, long double , long double );	510	void sincosl(long double, long double , long double );
		511	double cospi(double);
		512	float cospif(float);
		513	long double cospil(long double);
		514	double sinpi(double);
		515	float sinpif(float);
		516	long double sinpil(long double);
		517	double tanpi(double);
		518	float tanpif(float);
		519	long double tanpil(long double);
511	#endif	520	#endif
512		521
513	__END_DECLS	522	__END_DECLS




/*-
 * 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




/*-
 * 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);
}




/*-
 * 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




/*-
 * 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));
}




/*-
 * 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




/*-
 * 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));
}

Return to bug 218514

Lines 304-306 FBSD_1.5 { Link Here

(-)b/lib/msun/Symbol.map (+13 lines)
304	sincosf;	304	sincosf;
305	sincosl;	305	sincosl;
306	};	306	};
		307
		308	/* First added in 14.0-CURRENT */
		309	FBSD_1.7 {
		310	cospi;
		311	cospif;
		312	cospil;
		313	sinpi;
		314	sinpif;
		315	sinpil;
		316	tanpi;
		317	tanpif;
		318	tanpil;
		319	};

Added Link Here

(-)b/lib/msun/ld128/k_cospil.h (+42 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/k_cospi.c for implementation details.
29	*/
30
31	static inline long double
32	__kernel_cospil(long double x)
33	{
34	long double hi, lo;
35
36	hi = (double)x;
37	lo = x - hi;
38	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
39	hi *= pi_hi;
40	_2sumF(hi, lo);
41	return (__kernel_cosl(hi, lo));
42	}

Added Link Here

(-)b/lib/msun/ld128/k_sinpil.h (+42 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/k_sinpi.c for implementation details.
29	*/
30
31	static inline long double
32	__kernel_sinpil(long double x)
33	{
34	long double hi, lo;
35
36	hi = (double)x;
37	lo = x - hi;
38	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
39	hi *= pi_hi;
40	_2sumF(hi, lo);
41	return (__kernel_sinl(hi, lo, 1));
42	}

Added Link Here

(-)b/lib/msun/ld128/s_cospil.c (+109 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_cospi.c for implementation details.
29	*
30	* FIXME: This has not been compiled nor has it been tested for accuracy.
31	* FIXME: This should use bit twiddling.
32	*/
33
34	#include "math.h"
35	#include "math_private.h"
36
37	/*
38	* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
39	*/
40	static const long double
41	pi_hi = 3.14159265358979322702026593105983920e+00L,
42	pi_lo = 1.14423774522196636802434264184180742e-17L;
43
44	#include "k_cospil.h"
45	#include "k_sinpil.h"
46
47	volatile static const double vzero = 0;
48
49	long double
50	cospil(long double x)
51	{
52	long double ax, c, xf;
53	uint32_t ix;
54
55	ax = fabsl(x);
56
57	if (ax < 1) {
58	if (ax < 0.25) {
59	if (ax < 0x1p-60) {
60	if ((int)x == 0)
61	return (1);
62	}
63	return (__kernel_cospil(ax));
64	}
65
66	if (ax < 0.5)
67	c = __kernel_sinpil(0.5 - ax);
68	else if (ax < 0.75) {
69	if (ax == 0.5)
70	return (0);
71	c = -__kernel_sinpil(ax - 0.5);
72	} else
73	c = -__kernel_cospil(1 - ax);
74	return (c);
75	}
76
77	if (ax < 0x1p112) {
78	xf = floorl(ax);
79	ax -= xf;
80	if (x < 0.5) {
81	if (x < 0.25)
82	c = ax == 0 ? 1 : __kernel_cospil(ax);
83	else
84	c = __kernel_sinpil(0.5 - ax);
85	} else {
86	if (x < 0.75) {
87	if (ax == 0.5)
88	return (0);
89	c = -__kernel_sinpil(ax - 0.5);
90	} else
91	c = -__kernel_cospil(1 - ax);
92	}
93
94	if (xf > 0x1p50)
95	xf -= 0x1p50;
96	if (xf > 0x1p30)
97	xf -= 0x1p30;
98	ix = (uint32_t)xf;
99	return (ix & 1 ? -c : c);
100	}
101
102	if (isinf(x) \|\| isnan(x))
103	return (vzero / vzero);
104
105	/*
106	* \|x\| >= 0x1p112 is always an even integer, so return 1.
107	*/
108	return (1);
109	}

Added Link Here

(-)b/lib/msun/ld128/s_sinpil.c (+118 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_sinpi.c for implementation details.
29	*
30	* FIXME: This has not been compiled nor has it been tested for accuracy.
31	* FIXME: This should use bit twiddling.
32	*/
33
34	#include "math.h"
35	#include "math_private.h"
36
37	/*
38	* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
39	*/
40	static const long double
41	pi_hi = 3.14159265358979322702026593105983920e+00L,
42	pi_lo = 1.14423774522196636802434264184180742e-17L;
43
44	#include "k_cospil.h"
45	#include "k_sinpil.h"
46
47	volatile static const double vzero = 0;
48
49	long double
50	sinpil(long double x)
51	{
52	long double ax, hi, lo, s, xf, xhi, xlo;
53	uint32_t ix;
54
55	ax = fabsl(x);
56
57	if (ax < 1) {
58	if (ax < 0.25) {
59	if (ax < 0x1p-60) {
60	if (x == 0)
61	return (x);
62	hi = (double)x;
63	hi *= 0x1p113L;
64	lo = x * 0x1p113L - hi;
65	s = (pi_lo + pi_hi) * lo + pi_lo * lo +
66	pi_hi * hi;
67	return (s * 0x1p-113L);
68	}
69
70	s = __kernel_sinpil(ax);
71	return (copysignl(s, x));
72	}
73
74	if (ax < 0.5)
75	s = __kernel_cospil(0.5 - ax);
76	else if (ax < 0.75)
77	s = __kernel_cospil(ax - 0.5);
78	else
79	s = __kernel_sinpil(1 - ax);
80	return (copysignl(s, x));
81	}
82
83	if (ax < 0x1p112) {
84	xf = floorl(ax);
85	ax -= xf;
86	if (ax == 0) {
87	s = 0;
88	} else {
89	if (ax < 0.5) {
90	if (ax <= 0.25)
91	s = __kernel_sinpil(ax);
92	else
93	s = __kernel_cospil(0.5 - ax);
94	} else {
95	if (ax < 0.75)
96	s = __kernel_cospil(ax - 0.5);
97	else
98	s = __kernel_sinpil(1 - ax);
99	}
100
101	if (xf > 0x1p50)
102	xf -= 0x1p50;
103	if (xf > 0x1p30)
104	xf -= 0x1p30;
105	ix = (uint32_t)xf;
106	if (ix & 1) s = -s;
107	}
108	return (copysignl(s, x));
109	}
110
111	if (isinf(x) \|\| isnan(x))
112	return (vzero / vzero);
113
114	/*
115	* \|x\| >= 0x1p112 is always an integer, so return +-0.
116	*/
117	return (copysignl(0, x));
118	}

Added Link Here

(-)b/lib/msun/ld128/s_tanpil.c (+120 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_tanpi.c for implementation details.
29	*
30	* FIXME: This has not been compiled nor has it been tested for accuracy.
31	* FIXME: This should use bit twiddling.
32	*/
33
34	#include "math.h"
35	#include "math_private.h"
36
37	/*
38	* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
39	*/
40	static const long double
41	pi_hi = 3.14159265358979322702026593105983920e+00L,
42	pi_lo = 1.14423774522196636802434264184180742e-17L;
43
44	static inline long double
45	__kernel_tanpi(long double x)
46	{
47	long double hi, lo, t;
48
49	if (x < 0.25) {
50	hi = (double)x;
51	lo = x - hi;
52	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
53	hi *= pi_hi;
54	_2sumF(hi, lo);
55	t = __kernel_tanl(hi, lo, -1);
56	} else if (x > 0.25) {
57	x = 0.5 - x;
58	hi = (double)x;
59	lo = x - hi;
60	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
61	hi *= pi_hi;
62	_2sumF(hi, lo);
63	t = - __kernel_tanl(hi, lo, 1);
64	} else
65	t = 1;
66
67	return (t);
68	}
69
70	volatile static const double vzero = 0;
71
72	long double
73	tanpil(long double x)
74	{
75	long double ax, hi, lo, xf;
76	uint32_t ix;
77
78	ax = fabsl(ax);
79
80	if (ax < 1) {
81	if (ax < 0.5) {
82	if (ax < 0x1p-60) {
83	if (x == 0)
84	return (x);
85	hi = (double)x;
86	hi *= 0x1p113L
87	lo = x * 0x1p113L - hi;
88	t = (pi_lo + pi_hi) * lo + pi_lo * lo +
89	pi_hi * hi;
90	return (t * 0x1p-113L);
91	}
92	t = __kernel_tanpil(ax);
93	} else if (ax == 0.5)
94	return ((ax - ax) / (ax - ax));
95	else
96	t = -__kernel_tanpil(1 - ax);
97	return (copysignl(t, x));
98	}
99
100	if (ix < 0x1p112) {
101	xf = floorl(ax);
102	ax -= xf;
103	if (ax < 0.5)
104	t = ax == 0 ? 0 : __kernel_tanpil(ax);
105	else if (ax == 0.5)
106	return ((ax - ax) / (ax - ax));
107	else
108	t = -__kernel_tanpil(1 - ax);
109	return (copysignl(t, x));
110	}
111
112	/* x = +-inf or nan. */
113	if (isinf(x) \|\| isnan(x))
114	return (vzero / vzero);
115
116	/*
117	* \|x\| >= 0x1p53 is always an integer, so return +-0.
118	*/
119	return (copysignl(0, x));
120	}

Added Link Here

(-)b/lib/msun/ld80/k_cospil.h (+42 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/k_cospi.c for implementation details.
29	*/
30
31	static inline long double
32	__kernel_cospil(long double x)
33	{
34	long double hi, lo;
35
36	hi = (float)x;
37	lo = x - hi;
38	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
39	hi *= pi_hi;
40	_2sumF(hi, lo);
41	return (__kernel_cosl(hi, lo));
42	}

Added Link Here

(-)b/lib/msun/ld80/k_sinpil.h (+42 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/k_sinpi.c for implementation details.
29	*/
30
31	static inline long double
32	__kernel_sinpil(long double x)
33	{
34	long double hi, lo;
35
36	hi = (float)x;
37	lo = x - hi;
38	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
39	hi *= pi_hi;
40	_2sumF(hi, lo);
41	return (__kernel_sinl(hi, lo, 1));
42	}

Added Link Here

(-)b/lib/msun/ld80/s_cospil.c (+129 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_cospi.c for implementation details.
29	*/
30
31	#ifdef __i386__
32	#include <ieeefp.h>
33	#endif
34	#include <stdint.h>
35
36	#include "fpmath.h"
37	#include "math.h"
38	#include "math_private.h"
39
40	static const double
41	pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
42	pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
43
44	#include "k_cospil.h"
45	#include "k_sinpil.h"
46
47	volatile static const double vzero = 0;
48
49	long double
50	cospil(long double x)
51	{
52	long double ax, c;
53	uint64_t lx, m;
54	uint32_t j0;
55	uint16_t hx, ix;
56
57	EXTRACT_LDBL80_WORDS(hx, lx, x);
58	ix = hx & 0x7fff;
59	INSERT_LDBL80_WORDS(ax, ix, lx);
60
61	ENTERI();
62
63	if (ix < 0x3fff) { /* \|x\| < 1 */
64	if (ix < 0x3ffd) { /* \|x\| < 0.25 */
65	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
66	if ((int)x == 0)
67	RETURNI(1);
68	}
69	RETURNI(__kernel_cospil(ax));
70	}
71
72	if (ix < 0x3ffe) /* \|x\| < 0.5 */
73	c = __kernel_sinpil(0.5 - ax);
74	else if (lx < 0xc000000000000000ull) { /* \|x\| < 0.75 */
75	if (ax == 0.5)
76	RETURNI(0);
77	c = -__kernel_sinpil(ax - 0.5);
78	} else
79	c = -__kernel_cospil(1 - ax);
80	RETURNI(c);
81	}
82
83	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
84	/* Determine integer part of ax. */
85	j0 = ix - 0x3fff + 1;
86	if (j0 < 32) {
87	lx = (lx >> 32) << 32;
88	lx &= ~(((lx << 32)-1) >> j0);
89	} else {
90	m = (uint64_t)-1 >> (j0 + 1);
91	if (lx & m) lx &= ~m;
92	}
93	INSERT_LDBL80_WORDS(x, ix, lx);
94
95	ax -= x;
96	EXTRACT_LDBL80_WORDS(ix, lx, ax);
97
98	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
99	if (ix < 0x3ffd) /* \|x\| < 0.25 */
100	c = ix == 0 ? 1 : __kernel_cospil(ax);
101	else
102	c = __kernel_sinpil(0.5 - ax);
103
104	} else {
105	if (lx < 0xc000000000000000ull) { /* \|x\| < 0.75 */
106	if (ax == 0.5)
107	RETURNI(0);
108	c = -__kernel_sinpil(ax - 0.5);
109	} else
110	c = -__kernel_cospil(1 - ax);
111	}
112
113	if (j0 > 40)
114	x -= 0x1p40;
115	if (j0 > 30)
116	x -= 0x1p30;
117	j0 = (uint32_t)x;
118
119	RETURNI(j0 & 1 ? -c : c);
120	}
121
122	if (ix >= 0x7fff)
123	RETURNI(vzero / vzero);
124
125	/*
126	* \|x\| >= 0x1p63 is always an even integer, so return 1.
127	*/
128	RETURNI(1);
129	}

Added Link Here

(-)b/lib/msun/ld80/s_sinpil.c (+140 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_sinpi.c for implementation details.
29	*/
30
31	#ifdef __i386__
32	#include <ieeefp.h>
33	#endif
34	#include <stdint.h>
35
36	#include "fpmath.h"
37	#include "math.h"
38	#include "math_private.h"
39
40	static const union IEEEl2bits
41	pi_hi_u = LD80C(0xc90fdaa200000000, 1, 3.14159265346825122833e+00L),
42	pi_lo_u = LD80C(0x85a308d313198a2e, -33, 1.21542010130123852029e-10L);
43	#define pi_hi (pi_hi_u.e)
44	#define pi_lo (pi_lo_u.e)
45
46	#include "k_cospil.h"
47	#include "k_sinpil.h"
48
49	volatile static const double vzero = 0;
50
51	long double
52	sinpil(long double x)
53	{
54	long double ax, hi, lo, s;
55	uint64_t lx, m;
56	uint32_t j0;
57	uint16_t hx, ix;
58
59	EXTRACT_LDBL80_WORDS(hx, lx, x);
60	ix = hx & 0x7fff;
61	INSERT_LDBL80_WORDS(ax, ix, lx);
62
63	ENTERI();
64
65	if (ix < 0x3fff) { /* \|x\| < 1 */
66	if (ix < 0x3ffd) { /* \|x\| < 0.25 */
67	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
68	if (x == 0)
69	RETURNI(x);
70	INSERT_LDBL80_WORDS(hi, hx,
71	lx & 0xffffffff00000000ull);
72	hi *= 0x1p63L;
73	lo = x * 0x1p63L - hi;
74	s = (pi_lo + pi_hi) * lo + pi_lo * hi +
75	pi_hi * hi;
76	RETURNI(s * 0x1p-63L);
77	}
78	s = __kernel_sinpil(ax);
79	RETURNI((hx & 0x8000) ? -s : s);
80	}
81
82	if (ix < 0x3ffe) /* \|x\| < 0.5 */
83	s = __kernel_cospil(0.5 - ax);
84	else if (lx < 0xc000000000000000ull) /* \|x\| < 0.75 */
85	s = __kernel_cospil(ax - 0.5);
86	else
87	s = __kernel_sinpil(1 - ax);
88	RETURNI((hx & 0x8000) ? -s : s);
89	}
90
91	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
92	/* Determine integer part of ax. */
93	j0 = ix - 0x3fff + 1;
94	if (j0 < 32) {
95	lx = (lx >> 32) << 32;
96	lx &= ~(((lx << 32)-1) >> j0);
97	} else {
98	m = (uint64_t)-1 >> (j0 + 1);
99	if (lx & m) lx &= ~m;
100	}
101	INSERT_LDBL80_WORDS(x, ix, lx);
102
103	ax -= x;
104	EXTRACT_LDBL80_WORDS(ix, lx, ax);
105
106	if (ix == 0) {
107	s = 0;
108	} else {
109	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
110	if (ix < 0x3ffd) /* \|x\| < 0.25 */
111	s = __kernel_sinpil(ax);
112	else
113	s = __kernel_cospil(0.5 - ax);
114	} else {
115	/* \|x\| < 0.75 */
116	if (lx < 0xc000000000000000ull)
117	s = __kernel_cospil(ax - 0.5);
118	else
119	s = __kernel_sinpil(1 - ax);
120	}
121
122	if (j0 > 40)
123	x -= 0x1p40;
124	if (j0 > 30)
125	x -= 0x1p30;
126	j0 = (uint32_t)x;
127	if (j0 & 1) s = -s;
128	}
129	RETURNI((hx & 0x8000) ? -s : s);
130	}
131
132	/* x = +-inf or nan. */
133	if (ix >= 0x7fff)
134	RETURNI(vzero / vzero);
135
136	/*
137	* \|x\| >= 0x1p63 is always an integer, so return +-0.
138	*/
139	RETURNI(copysignl(0, x));
140	}

Added Link Here

(-)b/lib/msun/ld80/s_tanpil.c (+139 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_tanpi.c for implementation details.
29	*/
30
31	#ifdef __i386__
32	#include <ieeefp.h>
33	#endif
34	#include <stdint.h>
35
36	#include "fpmath.h"
37	#include "math.h"
38	#include "math_private.h"
39
40	static const double
41	pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
42	pi_lo = -2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
43
44	static inline long double
45	__kernel_tanpil(long double x)
46	{
47	long double hi, lo, t;
48
49	if (x < 0.25) {
50	hi = (float)x;
51	lo = x - hi;
52	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
53	hi *= pi_hi;
54	_2sumF(hi, lo);
55	t = __kernel_tanl(hi, lo, -1);
56	} else if (x > 0.25) {
57	x = 0.5 - x;
58	hi = (float)x;
59	lo = x - hi;
60	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
61	hi *= pi_hi;
62	_2sumF(hi, lo);
63	t = - __kernel_tanl(hi, lo, 1);
64	} else
65	t = 1;
66
67	return (t);
68	}
69
70	volatile static const double vzero = 0;
71
72	long double
73	tanpil(long double x)
74	{
75	long double ax, hi, lo, t;
76	uint64_t lx, m;
77	uint32_t j0;
78	uint16_t hx, ix;
79
80	EXTRACT_LDBL80_WORDS(hx, lx, x);
81	ix = hx & 0x7fff;
82	INSERT_LDBL80_WORDS(ax, ix, lx);
83
84	ENTERI();
85
86	if (ix < 0x3fff) { /* \|x\| < 1 */
87	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
88	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
89	if (x == 0)
90	RETURNI(x);
91	INSERT_LDBL80_WORDS(hi, hx,
92	lx & 0xffffffff00000000ull);
93	hi *= 0x1p63L;
94	lo = x * 0x1p63L - hi;
95	t = (pi_lo + pi_hi) * lo + pi_lo * hi +
96	pi_hi * hi;
97	RETURNI(t * 0x1p-63L);
98	}
99	t = __kernel_tanpil(ax);
100	} else if (ax == 0.5)
101	RETURNI((ax - ax) / (ax - ax));
102	else
103	t = -__kernel_tanpil(1 - ax);
104	RETURNI((hx & 0x8000) ? -t : t);
105	}
106
107	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
108	/* Determine integer part of ax. */
109	j0 = ix - 0x3fff + 1;
110	if (j0 < 32) {
111	lx = (lx >> 32) << 32;
112	lx &= ~(((lx << 32)-1) >> j0);
113	} else {
114	m = (uint64_t)-1 >> (j0 + 1);
115	if (lx & m) lx &= ~m;
116	}
117	INSERT_LDBL80_WORDS(x, ix, lx);
118
119	ax -= x;
120	EXTRACT_LDBL80_WORDS(ix, lx, ax);
121
122	if (ix < 0x3ffe) /* \|x\| < 0.5 */
123	t = ax == 0 ? 0 : __kernel_tanpil(ax);
124	else if (ax == 0.5)
125	RETURNI((ax - ax) / (ax - ax));
126	else
127	t = -__kernel_tanpil(1 - ax);
128	RETURNI((hx & 0x8000) ? -t : t);
129	}
130
131	/* x = +-inf or nan. */
132	if (ix >= 0x7fff)
133	RETURNI(vzero / vzero);
134
135	/*
136	* \|x\| >= 0x1p63 is always an integer, so return +-0.
137	*/
138	RETURNI(copysignl(0, x));
139	}

Added Link Here

(-)b/lib/msun/man/cospi.3 (+111 lines)
1	.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
2	.\" All rights reserved.
3	.\"
4	.\" Redistribution and use in source and binary forms, with or without
5	.\" modification, are permitted provided that the following conditions
6	.\" are met:
7	.\" 1. Redistributions of source code must retain the above copyright
8	.\" notice, this list of conditions and the following disclaimer.
9	.\" 2. Redistributions in binary form must reproduce the above copyright
10	.\" notice, this list of conditions and the following disclaimer in the
11	.\" documentation and/or other materials provided with the distribution.
12	.\"
13	.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
14	.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15	.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
16	.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
17	.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
18	.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
19	.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
20	.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
21	.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
22	.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
23	.\" SUCH DAMAGE.
24	.\"
25	.\" $FreeBSD$
26	.\"
27	.Dd April 1, 2017
28	.Dt COSPI 3
29	.Os
30	.Sh NAME
31	.Nm cospi ,
32	.Nm cospif ,
33	.Nm cospil
34	.Nd half\(encycle cosine functions
35	.Sh LIBRARY
36	.Lb libm
37	.Sh SYNOPSIS
38	.In math.h
39	.Ft double
40	.Fn cospi "double x"
41	.Ft float
42	.Fn cospif "float x"
43	.Ft long double
44	.Fn cospil "long double x"
45	.Sh DESCRIPTION
46	The
47	.Fn cospi ,
48	.Fn cospif ,
49	and
50	.Fn cospil
51	functions compute the cosine of
52	.Fa "\(*p \(mu x" .
53	and measure angles in half-cycles.
54	.Sh RETURN VALUES
55	The
56	.Fn cospi ,
57	.Fn cospif ,
58	and
59	.Fn cospil
60	functions returns
61	.Fn cos "\(*p \(mu x" .
62	If \(Bax\(Ba \*(Ge 2^(p - 1)
63	where p is the floating\(enpoint precision of
64	.Ar x ,
65	then the returned value is 1 and it has no significance.
66	.Sh SPECIAL VALUES
67	.Bl -tag
68	.It
69	.Fn cospi \*(Pm0
70	returns 1.
71	.It
72	.Fn cospi \*(Pmn/2
73	returns 0 for positive integers
74	.Ar n .
75	.It
76	.Fn cospi n
77	returns 1 for even integers
78	.Ar n .
79	.It
80	.Fn cospi n
81	returns \-1 for odd integers
82	.Ar n .
83	.It
84	.Fn cospi \*(Pm\(if
85	return an \*(Na and raises an FE_INVALID exception.
86	.It
87	.Fn cospi \*(Na
88	return an \*(Na and raises an FE_INVALID exception.
89	.El
90	.Sh SEE ALSO
91	.Xr cos 3 ,
92	.Xr fenv 3 ,
93	.Xr math 3 ,
94	.Xr sin 3 ,
95	.Xr sinpi 3 ,
96	.Xr tan 3 ,
97	.Xr tanpi 3
98	.Sh AUTHORS
99	The half\(encycle trignometric functions were written by
100	.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
101	.Sh STANDARDS
102	These functions conform to
103	IEEE Std 754\(tm\(en2008 ,
104	\(dqIEEE Standard for Floating-Point Arithmetic\(dq
105	and to
106	ISO/IEC TS 18661-4 ,
107	\(dqInformation technology \(em Programming languages, their environments,
108	and system software interfaces \(em Floating\(enpoint extensions for
109	C\(dq \(em Part 4: Supplementary functions.
110
111

Added Link Here

(-)b/lib/msun/man/sinpi.3 (+102 lines)
1	.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
2	.\" All rights reserved.
3	.\"
4	.\" Redistribution and use in source and binary forms, with or without
5	.\" modification, are permitted provided that the following conditions
6	.\" are met:
7	.\" 1. Redistributions of source code must retain the above copyright
8	.\" notice, this list of conditions and the following disclaimer.
9	.\" 2. Redistributions in binary form must reproduce the above copyright
10	.\" notice, this list of conditions and the following disclaimer in the
11	.\" documentation and/or other materials provided with the distribution.
12	.\"
13	.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
14	.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15	.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
16	.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
17	.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
18	.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
19	.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
20	.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
21	.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
22	.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
23	.\" SUCH DAMAGE.
24	.\"
25	.\" $FreeBSD$
26	.\"
27	.Dd April 1, 2017
28	.Dt SINPI 3
29	.Os
30	.Sh NAME
31	.Nm sinpi ,
32	.Nm sinpif ,
33	.Nm sinpil
34	.Nd half\(encycle sine functions
35	.Sh LIBRARY
36	.Lb libm
37	.Sh SYNOPSIS
38	.In math.h
39	.Ft double
40	.Fn sinpi "double x"
41	.Ft float
42	.Fn sinpif "float x"
43	.Ft long double
44	.Fn sinpil "long double x"
45	.Sh DESCRIPTION
46	The
47	.Fn sinpi ,
48	.Fn sinpif ,
49	and
50	.Fn sinpil
51	functions compute the sine of
52	.Fa "\(*p \(mu x" .
53	and measure angles in half-cycles.
54	.Sh RETURN VALUES
55	The
56	.Fn sinpi ,
57	.Fn sinpif ,
58	and
59	.Fn sinpil
60	functions returns
61	.Fn sin "\(*p \(mu x" .
62	If \(Bax\(Ba \*(Ge 2^(p - 1)
63	where p is the floating\(enpoint precision of
64	.Ar x ,
65	then the returned value is \*(Pm0 and it has no significance.
66	.Sh SPECIAL VALUES
67	.Bl -tag
68	.It
69	.Fn sinpi \*(Pm0
70	returns \*(Pm0.
71	.It
72	.Fn sinpi \*(Pmn
73	returns \*(Pm0 for positive integers
74	.Ar n .
75	.It
76	.Fn sinpi \*(Pm\(if
77	return an \*(Na and raises an FE_INVALID exception.
78	.It
79	.Fn sinpi \*(Na
80	return an \*(Na and raises an FE_INVALID exception.
81	.El
82	.Sh SEE ALSO
83	.Xr cos 3 ,
84	.Xr cospi 3 ,
85	.Xr fenv 3 ,
86	.Xr math 3 ,
87	.Xr sin 3 ,
88	.Xr tan 3 ,
89	.Xr tanpi 3
90	.Sh AUTHORS
91	The half\(encycle trignometric functions were written by
92	.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
93	.Sh STANDARDS
94	These functions conform to
95	IEEE Std 754\(tm\(en2008 ,
96	\(dqIEEE Standard for Floating-Point Arithmetic\(dq
97	and to
98	ISO/IEC TS 18661-4 ,
99	\(dqInformation technology \(em Programming languages, their environments,
100	and system software interfaces \(em Floating\(enpoint extensions for
101	C\(dq \(em Part 4: Supplementary functions.
102

Added Link Here

(-)b/lib/msun/man/tanpi.3 (+106 lines)
1	.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
2	.\" All rights reserved.
3	.\"
4	.\" Redistribution and use in source and binary forms, with or without
5	.\" modification, are permitted provided that the following conditions
6	.\" are met:
7	.\" 1. Redistributions of source code must retain the above copyright
8	.\" notice, this list of conditions and the following disclaimer.
9	.\" 2. Redistributions in binary form must reproduce the above copyright
10	.\" notice, this list of conditions and the following disclaimer in the
11	.\" documentation and/or other materials provided with the distribution.
12	.\"
13	.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
14	.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15	.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
16	.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
17	.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
18	.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
19	.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
20	.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
21	.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
22	.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
23	.\" SUCH DAMAGE.
24	.\"
25	.\" $FreeBSD$
26	.\"
27	.Dd April 1, 2017
28	.Dt TANPI 3
29	.Os
30	.Sh NAME
31	.Nm tanpi ,
32	.Nm tanpif ,
33	.Nm tanpil
34	.Nd half\(encycle tangent functions
35	.Sh LIBRARY
36	.Lb libm
37	.Sh SYNOPSIS
38	.In math.h
39	.Ft double
40	.Fn tanpi "double x"
41	.Ft float
42	.Fn tanpif "float x"
43	.Ft long double
44	.Fn tanpil "long double x"
45	.Sh DESCRIPTION
46	The
47	.Fn tanpi ,
48	.Fn tanpif ,
49	and
50	.Fn tanpil
51	functions compute the tangent of
52	.Fa "\(*p \(mu x"
53	and measure angles in half-cycles.
54	.Sh RETURN VALUES
55	The
56	.Fn tanpi ,
57	.Fn tanpif ,
58	and
59	.Fn tanpil
60	functions returns
61	.Fn tan "\(*p \(mu x" .
62	If \(Bax\(Ba \*(Ge 2^(p - 1)
63	where p is the floating\(enpoint precision of
64	.Ar x ,
65	then the returned value is \*(Pm0 and it has no significance.
66	.Sh SPECIAL VALUES
67	.Bl -tag
68	.It
69	.Fn tanpi \*(Pm0
70	returns \*(Pm0.
71	.It
72	.Fn tanpi \*(Pmn
73	returns \*(Pm0 for positive integers
74	.Ar n .
75	.It
76	.Fn tanpi \*(Pmn/2
77	returns \*(Na for n > 0 and raises an FE_INVALID exception.
78	.It
79	.Fn tanpi \*(Pm\(if
80	return an \*(Na and raises an FE_INVALID exception.
81	.It
82	.Fn tanpi \*(Na
83	return an \*(Na and raises an FE_INVALID exception.
84	.El
85	.Sh SEE ALSO
86	.Xr cos 3 ,
87	.Xr cospi 3 ,
88	.Xr fenv 3 ,
89	.Xr math 3 ,
90	.Xr sin 3 ,
91	.Xr sinpi 3 ,
92	.Xr tan 3 ,
93	.Sh AUTHORS
94	The half\(encycle trignometric functions were written by
95	.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
96	.Sh STANDARDS
97	These functions conform to
98	IEEE Std 754\(tm\(en2008 ,
99	\(dqIEEE Standard for Floating-Point Arithmetic\(dq
100	and to
101	ISO/IEC TS 18661-4 ,
102	\(dqInformation technology \(em Programming languages, their environments,
103	and system software interfaces \(em Floating\(enpoint extensions for
104	C\(dq \(em Part 4: Supplementary functions.
105
106

Added Link Here

(-)b/lib/msun/src/k_cospi.h (+44 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* The basic kernel for x in [0,0.25]. To use the kernel for cos(x), the
29	* argument to __kernel_cospi() must be multiplied by pi.
30	*/
31
32	static inline double
33	__kernel_cospi(double x)
34	{
35	double_t hi, lo;
36
37	hi = (float)x;
38	lo = x - hi;
39	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
40	hi *= pi_hi;
41	_2sumF(hi, lo);
42	return (__kernel_cos(hi, lo));
43	}
44

Added Link Here

(-)b/lib/msun/src/k_sinpi.h (+43 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* The basic kernel for x in [0,0.25]. To use the kernel for sin(x), the
29	* argument to __kernel_sinpi() must be multiplied by pi.
30	*/
31
32	static inline double
33	__kernel_sinpi(double x)
34	{
35	double_t hi, lo;
36
37	hi = (float)x;
38	lo = x - hi;
39	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
40	hi *= pi_hi;
41	_2sumF(hi, lo);
42	return (__kernel_sin(hi, lo, 1));
43	}

Added Link Here

(-)b/lib/msun/src/s_cospi.c (+151 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/**
28	* cospi(x) computes cos(pi*x) without multiplication by pi (almost). First,
29	* note that cospi(-x) = cospi(x), so the algorithm considers only \|x\|. The
30	* method used depends on the magnitude of x.
31	*
32	* 1. For small \|x\|, cospi(x) = 1 with FE_INEXACT raised where a sloppy
33	* threshold is used. The threshold is \|x\| < 0x1pN with N = -(P/2+M).
34	* P is the precision of the floating-point type and M = 2 to 4.
35	*
36	* 2. For \|x\| < 1, argument reduction is not required and sinpi(x) is
37	* computed by calling a kernel that leverages the kernels for sin(x)
38	* ans cos(x). See k_sinpi.c and k_cospi.c for details.
39	*
40	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
41	* \|x\| = j0 + r with j0 an integer and the remainder r satisfies
42	* 0 <= r < 1. With the given domain, a simplified inline floor(x)
43	* is used. Also, note the following identity
44	*
45	* cospi(x) = cos(pi*(j0+r))
46	* = cos(pij0) cos(pir) - sin(pij0) * sin(pi*r)
47	* = cos(pij0) cos(pi*r)
48	* = +-cospi(r)
49	*
50	* If j0 is even, then cos(pij0) = 1. If j0 is odd, then cos(pij0) = -1.
51	* cospi(r) is then computed via an appropriate kernel.
52	*
53	* 4. For \|x\| >= 0x1p(P-1), \|x\| is integral and cospi(x) = 1.
54	*
55	* 5. Special cases:
56	*
57	* cospi(+-0) = 1.
58	* cospi(n.5) = 0 for n an integer.
59	* cospi(+-inf) = nan. Raises the "invalid" floating-point exception.
60	* cospi(nan) = nan. Raises the "invalid" floating-point exception.
61	*/
62
63	#include "math.h"
64	#include "math_private.h"
65
66	static const double
67	pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
68	pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
69
70	#include "k_cospi.h"
71	#include "k_sinpi.h"
72
73	volatile static const double vzero = 0;
74
75	double
76	cospi(double x)
77	{
78	double ax, c;
79	uint32_t hx, ix, j0, lx;
80
81	EXTRACT_WORDS(hx, lx, x);
82	ix = hx & 0x7fffffff;
83	INSERT_WORDS(ax, ix, lx);
84
85	if (ix < 0x3ff00000) { /* \|x\| < 1 */
86	if (ix < 0x3fd00000) { /* \|x\| < 0.25 */
87	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
88	if ((int)ax == 0)
89	return (1);
90	}
91	return (__kernel_cospi(ax));
92	}
93
94	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
95	c = __kernel_sinpi(0.5 - ax);
96	else if (ix < 0x3fe80000){ /* \|x\| < 0.75 */
97	if (ax == 0.5)
98	return (0);
99	c = -__kernel_sinpi(ax - 0.5);
100	} else
101	c = -__kernel_cospi(1 - ax);
102	return (c);
103	}
104
105	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
106	/* Determine integer part of ax. */
107	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
108	if (j0 < 20) {
109	ix &= ~(0x000fffff >> j0);
110	lx = 0;
111	} else {
112	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
113	}
114	INSERT_WORDS(x, ix, lx);
115
116	ax -= x;
117	EXTRACT_WORDS(ix, lx, ax);
118
119
120	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
121	if (ix < 0x3fd00000) /* \|x\| < 0.25 */
122	c = ix == 0 ? 1 : __kernel_cospi(ax);
123	else
124	c = __kernel_sinpi(0.5 - ax);
125	} else {
126	if (ix < 0x3fe80000) { /* \|x\| < 0.75 */
127	if (ax == 0.5)
128	return (0);
129	c = -__kernel_sinpi(ax - 0.5);
130	} else
131	c = -__kernel_cospi(1 - ax);
132	}
133
134	if (j0 > 30)
135	x -= 0x1p30;
136	j0 = (uint32_t)x;
137	return (j0 & 1 ? -c : c);
138	}
139
140	if (ix >= 0x7f800000)
141	return (vzero / vzero);
142
143	/*
144	* \|x\| >= 0x1p52 is always an even integer, so return 1.
145	*/
146	return (1);
147	}
148
149	#if LDBL_MANT_DIG == 53
150	__weak_reference(cospi, cospil);
151	#endif

Added Link Here

(-)b/lib/msun/src/s_cospif.c (+109 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_cospi.c for implementation details.
29	*/
30	#define INLINE_KERNEL_SINDF
31	#define INLINE_KERNEL_COSDF
32
33	#include "math.h"
34	#include "math_private.h"
35	#include "k_cosf.c"
36	#include "k_sinf.c"
37
38	#define __kernel_cospif(x) (__kernel_cosdf(M_PI * (x)))
39	#define __kernel_sinpif(x) (__kernel_sindf(M_PI * (x)))
40
41	volatile static const float vzero = 0;
42
43	float
44	cospif(float x)
45	{
46	float ax, c;
47	uint32_t ix, j0;
48
49	GET_FLOAT_WORD(ix, x);
50	ix = ix & 0x7fffffff;
51	SET_FLOAT_WORD(ax, ix);
52
53	if (ix < 0x3f800000) { /* \|x\| < 1 */
54	if (ix < 0x3e800000) { /* \|x\| < 0.25 */
55	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
56	/* Raise inexact iff != 0. */
57	if ((int)ax == 0)
58	return (1);
59	}
60	return (__kernel_cospif(ax));
61	}
62
63	if (ix < 0x3f000000) /* \|x\| < 0.5 */
64	c = __kernel_sinpif(0.5F - ax);
65	else if (ix < 0x3f400000) { /* \|x\| < 0.75 */
66	if (ix == 0x3f000000)
67	return (0);
68	c = -__kernel_sinpif(ax - 0.5F);
69	} else
70	c = -__kernel_cospif(1 - ax);
71	return (c);
72	}
73
74	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
75	/* Determine integer part of ax. */
76	j0 = ((ix >> 23) & 0xff) - 0x7f;
77	ix &= ~(0x007fffff >> j0);
78	SET_FLOAT_WORD(x, ix);
79
80	ax -= x;
81	GET_FLOAT_WORD(ix, ax);
82
83	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
84	if (ix < 0x3e800000) /* \|x\| < 0.25 */
85	c = ix == 0 ? 1 : __kernel_cospif(ax);
86	else
87	c = __kernel_sinpif(0.5F - ax);
88	} else {
89	if (ix < 0x3f400000) { /* \|x\| < 0.75 */
90	if (ix == 0x3f000000)
91	return (0);
92	c = -__kernel_sinpif(ax - 0.5F);
93	} else
94	c = -__kernel_cospif(1 - ax);
95	}
96
97	j0 = (uint32_t)x;
98	return (j0 & 1 ? -c : c);
99	}
100
101	/* x = +-inf or nan. */
102	if (ix >= 0x7f800000)
103	return (vzero / vzero);
104
105	/*
106	* \|x\| >= 0x1p23 is always an even integer, so return 1.
107	*/
108	return (1);
109	}

Added Link Here

(-)b/lib/msun/src/s_sinpi.c (+168 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/**
28	* sinpi(x) computes sin(pi*x) without multiplication by pi (almost). First,
29	* note that sinpi(-x) = -sinpi(x), so the algorithm considers only \|x\| and
30	* includes reflection symmetry by considering the sign of x on output. The
31	* method used depends on the magnitude of x.
32	*
33	* 1. For small \|x\|, sinpi(x) = pi * x where a sloppy threshold is used. The
34	* threshold is \|x\| < 0x1pN with N = -(P/2+M). P is the precision of the
35	* floating-point type and M = 2 to 4. To achieve high accuracy, pi is
36	* decomposed into high and low parts with the high part containing a
37	* number of trailing zero bits. x is also split into high and low parts.
38	*
39	* 2. For \|x\| < 1, argument reduction is not required and sinpi(x) is
40	* computed by calling a kernel that leverages the kernels for sin(x)
41	* ans cos(x). See k_sinpi.c and k_cospi.c for details.
42	*
43	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
44	* \|x\| = j0 + r with j0 an integer and the remainder r satisfies
45	* 0 <= r < 1. With the given domain, a simplified inline floor(x)
46	* is used. Also, note the following identity
47	*
48	* sinpi(x) = sin(pi*(j0+r))
49	* = sin(pij0) cos(pir) + cos(pij0) * sin(pi*r)
50	* = cos(pij0) sin(pi*r)
51	* = +-sinpi(r)
52	*
53	* If j0 is even, then cos(pij0) = 1. If j0 is odd, then cos(pij0) = -1.
54	* sinpi(r) is then computed via an appropriate kernel.
55	*
56	* 4. For \|x\| >= 0x1p(P-1), \|x\| is integral and sinpi(x) = copysign(0,x).
57	*
58	* 5. Special cases:
59	*
60	* sinpi(+-0) = +-0
61	* sinpi(+-n) = +-0, for positive integers n.
62	* sinpi(+-inf) = nan. Raises the "invalid" floating-point exception.
63	* sinpi(nan) = nan. Raises the "invalid" floating-point exception.
64	*/
65
66	#include "math.h"
67	#include "math_private.h"
68
69	static const double
70	pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
71	pi_lo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
72
73	#include "k_cospi.h"
74	#include "k_sinpi.h"
75
76	volatile static const double vzero = 0;
77
78	double
79	sinpi(double x)
80	{
81	double ax, hi, lo, s;
82	uint32_t hx, ix, j0, lx;
83
84	EXTRACT_WORDS(hx, lx, x);
85	ix = hx & 0x7fffffff;
86	INSERT_WORDS(ax, ix, lx);
87
88	if (ix < 0x3ff00000) { /* \|x\| < 1 */
89	if (ix < 0x3fd00000) { /* \|x\| < 0.25 */
90	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
91	if (x == 0)
92	return (x);
93	/*
94	* To avoid issues with subnormal values,
95	* scale the computation and rescale on
96	* return.
97	*/
98	INSERT_WORDS(hi, hx, 0);
99	hi *= 0x1p53;
100	lo = x * 0x1p53 - hi;
101	s = (pi_lo + pi_hi) * lo + pi_lo * hi +
102	pi_hi * hi;
103	return (s * 0x1p-53);
104	}
105
106	s = __kernel_sinpi(ax);
107	return ((hx & 0x80000000) ? -s : s);
108	}
109
110	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
111	s = __kernel_cospi(0.5 - ax);
112	else if (ix < 0x3fe80000) /* \|x\| < 0.75 */
113	s = __kernel_cospi(ax - 0.5);
114	else
115	s = __kernel_sinpi(1 - ax);
116	return ((hx & 0x80000000) ? -s : s);
117	}
118
119	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
120	/* Determine integer part of ax. */
121	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
122	if (j0 < 20) {
123	ix &= ~(0x000fffff >> j0);
124	lx = 0;
125	} else {
126	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
127	}
128	INSERT_WORDS(x, ix, lx);
129
130	ax -= x;
131	EXTRACT_WORDS(ix, lx, ax);
132
133	if (ix == 0)
134	s = 0;
135	else {
136	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
137	if (ix < 0x3fd00000) /* \|x\| < 0.25 */
138	s = __kernel_sinpi(ax);
139	else
140	s = __kernel_cospi(0.5 - ax);
141	} else {
142	if (ix < 0x3fe80000) /* \|x\| < 0.75 */
143	s = __kernel_cospi(ax - 0.5);
144	else
145	s = __kernel_sinpi(1 - ax);
146	}
147
148	if (j0 > 30)
149	x -= 0x1p30;
150	j0 = (uint32_t)x;
151	if (j0 & 1) s = -s;
152	}
153
154	return ((hx & 0x80000000) ? -s : s);
155	}
156
157	if (ix >= 0x7f800000)
158	return (vzero / vzero);
159
160	/*
161	* \|x\| >= 0x1p52 is always an integer, so return +-0.
162	*/
163	return (copysign(0, x));
164	}
165
166	#if LDBL_MANT_DIG == 53
167	__weak_reference(sinpi, sinpil);
168	#endif

Added Link Here

(-)b/lib/msun/src/s_sinpif.c (+122 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_sinpi.c for implementation details.
29	*/
30
31	#define INLINE_KERNEL_SINDF
32	#define INLINE_KERNEL_COSDF
33
34	#include "math.h"
35	#include "math_private.h"
36	#include "k_cosf.c"
37	#include "k_sinf.c"
38
39	#define __kernel_cospif(x) (__kernel_cosdf(M_PI * (x)))
40	#define __kernel_sinpif(x) (__kernel_sindf(M_PI * (x)))
41
42	static const float
43	pi_hi = 3.14160156e+00F, /* 0x40491000 */
44	pi_lo = -8.90890988e-06F; /* 0xb715777a */
45
46	volatile static const float vzero = 0;
47
48	float
49	sinpif(float x)
50	{
51	float ax, hi, lo, s;
52	uint32_t hx, ix, j0;
53
54	GET_FLOAT_WORD(hx, x);
55	ix = hx & 0x7fffffff;
56	SET_FLOAT_WORD(ax, ix);
57
58	if (ix < 0x3f800000) { /* \|x\| < 1 */
59	if (ix < 0x3e800000) { /* \|x\| < 0.25 */
60	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
61	if (x == 0)
62	return (x);
63	SET_FLOAT_WORD(hi, hx & 0xffff0000);
64	hi *= 0x1p23F;
65	lo = x * 0x1p23F - hi;
66	s = (pi_lo + pi_hi) * lo + pi_lo * hi +
67	pi_hi * hi;
68	return (s * 0x1p-23F);
69	}
70
71	s = __kernel_sinpif(ax);
72	return ((hx & 0x80000000) ? -s : s);
73	}
74
75	if (ix < 0x3f000000) /* \|x\| < 0.5 */
76	s = __kernel_cospif(0.5F - ax);
77	else if (ix < 0x3f400000) /* \|x\| < 0.75 */
78	s = __kernel_cospif(ax - 0.5F);
79	else
80	s = __kernel_sinpif(1 - ax);
81	return ((hx & 0x80000000) ? -s : s);
82	}
83
84	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
85	/* Determine integer part of ax. */
86	j0 = ((ix >> 23) & 0xff) - 0x7f;
87	ix &= ~(0x007fffff >> j0);
88	SET_FLOAT_WORD(x, ix);
89
90	ax -= x;
91	GET_FLOAT_WORD(ix, ax);
92
93	if (ix == 0)
94	s = 0;
95	else {
96	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
97	if (ix < 0x3e800000) /* \|x\| < 0.25 */
98	s = __kernel_sinpif(ax);
99	else
100	s = __kernel_cospif(0.5F - ax);
101	} else {
102	if (ix < 0x3f400000) /* \|x\| < 0.75 */
103	s = __kernel_cospif(ax - 0.5F);
104	else
105	s = __kernel_sinpif(1 - ax);
106	}
107
108	j0 = (uint32_t)x;
109	s = (j0 & 1) ? -s : s;
110	}
111	return ((hx & 0x80000000) ? -s : s);
112	}
113
114	/* x = +-inf or nan. */
115	if (ix >= 0x7f800000)
116	return (vzero / vzero);
117
118	/*
119	* \|x\| >= 0x1p23 is always an integer, so return +-0.
120	*/
121	return (copysignf(0, x));
122	}

Added Link Here

(-)b/lib/msun/src/s_tanpi.c (+176 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/**
28	* tanpi(x) computes tan(pi*x) without multiplication by pi (almost). First,
29	* note that tanpi(-x) = -tanpi(x), so the algorithm considers only \|x\| and
30	* includes reflection symmetry by considering the sign of x on output. The
31	* method used depends on the magnitude of x.
32	*
33	* 1. For small \|x\|, tanpi(x) = pi * x where a sloppy threshold is used. The
34	* threshold is \|x\| < 0x1pN with N = -(P/2+M). P is the precision of the
35	* floating-point type and M = 2 to 4. To achieve high accuracy, pi is
36	* decomposed into high and low parts with the high part containing a
37	* number of trailing zero bits. x is also split into high and low parts.
38	*
39	* 2. For \|x\| < 1, argument reduction is not required and tanpi(x) is
40	* computed by a direct call to a kernel, which uses the kernel for
41	* tan(x). See below.
42	*
43	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
44	* \|x\| = j0 + r with j0 an integer and the remainder r satisfies
45	* 0 <= r < 1. With the given domain, a simplified inline floor(x)
46	* is used. Also, note the following identity
47	*
48	* tan(pij0) + tan(pir)
49	* tanpi(x) = tan(pi*(j0+r)) = ---------------------------- = tanpi(r)
50	* 1 - tan(pij0) tan(pi*r)
51	*
52	* So, after argument reduction, the kernel is again invoked.
53	*
54	* 4. For \|x\| >= 0x1p(P-1), \|x\| is integral and tanpi(x) = copysign(0,x).
55	*
56	* 5. Special cases:
57	*
58	* tanpi(+-0) = +-0
59	* tanpi(+-n) = +-0, for positive integers n.
60	* tanpi(+-n+1/4) = +-1, for positive integers n.
61	* tanpi(+-n+1/2) = NaN, for positive integers n.
62	* tanpi(+-inf) = NaN. Raises the "invalid" floating-point exception.
63	* tanpi(nan) = NaN. Raises the "invalid" floating-point exception.
64	*/
65
66	#include "math.h"
67	#include "math_private.h"
68
69	static const double
70	pi_hi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
71	pi_lo = -2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
72
73	/*
74	* The kernel for tanpi(x) multiplies x by an 80-bit approximation of
75	* pi, where the hi and lo parts are used with with kernel for tan(x).
76	*/
77	static inline double
78	__kernel_tanpi(double x)
79	{
80	double_t hi, lo, t;
81
82	if (x < 0.25) {
83	hi = (float)x;
84	lo = x - hi;
85	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
86	hi *= pi_hi;
87	_2sumF(hi, lo);
88	t = __kernel_tan(hi, lo, 1);
89	} else if (x > 0.25) {
90	x = 0.5 - x;
91	hi = (float)x;
92	lo = x - hi;
93	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
94	hi *= pi_hi;
95	_2sumF(hi, lo);
96	t = - __kernel_tan(hi, lo, -1);
97	} else
98	t = 1;
99
100	return (t);
101	}
102
103	volatile static const double vzero = 0;
104
105	double
106	tanpi(double x)
107	{
108	double ax, hi, lo, t;
109	uint32_t hx, ix, j0, lx;
110
111	EXTRACT_WORDS(hx, lx, x);
112	ix = hx & 0x7fffffff;
113	INSERT_WORDS(ax, ix, lx);
114
115	if (ix < 0x3ff00000) { /* \|x\| < 1 */
116	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
117	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
118	if (x == 0)
119	return (x);
120	/*
121	* To avoid issues with subnormal values,
122	* scale the computation and rescale on
123	* return.
124	*/
125	INSERT_WORDS(hi, hx, 0);
126	hi *= 0x1p53;
127	lo = x * 0x1p53 - hi;
128	t = (pi_lo + pi_hi) * lo + pi_lo * hi +
129	pi_hi * hi;
130	return (t * 0x1p-53);
131	}
132	t = __kernel_tanpi(ax);
133	} else if (ax == 0.5)
134	return ((ax - ax) / (ax - ax));
135	else
136	t = - __kernel_tanpi(1 - ax);
137	return ((hx & 0x80000000) ? -t : t);
138	}
139
140	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
141	/* Determine integer part of ax. */
142	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
143	if (j0 < 20) {
144	ix &= ~(0x000fffff >> j0);
145	lx = 0;
146	} else {
147	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
148	}
149	INSERT_WORDS(x,ix,lx);
150
151	ax -= x;
152	EXTRACT_WORDS(ix, lx, ax);
153
154	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
155	t = ax == 0 ? 0 : __kernel_tanpi(ax);
156	else if (ax == 0.5)
157	return ((ax - ax) / (ax - ax));
158	else
159	t = - __kernel_tanpi(1 - ax);
160
161	return ((hx & 0x80000000) ? -t : t);
162	}
163
164	/* x = +-inf or nan. */
165	if (ix >= 0x7f800000)
166	return (vzero / vzero);
167
168	/*
169	* \|x\| >= 0x1p52 is always an integer, so return +-0.
170	*/
171	return (copysign(0, x));
172	}
173
174	#if LDBL_MANT_DIG == 53
175	__weak_reference(tanpi, tanpil);
176	#endif

Added Link Here

(-)b/lib/msun/src/s_tanpif.c (+114 lines)
1	/*-
2	* Copyright (c) 2017 Steven G. Kargl
3	* All rights reserved.
4	*
5	* Redistribution and use in source and binary forms, with or without
6	* modification, are permitted provided that the following conditions
7	* are met:
8	* 1. Redistributions of source code must retain the above copyright
9	* notice unmodified, this list of conditions, and the following
10	* disclaimer.
11	* 2. Redistributions in binary form must reproduce the above copyright
12	* notice, this list of conditions and the following disclaimer in the
13	* documentation and/or other materials provided with the distribution.
14	*
15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25	*/
26
27	/*
28	* See ../src/s_tanpi.c for implementation details.
29	*/
30
31	#define INLINE_KERNEL_TANDF
32
33	#include "math.h"
34	#include "math_private.h"
35	#include "k_tanf.c"
36
37	static const float
38	pi_hi = 3.14160156e+00F, /* 0x40491000 */
39	pi_lo = -8.90890988e-06F; /* 0xb715777a */
40
41	static inline float
42	__kernel_tanpif(float x)
43	{
44	float t;
45
46	if (x < 0.25F)
47	t = __kernel_tandf(M_PI * x, 1);
48	else if (x > 0.25F)
49	t = -__kernel_tandf(M_PI * (0.5 - x), -1);
50	else
51	t = 1;
52
53	return (t);
54	}
55
56	volatile static const float vzero = 0;
57
58	float
59	tanpif(float x)
60	{
61	float ax, hi, lo, t;
62	uint32_t hx, ix, j0;
63
64	GET_FLOAT_WORD(hx, x);
65	ix = hx & 0x7fffffff;
66	SET_FLOAT_WORD(ax, ix);
67
68	if (ix < 0x3f800000) { /* \|x\| < 1 */
69	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
70	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
71	if (ix == 0)
72	return (x);
73	SET_FLOAT_WORD(hi, hx & 0xffff0000);
74	hi *= 0x1p23F;
75	lo = x * 0x1p23F - hi;
76	t = (pi_lo + pi_hi) * lo + pi_lo * hi +
77	pi_hi * hi;
78	return (t * 0x1p-23F);
79	}
80	t = __kernel_tanpif(ax);
81	} else if (ix == 0x3f000000)
82	return ((ax - ax) / (ax - ax));
83	else
84	t = - __kernel_tanpif(1 - ax);
85	return ((hx & 0x80000000) ? -t : t);
86	}
87
88	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
89	/* Determine integer part of ax. */
90	j0 = ((ix >> 23) & 0xff) - 0x7f;
91	ix &= ~(0x007fffff >> j0);
92	SET_FLOAT_WORD(x, ix);
93
94	ax -= x;
95	GET_FLOAT_WORD(ix, ax);
96
97	if (ix < 0x3f000000) /* \|x\| < 0.5 */
98	t = ix == 0 ? 0 : __kernel_tanpif(ax);
99	else if (ix == 0x3f000000)
100	return ((ax - ax) / (ax - ax));
101	else
102	t = - __kernel_tanpif(1 - ax);
103	return ((hx & 0x80000000) ? -t : t);
104	}
105
106	/* x = +-inf or nan. */
107	if (ix >= 0x7f800000)
108	return (vzero / vzero);
109
110	/*
111	* \|x\| >= 0x1p23 is always an integer, so return +-0.
112	*/
113	return (copysignf(0, x));
114	}