Attachment #182648 for bug #218514

Lines 59-64 Link Here

(-)lib/msun/Makefile (-5 / +23 lines)
59	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \	59	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
60	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \	60	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \
61	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \	61	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
		62	s_cospi.c s_cospif.c \
62	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \	63	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
63	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \	64	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
64	s_finite.c s_finitef.c \	65	s_finite.c s_finitef.c \
Lines 74-80 Link Here
74	s_rint.c s_rintf.c s_round.c s_roundf.c \	75	s_rint.c s_rintf.c s_round.c s_roundf.c \
75	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \	76	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
76	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \	77	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
77	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \	78	s_sinpi.c s_sinpif.c \
		79	s_tan.c s_tanf.c s_tanh.c s_tanhf.c \
		80	s_tanpi.c s_tanpif.c \
		81	s_tgammaf.c s_trunc.c s_truncf.c \
78	w_cabs.c w_cabsf.c w_drem.c w_dremf.c	82	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
79		83
80	# Location of fpmath.h and _fpmath.h	84	# Location of fpmath.h and _fpmath.h
Lines 100-110 Link Here
100	e_lgammal.c e_lgammal_r.c \	104	e_lgammal.c e_lgammal_r.c \
101	e_remainderl.c e_sinhl.c e_sqrtl.c \	105	e_remainderl.c e_sinhl.c e_sqrtl.c \
102	invtrig.c k_cosl.c k_sinl.c k_tanl.c \	106	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
103	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \	107	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c \
		108	s_cospil.c \
		109	s_cprojl.c \
104	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \	110	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
105	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \	111	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
106	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \	112	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
107	s_scalbnl.c s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c	113	s_scalbnl.c s_sinl.c \
		114	s_sinpil.c \
		115	s_tanhl.c s_tanl.c \
		116	s_tanpil.c s_truncl.c w_cabsl.c
108	.endif	117	.endif
109		118
110	# C99 complex functions	119	# C99 complex functions
Lines 131-143 Link Here
131		140
132	MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \	141	MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
133	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \	142	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
134	cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \	143	cimag.3 copysign.3 cos.3 cosh.3 \
		144	cospi.3 \
		145	csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
135	feclearexcept.3 feenableexcept.3 fegetenv.3 \	146	feclearexcept.3 feenableexcept.3 fegetenv.3 \
136	fegetround.3 fenv.3 floor.3 \	147	fegetround.3 fenv.3 floor.3 \
137	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \	148	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
138	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \	149	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
139	nextafter.3 remainder.3 rint.3 \	150	nextafter.3 remainder.3 rint.3 \
140	round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \	151	round.3 scalbn.3 signbit.3 sin.3 sinh.3 \
		152	sinpi.3 \
		153	sqrt.3 tan.3 tanh.3 \
		154	tanpi.3 \
		155	trunc.3 \
141	complex.3	156	complex.3
142		157
143	MLINKS+=acos.3 acosf.3 acos.3 acosl.3	158	MLINKS+=acos.3 acosf.3 acos.3 acosl.3
Lines 166-171 Link Here
166	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3	181	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
167	MLINKS+=cos.3 cosf.3 cos.3 cosl.3	182	MLINKS+=cos.3 cosf.3 cos.3 cosl.3
168	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3	183	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
		184	MLINKS+=cospi.3 cospif.3 cospi.3 cospil.3
169	MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3	185	MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
170	MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3	186	MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
171	MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \	187	MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \
Lines 216-225 Link Here
216	MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3	232	MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
217	MLINKS+=sin.3 sinf.3 sin.3 sinl.3	233	MLINKS+=sin.3 sinf.3 sin.3 sinl.3
218	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3	234	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
		235	MLINKS+=sinpi.3 sinpif.3 sinpi.3 sinpil.3
219	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \	236	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
220	sqrt.3 sqrtl.3	237	sqrt.3 sqrtl.3
221	MLINKS+=tan.3 tanf.3 tan.3 tanl.3	238	MLINKS+=tan.3 tanf.3 tan.3 tanl.3
222	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3	239	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
		240	MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3
223	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3	241	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
224		242
225	.include <src.opts.mk>	243	.include <src.opts.mk>

Lines 294-297 Link Here

(-)lib/msun/Symbol.map (+9 lines)
294	casinl;	294	casinl;
295	catanl;	295	catanl;
296	catanhl;	296	catanhl;
		297	cospi;
		298	cospif;
		299	cospil;
		300	sinpi;
		301	sinpif;
		302	sinpil;
		303	tanpi;
		304	tanpif;
		305	tanpil;
297	};	306	};




/*-
 * 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 * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_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 * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_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"

static const long double
pihi = 3.14159265358979322702026593105983920e+00L,
pilo = 1.14423774522196636802434264184180742e-17L;

#include "k_cospil.c"
#include "k_sinpil.c"

long double
cospil(long double x)
{
	volatile static const double vzero = 0;
	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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = 1;
	return (ax);
}




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

static const long double
pihi = 3.14159265358979322702026593105983920e+00L,
pilo = 1.14423774522196636802434264184180742e-17L;

#include "k_cospil.c"
#include "k_sinpil.c"

long double
sinpil(long double x)
{
	volatile static const double vzero = 0;
	long double ax, 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);
				ax = (double)x;
				x -= ax;
				s = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (s);
			}

			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = copysignl(0, x);
	return (ax);
}




/*-
 * 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.
 * FIXME: This should use a polynomial approximation in [0,0.25], but the
 * FIXME: increase in max ULP is probably small, so who cares.
 */

#include "math.h"
#include "math_private.h"

static const long double
pihi = 3.14159265358979322702026593105983920e+00L,
pilo = 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 * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_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 * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_2sumF(hi, lo);
		t = - __kernel_tanl(hi, lo, 1);
	} else
		t = 1;

	return (t);
}

long double
tanpil(long double x)
{
	volatile static const double vzero = 0;
	long double ax, xf;
	uint32_t ix;

	ax = fabsl(ax);

	if (ax < 1) {
		if (ax < 0.5) {
			if (ax < 0x1p-60) {
				if (x == 0)
					return (x);
				ax = (double)x;
				x -= ax;
				t = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (t);
			}
			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = copysignl(0, x);
	return (ax);
}




/*-
 * 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 * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_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 * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_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
pihi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */

#include "k_cospil.c"
#include "k_sinpil.c"

long double
cospil(long double x)
{
	volatile static const double vzero = 0;
	long double ax, c;
	uint32_t j0;
	uint64_t lx;
	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 {
			uint64_t 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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = 1;
	RETURNI(ax);
}




/*-
 * 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 double
pihi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */

#include "k_cospil.c"
#include "k_sinpil.c"

long double
sinpil(long double x)
{
	volatile static const double vzero = 0;
	long double ax, s;
	uint32_t j0;
	uint64_t lx;
	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);
				ax = (float)x;
				x -= ax;
				s = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				RETURNI(s);
			}
			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 {
			uint64_t 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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = copysignl(0, x);
	RETURNI(ax);
}




/*-
 * 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
pihi =  3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo = -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 * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_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 * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_2sumF(hi, lo);
		t = - __kernel_tanl(hi, lo, 1);
	} else
		t = 1;

	return (t);
}

long double
tanpil(long double x)
{
	volatile static const double vzero = 0;
	long double ax, t;
	uint32_t j0;
	uint64_t lx;
	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);
				ax = (float)x;
				x -= ax;
				t = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				RETURNI(t);
			}
			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 {
			uint64_t 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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = copysignl(0, x);
	RETURNI(ax);
}




.\" 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 exploit the kernel for cos(x),
 * the argument to __kernel_cospi() must be multiply by pi.  As pi is an
 * irrational number, this allows the splitting of pi*x into high and low
 * parts.
 */

static inline double
__kernel_cospi(double x)
{
	double hi, lo;
	hi = (float)x;
	lo = x - hi;
	lo = lo * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_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 exploit the kernel for sin(x),
 * the argument to __kernel_sinpi() must be multiply by pi.  As pi is an
 * irrational number, this allows the splitting of pi*x into high and low
 * parts.
 */

static inline double
__kernel_sinpi(double x)
{
	double hi, lo;
	hi = (float)x;
	lo = x - hi;
	lo = lo * (pilo + pihi) + hi * pilo;
	hi *= pihi;
	_2sumF(hi, lo);
	return (__kernel_sin(hi, lo, 1));
}

Lines 500-505 Link Here

(-)lib/msun/src/math.h (+9 lines)
500		500
501	#if __BSD_VISIBLE	501	#if __BSD_VISIBLE
502	long double lgammal_r(long double, int *);	502	long double lgammal_r(long double, int *);
		503	double cospi(double);
		504	float cospif(float);
		505	long double cospil(long double);
		506	double sinpi(double);
		507	float sinpif(float);
		508	long double sinpil(long double);
		509	double tanpi(double);
		510	float tanpif(float);
		511	long double tanpil(long double);
503	#endif	512	#endif
504		513
505	__END_DECLS	514	__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
pihi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */

#include "k_cospi.c"
#include "k_sinpi.c"

double
cospi(double x)
{
	volatile static const double vzero = 0;
	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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID?
	 */
	if (ax + 1 > 1)
		ax = 1;
	return (ax);
}

#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)   ((float)__kernel_cosdf(M_PI * (x)))
#define  __kernel_sinpif(x)   ((float)__kernel_sindf(M_PI * (x)))

float
cospif(float x)
{
	volatile static const float vzero = 0;
	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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = 1;
	return (ax);
}




/*-
 * 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
pihi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */

#include "k_cospi.c"
#include "k_sinpi.c"

double
sinpi(double x)
{
	volatile static const double vzero = 0;
	double ax, 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);
				INSERT_WORDS(ax, hx, 0);
				x -= ax;
				s = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (s);
			}

			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID?
	 */
	if (ax + 1 > 1)
		ax = copysign(0, x);
	return (ax);
}

#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)   ((float)__kernel_cosdf(M_PI * (x)))
#define  __kernel_sinpif(x)   ((float)__kernel_sindf(M_PI * (x)))

static const float pihi =  3.14160156e+00f;	/* 0x40491000 */
static const float pilo = -8.90890988e-06f;	/* 0xb715777a */

float
sinpif(float x)
{
	volatile static const float vzero = 0;
	float ax, 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(ax, hx & 0xffff0000);
				x -= ax;
				s = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (s);
			}

			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID?
	 */
	if (ax + 1 > 1)
		ax = copysignf(0, x);
	return (ax);
}




/*-
 * 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 
pihi  =  3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
pilo  = -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 hi, lo, t;

	if (x < 0.25) {
		hi = (float)x;
		lo = x - hi;
		lo = lo * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_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 * (pilo + pihi) + hi * pilo;
		hi *= pihi;
		_2sumF(hi, lo);
		t = - __kernel_tan(hi, lo, -1);
	} else
		t = 1;

	return (t);
}

double
tanpi(double x)
{
	volatile static const double vzero = 0;
	double ax, 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);
				ax = (float)x;
				x -= ax;
				t = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (t);
			}
			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID?
	 */
	if (ax + 1 > 1)
		ax = copysign(0, x);
	return (ax);
}

#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
pihi =  3.14160156e+00f,	/* 0x40491000 */
pilo = -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);
}

float
tanpif(float x)
{
	volatile static const float vzero = 0;
	float ax, 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(ax, hx & 0xffff0000);
				x -= ax;
				t = pilo * x + pihi * x + pilo * ax
				    + pihi * ax;
				return (t);
			}
			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.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	if (ax + 1 > 1)
		ax = copysignf(0, x);
	return (ax);
}

Return to bug 218514

Line 0 Link Here

(-)lib/msun/ld128/k_cospil.c (+41 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	hi = (double)x;
	36	lo = x - hi;
	37	lo = lo * (pilo + pihi) + hi * pilo;
	38	hi *= pihi;
	39	_2sumF(hi, lo);
	40	return (__kernel_cosl(hi, lo));
	41	}

Line 0 Link Here

(-)lib/msun/ld128/k_sinpil.c (+41 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	hi = (double)x;
	36	lo = x - hi;
	37	lo = lo * (pilo + pihi) + hi * pilo;
	38	hi *= pihi;
	39	_2sumF(hi, lo);
	40	return (__kernel_sinl(hi, lo, 1));
	41	}

Line 0 Link Here

(-)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	static const long double
		38	pihi = 3.14159265358979322702026593105983920e+00L,
		39	pilo = 1.14423774522196636802434264184180742e-17L;
		40
		41	#include "k_cospil.c"
		42	#include "k_sinpil.c"
		43
		44	long double
		45	cospil(long double x)
		46	{
		47	volatile static const double vzero = 0;
		48	long double ax, c, xf;
		49	uint32_t ix;
		50
		51	ax = fabsl(x);
		52
		53	if (ax < 1) {
		54	if (ax < 0.25) {
		55	if (ax < 0x1p-60) {
		56	if ((int)x == 0)
		57	return (1);
		58	}
		59	return (__kernel_cospil(ax));
		60	}
		61
		62	if (ax < 0.5)
		63	c = __kernel_sinpil(0.5 - ax);
		64	else if (ax < 0.75) {
65	if (ax == 0.5)
66	return (0);
67	c = -__kernel_sinpil(ax - 0.5);
68	} else
69	c = -__kernel_cospil(1 - ax);
70	return (c);
71	}
72
73	if (ax < 0x1p112) {
74
75	xf = floorl(ax);
76
77	ax -= xf;
78
79	if (x < 0.5) {
80	if (x < 0.25)
81	c = ax == 0 ? 1 : __kernel_cospil(ax);
82	else
83	c = __kernel_sinpil(0.5 - ax);
84	} else {
85	if (x < 0.75) {
86	if (ax == 0.5)
87	return (0);
88	c = -__kernel_sinpil(ax - 0.5);
89	} else
90	c = -__kernel_cospil(1 - ax);
91	}
92
93	if (xf > 0x1p50) xf -= 0x1p50;
94	if (xf > 0x1p30) xf -= 0x1p30;
95	ix = (uint32_t)xf;
96	return (ix & 1 ? -c : c);
97	}
98
99	if (isinf(x) \|\| isnan(x))
100	return (vzero / vzero);
101
102	/*
103	* \|x\| >= 0x1p112 is always an even integer, so return 1.
104	* FIXME: should this raise FE_INEXACT or FE_INVALID.
105	*/
106	if (ax + 1 > 1)
107	ax = 1;
108	return (ax);
109	}

Line 0 Link Here

(-)lib/msun/ld128/s_sinpil.c (+117 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	static const long double
		38	pihi = 3.14159265358979322702026593105983920e+00L,
		39	pilo = 1.14423774522196636802434264184180742e-17L;
		40
		41	#include "k_cospil.c"
		42	#include "k_sinpil.c"
		43
		44	long double
		45	sinpil(long double x)
		46	{
		47	volatile static const double vzero = 0;
		48	long double ax, s, xf, xhi, xlo;
		49	uint32_t ix;
		50
		51	ax = fabsl(x);
		52
		53	if (ax < 1) {
		54	if (ax < 0.25) {
		55	if (ax < 0x1p-60) {
		56	if (x == 0)
		57	return (x);
		58	ax = (double)x;
		59	x -= ax;
		60	s = pilo * x + pihi * x + pilo * ax
		61	+ pihi * ax;
		62	return (s);
		63	}
		64
65	s = __kernel_sinpil(ax);
66	return (copysignl(s, x));
67	}
68
69	if (ax < 0.5)
70	s = __kernel_cospil(0.5 - ax);
71	else if (ax < 0.75)
72	s = __kernel_cospil(ax - 0.5);
73	else
74	s = __kernel_sinpil(1 - ax);
75	return (copysignl(s, x));
76	}
77
78	if (ax < 0x1p112) {
79
80	xf = floorl(ax);
81
82	ax -= xf;
83
84	if (ax == 0) {
85	s = 0;
86	} else {
87	if (ax < 0.5) {
88	if (ax <= 0.25)
89	s = __kernel_sinpil(ax);
90	else
91	s = __kernel_cospil(0.5 - ax);
92	} else {
93	if (ax < 0.75)
94	s = __kernel_cospil(ax - 0.5);
95	else
96	s = __kernel_sinpil(1 - ax);
97	}
98
99	if (xf > 0x1p50) xf -= 0x1p50;
100	if (xf > 0x1p30) xf -= 0x1p30;
101	ix = (uint32_t)xf;
102	if (ix & 1) s = -s;
103	}
104	return (copysignl(s, x));
105	}
106
107	if (isinf(x) \|\| isnan(x))
108	return (vzero / vzero);
109
110	/*
111	* \|x\| >= 0x1p112 is always an integer, so return +-0.
112	* FIXME: should this raise FE_INEXACT or FE_INVALID.
113	*/
114	if (ax + 1 > 1)
115	ax = copysignl(0, x);
116	return (ax);
117	}

Line 0 Link Here

(-)lib/msun/ld128/s_tanpil.c (+123 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	* FIXME: This should use a polynomial approximation in [0,0.25], but the
		33	* FIXME: increase in max ULP is probably small, so who cares.
		34	*/
		35
		36	#include "math.h"
		37	#include "math_private.h"
		38
		39	static const long double
		40	pihi = 3.14159265358979322702026593105983920e+00L,
		41	pilo = 1.14423774522196636802434264184180742e-17L;
		42
		43	static inline long double
		44	__kernel_tanpi(long double x)
		45	{
		46	long double hi, lo, t;
		47
		48	if (x < 0.25) {
		49	hi = (double)x;
		50	lo = x - hi;
		51	lo = lo * (pilo + pihi) + hi * pilo;
		52	hi *= pihi;
		53	_2sumF(hi, lo);
		54	t = __kernel_tanl(hi, lo, -1);
		55	} else if (x > 0.25) {
		56	x = 0.5 - x;
		57	hi = (double)x;
		58	lo = x - hi;
		59	lo = lo * (pilo + pihi) + hi * pilo;
		60	hi *= pihi;
		61	_2sumF(hi, lo);
		62	t = - __kernel_tanl(hi, lo, 1);
		63	} else
		64	t = 1;
65
66	return (t);
67	}
68
69	long double
70	tanpil(long double x)
71	{
72	volatile static const double vzero = 0;
73	long double ax, xf;
74	uint32_t ix;
75
76	ax = fabsl(ax);
77
78	if (ax < 1) {
79	if (ax < 0.5) {
80	if (ax < 0x1p-60) {
81	if (x == 0)
82	return (x);
83	ax = (double)x;
84	x -= ax;
85	t = pilo * x + pihi * x + pilo * ax
86	+ pihi * ax;
87	return (t);
88	}
89	t = __kernel_tanpil(ax);
90	} else if (ax == 0.5)
91	return ((ax - ax) / (ax - ax));
92	else
93	t = -__kernel_tanpil(1 - ax);
94	return (copysignl(t, x));
95	}
96
97	if (ix < 0x1p112) {
98
99	xf = floorl(ax);
100
101	ax -= xf;
102
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	* FIXME: should this raise FE_INEXACT or FE_INVALID.
119	*/
120	if (ax + 1 > 1)
121	ax = copysignl(0, x);
122	return (ax);
123	}

Line 0 Link Here

(-)lib/msun/ld80/k_cospil.c (+41 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	hi = (float)x;
	36	lo = x - hi;
	37	lo = lo * (pilo + pihi) + hi * pilo;
	38	hi *= pihi;
	39	_2sumF(hi, lo);
	40	return (__kernel_cosl(hi, lo));
	41	}

Line 0 Link Here

(-)lib/msun/ld80/k_sinpil.c (+41 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	hi = (float)x;
	36	lo = x - hi;
	37	lo = lo * (pilo + pihi) + hi * pilo;
	38	hi *= pihi;
	39	_2sumF(hi, lo);
	40	return (__kernel_sinl(hi, lo, 1));
	41	}

Line 0 Link Here

(-)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	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
		42	pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
		43
		44	#include "k_cospil.c"
		45	#include "k_sinpil.c"
		46
		47	long double
		48	cospil(long double x)
		49	{
		50	volatile static const double vzero = 0;
		51	long double ax, c;
		52	uint32_t j0;
		53	uint64_t lx;
		54	uint16_t hx, ix;
		55
		56	EXTRACT_LDBL80_WORDS(hx, lx, x);
		57	ix = hx & 0x7fff;
		58	INSERT_LDBL80_WORDS(ax, ix, lx);
		59
		60	ENTERI();
		61
		62	if (ix < 0x3fff) { /* \|x\| < 1 */
		63	if (ix < 0x3ffd) { /* \|x\| < 0.25 */
		64	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
65	if ((int)x == 0)
66	RETURNI(1);
67	}
68	RETURNI(__kernel_cospil(ax));
69	}
70
71	if (ix < 0x3ffe) /* \|x\| < 0.5 */
72	c = __kernel_sinpil(0.5 - ax);
73	else if (lx < 0xc000000000000000ull) { /* \|x\| < 0.75 */
74	if (ax == 0.5)
75	RETURNI(0);
76	c = -__kernel_sinpil(ax - 0.5);
77	} else
78	c = -__kernel_cospil(1 - ax);
79	RETURNI(c);
80	}
81
82	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
83	/* Determine integer part of ax. */
84	j0 = ix - 0x3fff + 1;
85	if (j0 < 32) {
86	lx = (lx >> 32) << 32;
87	lx &= ~(((lx << 32)-1) >> j0);
88	} else {
89	uint64_t m = (uint64_t)-1 >> (j0 + 1);
90	if (lx & m) lx &= ~m;
91	}
92	INSERT_LDBL80_WORDS(x, ix, lx);
93
94	ax -= x;
95	EXTRACT_LDBL80_WORDS(ix, lx, ax);
96
97	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
98	if (ix < 0x3ffd) /* \|x\| < 0.25 */
99	c = ix == 0 ? 1 : __kernel_cospil(ax);
100	else
101	c = __kernel_sinpil(0.5 - ax);
102
103	} else {
104	if (lx < 0xc000000000000000ull) { /* \|x\| < 0.75 */
105	if (ax == 0.5)
106	RETURNI(0);
107	c = -__kernel_sinpil(ax - 0.5);
108	} else
109	c = -__kernel_cospil(1 - ax);
110	}
111
112	if (j0 > 40) x -= 0x1p40;
113	if (j0 > 30) x -= 0x1p30;
114	j0 = (uint32_t)x;
115
116	RETURNI(j0 & 1 ? -c : c);
117	}
118
119	if (ix >= 0x7fff)
120	RETURNI(vzero / vzero);
121
122	/*
123	* \|x\| >= 0x1p63 is always an even integer, so return 1.
124	* FIXME: should this raise FE_INEXACT or FE_INVALID.
125	*/
126	if (ax + 1 > 1)
127	ax = 1;
128	RETURNI(ax);
129	}

Line 0 Link Here

(-)lib/msun/ld80/s_sinpil.c (+136 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 double
		41	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
		42	pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
		43
		44	#include "k_cospil.c"
		45	#include "k_sinpil.c"
		46
		47	long double
		48	sinpil(long double x)
		49	{
		50	volatile static const double vzero = 0;
		51	long double ax, s;
		52	uint32_t j0;
		53	uint64_t lx;
		54	uint16_t hx, ix;
		55
		56	EXTRACT_LDBL80_WORDS(hx, lx, x);
		57	ix = hx & 0x7fff;
		58	INSERT_LDBL80_WORDS(ax, ix, lx);
		59
		60	ENTERI();
		61
		62	if (ix < 0x3fff) { /* \|x\| < 1 */
		63	if (ix < 0x3ffd) { /* \|x\| < 0.25 */
		64	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
65	if (x == 0)
66	RETURNI(x);
67	ax = (float)x;
68	x -= ax;
69	s = pilo * x + pihi * x + pilo * ax
70	+ pihi * ax;
71	RETURNI(s);
72	}
73	s = __kernel_sinpil(ax);
74	RETURNI((hx & 0x8000) ? -s : s);
75	}
76
77	if (ix < 0x3ffe) /* \|x\| < 0.5 */
78	s = __kernel_cospil(0.5 - ax);
79	else if (lx < 0xc000000000000000ull) /* \|x\| < 0.75 */
80	s = __kernel_cospil(ax - 0.5);
81	else
82	s = __kernel_sinpil(1 - ax);
83	RETURNI((hx & 0x8000) ? -s : s);
84	}
85
86	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
87	/* Determine integer part of ax. */
88	j0 = ix - 0x3fff + 1;
89	if (j0 < 32) {
90	lx = (lx >> 32) << 32;
91	lx &= ~(((lx << 32)-1) >> j0);
92	} else {
93	uint64_t m = (uint64_t)-1 >> (j0 + 1);
94	if (lx & m) lx &= ~m;
95	}
96	INSERT_LDBL80_WORDS(x, ix, lx);
97
98	ax -= x;
99	EXTRACT_LDBL80_WORDS(ix, lx, ax);
100
101	if (ix == 0) {
102	s = 0;
103	} else {
104	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
105	if (ix < 0x3ffd) /* \|x\| < 0.25 */
106	s = __kernel_sinpil(ax);
107	else
108	s = __kernel_cospil(0.5 - ax);
109	} else {
110	/* \|x\| < 0.75 */
111	if (lx < 0xc000000000000000ull)
112	s = __kernel_cospil(ax - 0.5);
113	else
114	s = __kernel_sinpil(1 - ax);
115	}
116
117	if (j0 > 40) x -= 0x1p40;
118	if (j0 > 30) x -= 0x1p30;
119	j0 = (uint32_t)x;
120	if (j0 & 1) s = -s;
121	}
122	RETURNI((hx & 0x8000) ? -s : s);
123	}
124
125	/* x = +-inf or nan. */
126	if (ix >= 0x7fff)
127	RETURNI(vzero / vzero);
128
129	/*
130	* \|x\| >= 0x1p63 is always an integer, so return +-0.
131	* FIXME: should this raise FE_INEXACT or FE_INVALID.
132	*/
133	if (ax + 1 > 1)
134	ax = copysignl(0, x);
135	RETURNI(ax);
136	}

Line 0 Link Here

(-)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	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
		42	pilo = -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 * (pilo + pihi) + hi * pilo;
		53	hi *= pihi;
		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 * (pilo + pihi) + hi * pilo;
		61	hi *= pihi;
		62	_2sumF(hi, lo);
		63	t = - __kernel_tanl(hi, lo, 1);
		64	} else
65	t = 1;
66
67	return (t);
68	}
69
70	long double
71	tanpil(long double x)
72	{
73	volatile static const double vzero = 0;
74	long double ax, t;
75	uint32_t j0;
76	uint64_t lx;
77	uint16_t hx, ix;
78
79	EXTRACT_LDBL80_WORDS(hx, lx, x);
80	ix = hx & 0x7fff;
81	INSERT_LDBL80_WORDS(ax, ix, lx);
82
83	ENTERI();
84
85	if (ix < 0x3fff) { /* \|x\| < 1 */
86	if (ix < 0x3ffe) { /* \|x\| < 0.5 */
87	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
88	if (x == 0)
89	RETURNI(x);
90	ax = (float)x;
91	x -= ax;
92	t = pilo * x + pihi * x + pilo * ax
93	+ pihi * ax;
94	RETURNI(t);
95	}
96	t = __kernel_tanpil(ax);
97	} else if (ax == 0.5)
98	RETURNI((ax - ax) / (ax - ax));
99	else
100	t = -__kernel_tanpil(1 - ax);
101	RETURNI((hx & 0x8000) ? -t : t);
102	}
103
104	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
105	/* Determine integer part of ax. */
106	j0 = ix - 0x3fff + 1;
107	if (j0 < 32) {
108	lx = (lx >> 32) << 32;
109	lx &= ~(((lx << 32)-1) >> j0);
110	} else {
111	uint64_t m = (uint64_t)-1 >> (j0 + 1);
112	if (lx & m) lx &= ~m;
113	}
114	INSERT_LDBL80_WORDS(x, ix, lx);
115
116	ax -= x;
117	EXTRACT_LDBL80_WORDS(ix, lx, ax);
118
119	if (ix < 0x3ffe) /* \|x\| < 0.5 */
120	t = ax == 0 ? 0 : __kernel_tanpil(ax);
121	else if (ax == 0.5)
122	RETURNI((ax - ax) / (ax - ax));
123	else
124	t = -__kernel_tanpil(1 - ax);
125	RETURNI((hx & 0x8000) ? -t : t);
126	}
127
128	/* x = +-inf or nan. */
129	if (ix >= 0x7fff)
130	RETURNI(vzero / vzero);
131
132	/*
133	* \|x\| >= 0x1p63 is always an integer, so return +-0.
134	* FIXME: should this raise FE_INEXACT or FE_INVALID.
135	*/
136	if (ax + 1 > 1)
137	ax = copysignl(0, x);
138	RETURNI(ax);
139	}

Line 0 Link Here

(-)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

Line 0 Link Here

(-)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

Line 0 Link Here

(-)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

Line 0 Link Here

(-)lib/msun/src/k_cospi.c (+45 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 exploit the kernel for cos(x),
	29	* the argument to __kernel_cospi() must be multiply by pi. As pi is an
	30	* irrational number, this allows the splitting of pi*x into high and low
	31	* parts.
	32	*/
	33
	34	static inline double
	35	__kernel_cospi(double x)
	36	{
	37	double hi, lo;
	38	hi = (float)x;
	39	lo = x - hi;
	40	lo = lo * (pilo + pihi) + hi * pilo;
	41	hi *= pihi;
	42	_2sumF(hi, lo);
	43	return (__kernel_cos(hi, lo));
	44	}
	45

Line 0 Link Here

(-)lib/msun/src/k_sinpi.c (+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 exploit the kernel for sin(x),
	29	* the argument to __kernel_sinpi() must be multiply by pi. As pi is an
	30	* irrational number, this allows the splitting of pi*x into high and low
	31	* parts.
	32	*/
	33
	34	static inline double
	35	__kernel_sinpi(double x)
	36	{
	37	double hi, lo;
	38	hi = (float)x;
	39	lo = x - hi;
	40	lo = lo * (pilo + pihi) + hi * pilo;
	41	hi *= pihi;
	42	_2sumF(hi, lo);
	43	return (__kernel_sin(hi, lo, 1));
	44	}

Line 0 Link Here

(-)lib/msun/src/s_cospi.c (+152 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	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
68	pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
69
70	#include "k_cospi.c"
71	#include "k_sinpi.c"
72
73	double
74	cospi(double x)
75	{
76	volatile static const double vzero = 0;
77	double ax, c;
78	uint32_t hx, ix, j0, lx;
79
80	EXTRACT_WORDS(hx, lx, x);
81	ix = hx & 0x7fffffff;
82	INSERT_WORDS(ax, ix, lx);
83
84	if (ix < 0x3ff00000) { /* \|x\| < 1 */
85	if (ix < 0x3fd00000) { /* \|x\| < 0.25 */
86	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
87	if ((int)ax == 0)
88	return (1);
89	}
90	return (__kernel_cospi(ax));
91	}
92
93	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
94	c = __kernel_sinpi(0.5 - ax);
95	else if (ix < 0x3fe80000){ /* \|x\| < 0.75 */
96	if (ax == 0.5)
97	return (0);
98	c = -__kernel_sinpi(ax - 0.5);
99	} else
100	c = -__kernel_cospi(1 - ax);
101	return (c);
102	}
103
104	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
105	/* Determine integer part of ax. */
106	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
107	if (j0 < 20) {
108	ix &= ~(0x000fffff >> j0);
109	lx = 0;
110	} else {
111	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
112	}
113	INSERT_WORDS(x, ix, lx);
114
115	ax -= x;
116	EXTRACT_WORDS(ix, lx, ax);
117
118
119	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
120	if (ix < 0x3fd00000) /* \|x\| < 0.25 */
121	c = ix == 0 ? 1 : __kernel_cospi(ax);
122	else
123	c = __kernel_sinpi(0.5 - ax);
124	} else {
125	if (ix < 0x3fe80000) { /* \|x\| < 0.75 */
126	if (ax == 0.5)
127	return (0);
128	c = -__kernel_sinpi(ax - 0.5);
129	} else
130	c = -__kernel_cospi(1 - ax);
131	}
132
133	if (j0 > 30) x -= 0x1p30;
134	j0 = (uint32_t)x;
135	return (j0 & 1 ? -c : c);
136	}
137
138	if (ix >= 0x7f800000)
139	return (vzero / vzero);
140
141	/*
142	* \|x\| >= 0x1p52 is always an even integer, so return 1.
143	* FIXME: should this raise FE_INEXACT or FE_INVALID?
144	*/
145	if (ax + 1 > 1)
146	ax = 1;
147	return (ax);
148	}
149
150	#if LDBL_MANT_DIG == 53
151	__weak_reference(cospi, cospil);
152	#endif

Line 0 Link Here

(-)lib/msun/src/s_cospif.c (+111 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) ((float)__kernel_cosdf(M_PI * (x)))
		39	#define __kernel_sinpif(x) ((float)__kernel_sindf(M_PI * (x)))
		40
		41	float
		42	cospif(float x)
		43	{
		44	volatile static const float vzero = 0;
		45	float ax, c;
		46	uint32_t ix, j0;
		47
		48	GET_FLOAT_WORD(ix, x);
		49	ix = ix & 0x7fffffff;
		50	SET_FLOAT_WORD(ax, ix);
		51
		52	if (ix < 0x3f800000) { /* \|x\| < 1 */
		53	if (ix < 0x3e800000) { /* \|x\| < 0.25 */
		54	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
		55	/* Raise inexact iff != 0. */
		56	if ((int)ax == 0)
		57	return (1);
		58	}
		59	return (__kernel_cospif(ax));
		60	}
		61
		62	if (ix < 0x3f000000) /* \|x\| < 0.5 */
		63	c = __kernel_sinpif(0.5f - ax);
		64	else if (ix < 0x3f400000) { /* \|x\| < 0.75 */
65	if (ix == 0x3f000000)
66	return (0);
67	c = -__kernel_sinpif(ax - 0.5f);
68	} else
69	c = -__kernel_cospif(1 - ax);
70	return (c);
71	}
72
73	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
74	/* Determine integer part of ax. */
75	j0 = ((ix >> 23) & 0xff) - 0x7f;
76	ix &= ~(0x007fffff >> j0);
77	SET_FLOAT_WORD(x, ix);
78
79	ax -= x;
80	GET_FLOAT_WORD(ix, ax);
81
82	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
83	if (ix < 0x3e800000) /* \|x\| < 0.25 */
84	c = ix == 0 ? 1 : __kernel_cospif(ax);
85	else
86	c = __kernel_sinpif(0.5f - ax);
87	} else {
88	if (ix < 0x3f400000) { /* \|x\| < 0.75 */
89	if (ix == 0x3f000000)
90	return (0);
91	c = -__kernel_sinpif(ax - 0.5f);
92	} else
93	c = -__kernel_cospif(1 - ax);
94	}
95
96	j0 = (uint32_t)x;
97	return (j0 & 1 ? -c : c);
98	}
99
100	/* x = +-inf or nan. */
101	if (ix >= 0x7f800000)
102	return (vzero / vzero);
103
104	/*
105	* \|x\| >= 0x1p23 is always an even integer, so return 1.
106	* FIXME: should this raise FE_INEXACT or FE_INVALID.
107	*/
108	if (ax + 1 > 1)
109	ax = 1;
110	return (ax);
111	}

Line 0 Link Here

(-)lib/msun/src/s_sinpi.c (+163 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	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
71	pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
72
73	#include "k_cospi.c"
74	#include "k_sinpi.c"
75
76	double
77	sinpi(double x)
78	{
79	volatile static const double vzero = 0;
80	double ax, s;
81	uint32_t hx, ix, j0, lx;
82
83	EXTRACT_WORDS(hx, lx, x);
84	ix = hx & 0x7fffffff;
85	INSERT_WORDS(ax, ix, lx);
86
87	if (ix < 0x3ff00000) { /* \|x\| < 1 */
88	if (ix < 0x3fd00000) { /* \|x\| < 0.25 */
89	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
90	if (x == 0)
91	return (x);
92	INSERT_WORDS(ax, hx, 0);
93	x -= ax;
94	s = pilo * x + pihi * x + pilo * ax
95	+ pihi * ax;
96	return (s);
97	}
98
99	s = __kernel_sinpi(ax);
100	return ((hx & 0x80000000) ? -s : s);
101	}
102
103	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
104	s = __kernel_cospi(0.5 - ax);
105	else if (ix < 0x3fe80000) /* \|x\| < 0.75 */
106	s = __kernel_cospi(ax - 0.5);
107	else
108	s = __kernel_sinpi(1 - ax);
109	return ((hx & 0x80000000) ? -s : s);
110	}
111
112	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
113	/* Determine integer part of ax. */
114	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
115	if (j0 < 20) {
116	ix &= ~(0x000fffff >> j0);
117	lx = 0;
118	} else {
119	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
120	}
121	INSERT_WORDS(x, ix, lx);
122
123	ax -= x;
124	EXTRACT_WORDS(ix, lx, ax);
125
126	if (ix == 0)
127	s = 0;
128	else {
129	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
130	if (ix < 0x3fd00000) /* \|x\| < 0.25 */
131	s = __kernel_sinpi(ax);
132	else
133	s = __kernel_cospi(0.5 - ax);
134	} else {
135	if (ix < 0x3fe80000) /* \|x\| < 0.75 */
136	s = __kernel_cospi(ax - 0.5);
137	else
138	s = __kernel_sinpi(1 - ax);
139	}
140
141	if (j0 > 30) x -= 0x1p30;
142	j0 = (uint32_t)x;
143	if (j0 & 1) s = -s;
144	}
145
146	return ((hx & 0x80000000) ? -s : s);
147	}
148
149	if (ix >= 0x7f800000)
150	return (vzero / vzero);
151
152	/*
153	* \|x\| >= 0x1p52 is always an integer, so return +-0.
154	* FIXME: should this raise FE_INEXACT or FE_INVALID?
155	*/
156	if (ax + 1 > 1)
157	ax = copysign(0, x);
158	return (ax);
159	}
160
161	#if LDBL_MANT_DIG == 53
162	__weak_reference(sinpi, sinpil);
163	#endif

Line 0 Link Here

(-)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) ((float)__kernel_cosdf(M_PI * (x)))
		40	#define __kernel_sinpif(x) ((float)__kernel_sindf(M_PI * (x)))
		41
		42	static const float pihi = 3.14160156e+00f; /* 0x40491000 */
		43	static const float pilo = -8.90890988e-06f; /* 0xb715777a */
		44
		45	float
		46	sinpif(float x)
		47	{
		48	volatile static const float vzero = 0;
		49	float ax, s;
		50	uint32_t hx, ix, j0;
		51
		52	GET_FLOAT_WORD(hx, x);
		53	ix = hx & 0x7fffffff;
		54	SET_FLOAT_WORD(ax, ix);
		55
		56	if (ix < 0x3f800000) { /* \|x\| < 1 */
		57	if (ix < 0x3e800000) { /* \|x\| < 0.25 */
		58	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
		59	if (x == 0)
		60	return (x);
		61	SET_FLOAT_WORD(ax, hx & 0xffff0000);
		62	x -= ax;
		63	s = pilo * x + pihi * x + pilo * ax
		64	+ pihi * ax;
65	return (s);
66	}
67
68	s = __kernel_sinpif(ax);
69	return ((hx & 0x80000000) ? -s : s);
70	}
71
72	if (ix < 0x3f000000) /* \|x\| < 0.5 */
73	s = __kernel_cospif(0.5f - ax);
74	else if (ix < 0x3f400000) /* \|x\| < 0.75 */
75	s = __kernel_cospif(ax - 0.5f);
76	else
77	s = __kernel_sinpif(1 - ax);
78	return ((hx & 0x80000000) ? -s : s);
79	}
80
81	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
82	/* Determine integer part of ax. */
83	j0 = ((ix >> 23) & 0xff) - 0x7f;
84	ix &= ~(0x007fffff >> j0);
85	SET_FLOAT_WORD(x, ix);
86
87	ax -= x;
88	GET_FLOAT_WORD(ix, ax);
89
90	if (ix == 0)
91	s = 0;
92	else {
93	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
94	if (ix < 0x3e800000) /* \|x\| < 0.25 */
95	s = __kernel_sinpif(ax);
96	else
97	s = __kernel_cospif(0.5f - ax);
98	} else {
99	if (ix < 0x3f400000) /* \|x\| < 0.75 */
100	s = __kernel_cospif(ax - 0.5f);
101	else
102	s = __kernel_sinpif(1 - ax);
103	}
104
105	j0 = (uint32_t)x;
106	s = (j0 & 1) ? -s : s;
107	}
108	return ((hx & 0x80000000) ? -s : s);
109	}
110
111	/* x = +-inf or nan. */
112	if (ix >= 0x7f800000)
113	return (vzero / vzero);
114
115	/*
116	* \|x\| >= 0x1p23 is always an integer, so return +-0.
117	* FIXME: should this raise FE_INEXACT or FE_INVALID?
118	*/
119	if (ax + 1 > 1)
120	ax = copysignf(0, x);
121	return (ax);
122	}

Line 0 Link Here

(-)lib/msun/src/s_tanpi.c (+172 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	pihi = 3.1415926814079285e+00, /* 0x400921fb 0x58000000 */
71	pilo = -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 hi, lo, t;
81
82	if (x < 0.25) {
83	hi = (float)x;
84	lo = x - hi;
85	lo = lo * (pilo + pihi) + hi * pilo;
86	hi *= pihi;
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 * (pilo + pihi) + hi * pilo;
94	hi *= pihi;
95	_2sumF(hi, lo);
96	t = - __kernel_tan(hi, lo, -1);
97	} else
98	t = 1;
99
100	return (t);
101	}
102
103	double
104	tanpi(double x)
105	{
106	volatile static const double vzero = 0;
107	double ax, t;
108	uint32_t hx, ix, j0, lx;
109
110	EXTRACT_WORDS(hx, lx, x);
111	ix = hx & 0x7fffffff;
112	INSERT_WORDS(ax, ix, lx);
113
114	if (ix < 0x3ff00000) { /* \|x\| < 1 */
115	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
116	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
117	if (x == 0)
118	return (x);
119	ax = (float)x;
120	x -= ax;
121	t = pilo * x + pihi * x + pilo * ax
122	+ pihi * ax;
123	return (t);
124	}
125	t = __kernel_tanpi(ax);
126	} else if (ax == 0.5)
127	return ((ax - ax) / (ax - ax));
128	else
129	t = - __kernel_tanpi(1 - ax);
130	return ((hx & 0x80000000) ? -t : t);
131	}
132
133	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
134	/* Determine integer part of ax. */
135	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
136	if (j0 < 20) {
137	ix &= ~(0x000fffff >> j0);
138	lx = 0;
139	} else {
140	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
141	}
142	INSERT_WORDS(x,ix,lx);
143
144	ax -= x;
145	EXTRACT_WORDS(ix, lx, ax);
146
147	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
148	t = ax == 0 ? 0 : __kernel_tanpi(ax);
149	else if (ax == 0.5)
150	return ((ax - ax) / (ax - ax));
151	else
152	t = - __kernel_tanpi(1 - ax);
153
154	return ((hx & 0x80000000) ? -t : t);
155	}
156
157	/* x = +-inf or nan. */
158	if (ix >= 0x7f800000)
159	return (vzero / vzero);
160
161	/*
162	* \|x\| >= 0x1p52 is always an integer, so return +-0.
163	* FIXME: should this raise FE_INEXACT or FE_INVALID?
164	*/
165	if (ax + 1 > 1)
166	ax = copysign(0, x);
167	return (ax);
168	}
169
170	#if LDBL_MANT_DIG == 53
171	__weak_reference(tanpi, tanpil);
172	#endif

Line 0 Link Here

(-)lib/msun/src/s_tanpif.c (+115 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	pihi = 3.14160156e+00f, /* 0x40491000 */
		39	pilo = -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	float
		57	tanpif(float x)
		58	{
		59	volatile static const float vzero = 0;
		60	float ax, t;
		61	uint32_t hx, ix, j0;
		62
		63	GET_FLOAT_WORD(hx, x);
		64	ix = hx & 0x7fffffff;
65	SET_FLOAT_WORD(ax, ix);
66
67	if (ix < 0x3f800000) { /* \|x\| < 1 */
68	if (ix < 0x3f000000) { /* \|x\| < 0.5 */
69	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
70	if (ix == 0)
71	return (x);
72	SET_FLOAT_WORD(ax, hx & 0xffff0000);
73	x -= ax;
74	t = pilo * x + pihi * x + pilo * ax
75	+ pihi * ax;
76	return (t);
77	}
78	t = __kernel_tanpif(ax);
79	} else if (ix == 0x3f000000)
80	return ((ax - ax) / (ax - ax));
81	else
82	t = - __kernel_tanpif(1 - ax);
83	return ((hx & 0x80000000) ? -t : t);
84	}
85
86	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
87	/* Determine integer part of ax. */
88	j0 = ((ix >> 23) & 0xff) - 0x7f;
89	ix &= ~(0x007fffff >> j0);
90	SET_FLOAT_WORD(x, ix);
91
92	ax -= x;
93	GET_FLOAT_WORD(ix, ax);
94
95	if (ix < 0x3f000000) /* \|x\| < 0.5 */
96	t = ix == 0 ? 0 : __kernel_tanpif(ax);
97	else if (ix == 0x3f000000)
98	return ((ax - ax) / (ax - ax));
99	else
100	t = - __kernel_tanpif(1 - ax);
101	return ((hx & 0x80000000) ? -t : t);
102	}
103
104	/* x = +-inf or nan. */
105	if (ix >= 0x7f800000)
106	return (vzero / vzero);
107
108	/*
109	* \|x\| >= 0x1p23 is always an integer, so return +-0.
110	* FIXME: should this raise FE_INEXACT or FE_INVALID.
111	*/
112	if (ax + 1 > 1)
113	ax = copysignf(0, x);
114	return (ax);
115	}