Attachment #181622 for bug #218514

Lines 58-64 Link Here

(-)Makefile (-7 / +18 lines)
58	k_tan.c k_tanf.c \	58	k_tan.c k_tanf.c \
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 s_cospi.c s_cospif.c \
62	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \	62	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 \	63	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 \	64	s_finite.c s_finitef.c \
Lines 73-80 Link Here
73	s_nexttowardf.c s_remquo.c s_remquof.c \	73	s_nexttowardf.c s_remquo.c s_remquof.c \
74	s_rint.c s_rintf.c s_round.c s_roundf.c \	74	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 \	75	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 \	76	s_signgam.c s_significand.c s_significandf.c s_sin.c \
77	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \	77	s_sinf.c \
		78	s_sinpi.c s_sinpif.c \
		79	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tanpi.c s_tanpif.c \
		80	s_tgammaf.c s_trunc.c s_truncf.c \
78	w_cabs.c w_cabsf.c w_drem.c w_dremf.c	81	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
79		82
80	# Location of fpmath.h and _fpmath.h	83	# Location of fpmath.h and _fpmath.h
Lines 100-110 Link Here
100	e_lgammal.c e_lgammal_r.c \	103	e_lgammal.c e_lgammal_r.c \
101	e_remainderl.c e_sinhl.c e_sqrtl.c \	104	e_remainderl.c e_sinhl.c e_sqrtl.c \
102	invtrig.c k_cosl.c k_sinl.c k_tanl.c \	105	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 \	106	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cospil.c \
		107	s_cprojl.c \
104	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \	108	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 \	109	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 \	110	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	111	s_scalbnl.c s_sinl.c \
		112	s_sinpil.c s_tanhl.c s_tanl.c s_tanpil.c \
		113	s_truncl.c w_cabsl.c
108	.endif	114	.endif
109		115
110	# C99 complex functions	116	# C99 complex functions
Lines 131-143 Link Here
131		137
132	MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \	138	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 \	139	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 \	140	cimag.3 copysign.3 cos.3 cosh.3 cospi.3 \
		141	csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
135	feclearexcept.3 feenableexcept.3 fegetenv.3 \	142	feclearexcept.3 feenableexcept.3 fegetenv.3 \
136	fegetround.3 fenv.3 floor.3 \	143	fegetround.3 fenv.3 floor.3 \
137	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \	144	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 \	145	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
139	nextafter.3 remainder.3 rint.3 \	146	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 \	147	round.3 scalbn.3 signbit.3 sin.3 sinh.3 sinpi.3 \
		148	sqrt.3 tan.3 tanh.3 tanpi.3 trunc.3 \
141	complex.3	149	complex.3
142		150
143	MLINKS+=acos.3 acosf.3 acos.3 acosl.3	151	MLINKS+=acos.3 acosf.3 acos.3 acosl.3
Lines 166-171 Link Here
166	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3	174	MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
167	MLINKS+=cos.3 cosf.3 cos.3 cosl.3	175	MLINKS+=cos.3 cosf.3 cos.3 cosl.3
168	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3	176	MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
		177	MLINKS+=cospi.3 cospif.3 cospi.3 cospil.3
169	MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3	178	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	179	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 \	180	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	225	MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
217	MLINKS+=sin.3 sinf.3 sin.3 sinl.3	226	MLINKS+=sin.3 sinf.3 sin.3 sinl.3
218	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3	227	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
		228	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 \	229	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
220	sqrt.3 sqrtl.3	230	sqrt.3 sqrtl.3
221	MLINKS+=tan.3 tanf.3 tan.3 tanl.3	231	MLINKS+=tan.3 tanf.3 tan.3 tanl.3
222	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3	232	MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
		233	MLINKS+=tanpi.3 tanpif.3 tanpi.3 tanpil.3
223	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3	234	MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
224		235
225	.include <src.opts.mk>	236	.include <src.opts.mk>

Lines 294-297 Link Here

(-)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.
 *
 * FIXME:  This has not been compiled nor has it been tested for accuracy.
 * FIXME:  This should use bit twiddling.
 */

static const long double
c0   =  1,
c1hi = -4.93480220054467930941724549993807541e+00L,
c1lo = -1.61159244006703024615160135944400084e-34L,
c02  =  4.05871212641676821818501386202937963e+00L,
c03  = -1.33526276885458949587530478285057448e+00L,
c04  =  2.35330630358893204541879352772399571e-01L,
c05  = -2.58068913900140600125982936922503890e-02L,
c06  =  1.92957430940392304790328456072058283e-03L,
c07  = -1.04638104924845707113757250204498056e-04L,
c08  =  4.30306958703294680806457740580669010e-06L,
c09  = -1.38789524622131344829253970500107607e-07L,
c10  =  3.60473079732244660718578413933676910e-09L,
c11  = -7.70070692297072534641720608808474790e-11L,
c12  =  1.37684487309360112303144478033041762e-12L,
c13  = -2.07957304470476734518573046614539573e-14L;

static inline long double
__kernel_cospi(long double x)
{
	union IEEEl2bits u;
	long double c, chi, clo, x2, x2hi, x2lo;

	x2 = x * x;

	c = c09 + (c10 + (c11 + (c12 + c13 * x2) * x2) * x2) * x2;
	c = c05 + (c06 + (c07 + (c08 + c   * x2) * x2) * x2) * x2;
	c =       (c02 + (c03 + (c04 + c   * x2) * x2) * x2) * x2 + c1hi;

	if (x2 < 0x1p-6)		/* |x| < 0x1p-6 */
		 return (c0 + c * x2);

	u.e = x2;
	u.bits.manl = 0;
	x2hi = u.e;
	x2lo = x2 - x2hi;

	u.e = c;
	u.bits.manl = 0;
	chi = u.e;
	clo = c - chi + c1lo;

	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);

	return (c);
}




/*-
 * 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.
 *
 * FIXME:  This has not been compiled nor has it been tested for accuracy.
 * FIXME:  This should use bit twiddling.
 *
 * Compute sinpi(x) = sin(pi*x) in the interval [0,0.25]
 */

static const long double 
s0lo =  3.14159265358979322702026593105983920e+00L,
s0md =  1.14423774522196636802274144823600323e-17L,
s0lo = -5.10459839933399352423791604618557491e-52L,
s01  = -5.16771278004997002924605251118356578e+00L,
s02  =  2.55016403987734544385617758369521912e+00L,
s03  = -5.99264529320792076887739383518724498e-01L,
s04  =  8.21458866111282287988023605285125252e-02L,
s05  = -7.37043094571435077725854643888422319e-03L,
s06  =  4.66302805767612564382677987452304496e-04L,
s07  = -2.19153534478302140522509092972929589e-05L,
s08  =  7.95205400147494466063630770553231865e-07L,
s09  = -2.29484289960219849488354905723109270e-08L,
s10  =  5.39266447760508581521426149354268607e-10L,
s11  = -1.05182947998557665557194979309998767e-11L,
s12  =  1.72036427603492093074818577722479034e-13L;

static inline long double
__kernel_sinpi(long double x)
{
	union IEEEl2bits u;
	long double hi, lo, sm, x2, xhi, xlo;

	/* Split x into high and low parts. */
	u.e = x;
	u.bits.manl = 0;
	xhi = u.e;
	xlo = x - xhi;

	x2 = x * x;
	sm = x2 * (x2 * (x2 * (x2 * s12 + s11) + s10) + s09) + s08;
	sm = x2 * (x2 * (x2 * (x2 * sm  + s07) + s06) + s05) + s04;
	sm = x2 * (x2 * (x2 * sm  + s03) + s02) + s01;
	lo = xlo * (x + xhi) * sm;
	hi = xhi * xhi * sm;
	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;

	return (lo + xhi * s0md + xhi * s0hi);
}




/*-
 * 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"
#include "k_cospi.c"
#include "k_sinpi.c"

static inline long double
__compute_cospil(long double x)
{

	if (x < 0.5) {
		if (x < 0.25)
			return (__kernel_cospil(x));
		return (__kernel_sinpil(0.5 - x));
	} else {
		if (x < 0.75)
			return (-__kernel_sinpil(x - 0.5));
		return (-__kernel_cospil(1 - x));
	}
}

long double
cospil(long double x)
{
	const double huge = 1e+300;
	long double ax, xf;
	uint32_t ix;

	ax = fabsl(x);

	if (ax < 1) {			/* |x| < 1 */
		if (ix < 0x1p-60) {	/* |x| < 0x1p-60 */
			if (huge * x == 0)
				return (1);
			return (1);
		}
		if (ax == 0.5)
			return (0);
		return (__compute_cospil(ax));
	}

	if (ax < 0x1p112) {		/* 1 <= |x| < 0x1p112 */

		xf = floorl(ax);

		ax -= x;
		if (ax == 0.5)
			return (0);

		ax = ax == 0 ? 1 : __compute_cospil(ax);

		if (xf > 0x1p50) xf -= 0x1p50;
		if (xf > 0x1p30) xf -= 0x1p30;
		ix = (uint32_t)xf;

		return (ix & 1 ? -ax : ax);
	}

	if (isinf(x) || isnan(x))
		return (x - x);

	/*
	 * |x| >= 0x1p112 is always an even integer, so return 1.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (1);
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_sinpi.c for implementation details.
 *
 * FIXME:  This has not been compiled nor has it been tested for accuracy.
 * FIXME:  This should use bit twiddling.
 */

#include "math.h"
#include "math_private.h"
#include "k_cospil.c"
#include "k_sinpil.c"

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

static inline long double
__compute_sinpil(long double x)
{

	if (x < 0.5) {
		if (x <= 0.25)
			return (__kernel_sinpil(x));
		return (__kernel_cospil(0.5 - x));
	} else {
		if (x < 0.75)
			return (__kernel_cospil(x - 0.5));
		return (__kernel_sinpil(1 - x));
	}
}

long double
sinpi(long double x)
{
	long double ax, xf;
	uint32_t ix;

	ax = fabsl(x);

	if (ax < 1) {			/* |x| < 1 */
		if (ax < 0x1p-60) {	/* |x| < 0x1p-60 */
			if (x == 0)
				return (x);
			return (pilo * x + pihi * x);
		}
		ax = __compute_sinpil(ax);
		return (copysignl(ax, x));
	}

	if (ax < 0x1p112) {		/* 1 <= |x| < 0x1p112 */

		xf = floorl(ax);

		ax -= xf;
		if (ax == 0) {
			ax = 0;
		} else {
			ax = __compute_sinpil(ax);
			if (xf > 0x1p50) xf -= 0x1p50;
			if (xf > 0x1p30) xf -= 0x1p30;
			ix = (uint32_t)xf;
			if (ix & 1) ax = -ax;
		}
		return (copysignl(ax, x));
	}

	if (isinf(x) || isnan(x))
		return (x - x);

	/*
	 * |x| >= 0x1p112 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysignl(0, x));
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_tanpi.c for implementation details.
 *
 * FIXME: This has not been compiled nor has it been tested for accuracy.
 * FIXME: This should use bit twiddling.
 * 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;
static const long double pilo = 1.14423774522196636802434264184180742e-17L;

static inline long double
__k_tanpi(long double x)
{

	return (sinpil(x) / cospil(x));
}

static inline long double
__compute_tanpi(long double x)
{

	return (x < 0.5 ? __k_tanpil(x) : -__k_tanpil(1 - x));
}

long double
tanpil(long double x)
{
	long double ax, xf;
	uint32_t ix;

	ax = fabsl(ax);

	if (ax < 1) {			/* |x| < 1 */
		if (ax < 0x1p-60) {	/* |x| < 0x1p-60 */
			if (x == 0)
				return (x);
			return (pilo * x + pihi * x);
		}
		if (ax == 0.5)
			return ((x - x) / (x - x));
		ax = __compute_tanpil(ax);
		return (copysignl(ax, x));
	}

	if (ix < 0x1p112) {		/* 1 <= |x| < 0x1p52 */

		xf = floorl(ax);

		ax -= xf;

		if (ax == 0.5)
			return ((x - x) / (x - x));

		ax = ax == 0 ? 0 : __compute_tanpil(ax);
		return (copysignl(ax, x));
	}

	/* x = +-inf or nan. */
	if (isinf(x) || isnan(x))
		return (x - x);

	/*
	 * |x| >= 0x1p53 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysignl(0, x));
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/k_cospi.c for implementation details.
 */

static const union IEEEl2bits
c1hiu = LD80C(0x9de9e64df22ef2d2,   2, -4.93480220054467930927e+00L),
c1lou = LD80C(0xadc4baf0dd791c3e, -63, -1.47187625629472133798e-19L),
c2u   = LD80C(0x81e0f840dad61d9b,   2,  4.05871212641676821836e+00L);
#define	c1hi	(c1hiu.e)
#define	c1lo	(c1lou.e)
#define	c2	(c2u.e)

static const double
c3 = -1.3352627688545895e+00, /* 0xbff55d3c, 0x7e3cbffa */
c4 =  2.3533063035889312e-01, /* 0x3fce1f50, 0x6891bab8 */
c5 = -2.5806891390008115e-02, /* 0xbf9a6d1f, 0x2a2043da */
c6 =  1.9295743091289794e-03, /* 0x3f5f9d38, 0xa362e3d1 */
c7 = -1.0463809745687055e-04, /* 0xbf1b6e24, 0xd372e145 */
c8 =  4.3029513401606985e-06, /* 0x3ed20c42, 0x4210907c */
c9 = -1.3777927680211844e-07; /* 0xbe827e0f, 0x55d52bbb */

static const double c0 = 1;

static inline long double
__kernel_cospil(long double x)
{
  	long double c, chi, clo, x2, x2hi, x2lo;
	uint64_t lx;
	uint16_t ix;

	x2 = x * x;
	c = (c2 + (c3 + (c4 + (c5 + (c6 + (c7 + (c8 + c9 * x2) * x2)
	    * x2) * x2) * x2) * x2) * x2) * x2 + c1hi;

	EXTRACT_LDBL80_WORDS(ix, lx, x2);
	if (ix < 0x3ff9)		/* |x| < 0x1p-6 */
		 return (c0 + c * x2);

	INSERT_LDBL80_WORDS(x2hi, ix, (lx >> 32) << 32);
	x2lo = x2 - x2hi;

	EXTRACT_LDBL80_WORDS(ix, lx, c);
	INSERT_LDBL80_WORDS(chi, ix, (lx >> 40) << 40);
	clo = c - chi + c1lo;

	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);

	return (c);
}




/*-
 * 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 const union IEEEl2bits
s0hiu = LD80C(0xc90fdaa200000000,   1,  3.14159265346825122833e+00L),
s0lou = LD80C(0x85a308d313167b73, -33,  1.21542010130121322891e-10L),
s0mdu = LD80C(0xbe32e3e32bca9d07, -98, -4.68876217131824210625e-30L),
s1u   = LD80C(0xa55de7312df295f5,   2, -5.16771278004997002909e+00L),
s2u   = LD80C(0xa335e33bad570e85,   1,  2.55016403987734544098e+00L);

#define	s0hi	(s0hiu.e)
#define	s0md	(s0mdu.e)
#define	s0lo	(s0lou.e)
#define	s1	(s1u.e)
#define	s2	(s2u.e)

static const double
s3 = -5.9926452932079166e-01,	/* 0xbfe32d2c 0xce62bd82 */
s4 =  8.2145886611090388e-02,	/* 0x3fb50783 0x487edcdb */
s5 = -7.3704309439650085e-03,	/* 0xbf7e3074 0xfdc9c0c9 */
s6 =  4.6630275825043123e-04,	/* 0x3f3e8f43 0x18c2a7a0 */
s7 = -2.1914601021339888e-05,	/* 0xbef6faa7 0xea41adee */
s8 =  7.8877623841310396e-07;	/* 0x3eaa7789 0x4aacabad */

static inline long double
__kernel_sinpil(long double x)
{
	long double hi, lo, sm, x2, xhi, xlo;
	uint64_t lx;
	uint16_t ix;

	EXTRACT_LDBL80_WORDS(ix, lx, x);
	INSERT_LDBL80_WORDS(xhi, ix, (lx >> 32) << 32);
	xlo = x - xhi;

	x2 = x * x;
	sm = x2 * (x2 * (x2 * (x2 * (x2 * (x2 * (x2 * s8 + s7) + s6)
	    + s5) + s4) + s3) + s2) + s1;
	lo = xlo * (x + xhi) * sm;
	hi = xhi * xhi * sm;
	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;

	return (lo + xhi * s0md + xhi * s0hi);
}




/*-
 * 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"
#include "k_cospil.c"
#include "k_sinpil.c"

static inline long double
__compute_cospil(long double x)
{

	if (x < 0.5) {
		if (x < 0.25)
			return (__kernel_cospil(x));
		return (__kernel_sinpil(0.5 - x));
	} else {
		if (x < 0.75)
			return (-__kernel_sinpil(x - 0.5));
		return (-__kernel_cospil(1 - x));
	}
}

long double
cospil(long double x)
{
	const double huge = 1e+300;
	long double ax;
	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 < 0x3fdd) {	/* |x| < 0x1p-34 */
			if (huge * x == 0)
				RETURNI(1);
			RETURNI(1);
		}
		if (ix == 0x3ffe && lx == 0x8000000000000000ull)
			RETURNI(0);
		RETURNI(__compute_cospil(ax));
	}

	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;
		if (ax == 0.5)
			RETURNI(0);

		ax = ax == 0 ? 1 : __compute_cospil(ax);

		if (j0 > 40) x -= 0x1p40;
		if (j0 > 30) x -= 0x1p30;
		j0 = (uint32_t)x;

		return (j0 & 1 ? -ax : ax);
	}

	if (ix >= 0x7fff)
		return (x - x);

	/*
	 * |x| >= 0x1p52 is always an even integer, so return 1.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (1);
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_sinpi.c for implementation details.
 */

#ifdef __i386__
#include <ieeefp.h>
#endif
#include <stdint.h>

#include "fpmath.h"
#include "math.h"
#include "math_private.h"
#include "k_cospil.c"
#include "k_sinpil.c"

static inline long double
__compute_sinpil(long double x)
{

	if (x <= 0.5) {
		if (x <= 0.25)
			return (__kernel_sinpil(x));
		return (__kernel_cospil(0.5 - x));
	} else {
		if (x <= 0.75)
			return (__kernel_cospil(x - 0.5));
		return (__kernel_sinpil(1 - x));
	}
}

static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */

long double
sinpil(long double x)
{
	long double ax;
	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 < 0x3fdc) {	/* |x| < 0x1p-35 */
			if (x == 0)
				RETURNI(x);
			INSERT_LDBL80_WORDS(ax, hx, (lx >> 32) << 32);
			x -= ax;
			RETURNI(pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		ax = __compute_sinpil(ax);
		RETURNI((hx & 0x8000) ? -ax : ax);
	}

	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;
		if (ax == 0) {
			ax = 0;
		} else {
			ax = __compute_sinpil(ax);
			if (j0 > 40) x -= 0x1p40;
			if (j0 > 30) x -= 0x1p30;
			j0 = (uint32_t)x;
			if (j0 & 1) ax = -ax;
		}
		RETURNI((hx & 0x8000) ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7fff)
		RETURNI(x - x);

	/*
	 * |x| >= 0x1p63 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	RETURNI(copysignl(0, x));
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_tanpi.c for implementation details.
 */

#ifdef __i386__
#include <ieeefp.h>
#endif
#include <stdint.h>

#include "fpmath.h"
#include "math.h"
#include "math_private.h"
#include "k_cospil.c"
#include "k_sinpil.c"

static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */

static const union IEEEl2bits
t00hiu = LD80C(0xc90fdaa22168c235,   1,  3.14159265358979323851e+00L),
t00lou = LD80C(0xecdfa9e9520ede79, -65, -5.01599541779735313611e-20L),
t01hiu = LD80C(0xa55de7312df295f5,   3,  1.03354255600999400582e+01L),
t01lou = LD80C(0x8d51cc1ca8412de2, -62,  2.39404580332694307957e-19L),
t02hiu = LD80C(0xa335e33bad570ec7,   5,  4.08026246380375272847e+01L),
t02lou = LD80C(0xe72db1e59703a7ba, -60, -1.56652642070969805771e-18L),
t03u   = LD80C(0xa2fffcda474ac1d3,   7,  1.62999951975255277944e+02L),
t04u   = LD80C(0xa2fa3971d7988829,   9,  6.51909756145994949750e+02L),
t05u   = LD80C(0xa2f99792aec80f33,  11,  2.60759950512193695427e+03L),
t06u   = LD80C(0xa2f985abafb74d99,  13,  1.04303805377441052764e+04L),
t07u   = LD80C(0xa2f983819a4712aa,  15,  4.17215136963294523085e+04L);
#define	t00	(t00hiu.e)
#define	t00lo	(t00lou.e)
#define	t01	(t01hiu.e)
#define	t01lo	(t01lou.e)
#define	t02	(t02hiu.e)
#define	t02lo	(t02lou.e)
#define	t03	(t03u.e)
#define	t04	(t04u.e)
#define	t05	(t05u.e)
#define	t06	(t06u.e)
#define	t07	(t07u.e)

static const double
t08 =  1.6688612339742415e+05, /* 0x41045f30, 0xfcb7c9e9 */
t09 =  6.6753884815733030e+05, /* 0x41245f25, 0xb241ad77 */
t10 =  2.6704954502371075e+06, /* 0x41445fcf, 0xb9a15e9a */
t11 =  1.0666012381827524e+07, /* 0x41645803, 0x8c37ee5b */
t12 =  4.3253268125143178e+07, /* 0x41849ff0, 0xa1004b11 */
t13 =  1.5587209092892912e+08, /* 0x41a294d6, 0xb5db9c99 */
t14 =  1.0162479548145952e+09, /* 0x41ce495b, 0x496844a8 */
t15 = -2.9808831901354914e+09, /* 0xc1e63595, 0x5ec455f2 */
t16 =  8.5833091951381027e+10, /* 0x4233fc0d, 0x0b6f618b */
t17 = -6.8395913998321326e+11, /* 0xc263e7e4, 0x87d1e6d3 */
t18 =  5.1641660627562871e+12, /* 0x4292c981, 0x228a9126 */
t19 = -2.1258706455671141e+13, /* 0xc2b355ad, 0xa58c7724 */
t20 =  5.2041706779955102e+13; /* 0x42c7aa73, 0xb91a998d */

static inline long double
__kernel_tanpil(long double x)
{
	long double t, x2, xhi, xlo;
	uint64_t lx;
	uint16_t ix;

	EXTRACT_LDBL80_WORDS(ix, lx, x);
	INSERT_LDBL80_WORDS(xhi, ix, (lx >> 32) << 32);
	xlo = x - xhi;
	xlo *= (xlo + xhi + xhi);
	xhi *= xhi;
	x2 = x * x;

	if (x < 0.25) {
		t = x2 * (x2 * (x2 * (x2 * t20 + t19) + t18) + t17) + t16;
		t = x2 * (x2 * (x2 * (x2 * t   + t15) + t14) + t13) + t12;
		t = x2 * (x2 * (x2 * (x2 * t   + t11) + t10) + t09) + t08;
		t = x2 * (x2 * (x2 * (x2 * t   + t07) + t06) + t05) + t04;
		t = t02lo + x2 * (x2 * t + t03) + t02;
		t = t01lo + xlo * t + xhi * t + t01;
		t = t00lo + xlo * t + xhi * t + t00;
		return (t * x);
	} else if (x > 0.25) {
		x = 0.5 - x;
		t = __kernel_cospil(x) / __kernel_sinpil(x);
	} else
		t = 1;
	return (t);
}

static inline long double
__compute_tanpil(long double x)
{

	return (x < 0.5 ? __kernel_tanpil(x) : -__kernel_tanpil(1 - x));
}

long double
tanpil(long double x)
{
	long double ax;
	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 < 0x3fdd) {	/* |x| < 0x1p-34 */
			if (x == 0)
				RETURNI(x);
			INSERT_LDBL80_WORDS(ax, ix, (lx >> 32) << 32);
			x -= ax;
			RETURNI(pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		if (ix == 0x3ffe && lx == 0x8000000000000000ull)
			RETURNI((x - x) / (x - x));
		ax = __compute_tanpil(ax);
		RETURNI((hx & 0x8000) ? -ax : ax);
	}

	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;

		if (ax == 0.5)
			RETURNI((x - x) / (x - x));

		ax = ax == 0 ? 0 : __compute_tanpil(ax);
		RETURNI((hx & 0x8000) ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7fff)
		RETURNI(x - x);

	/*
	 * |x| >= 0x1p53 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	RETURNI(copysignl(0, x));
}




.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" All rights reserved.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" are met:
.\" 1. Redistributions of source code must retain the above copyright
.\"    notice, this list of conditions and the following disclaimer.
.\" 2. Redistributions in binary form must reproduce the above copyright
.\"    notice, this list of conditions and the following disclaimer in the
.\"    documentation and/or other materials provided with the distribution.
.\"
.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
.\" ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt COSPI 3
.Os
.Sh NAME
.Nm cospi ,
.Nm cospif ,
.Nm cospil
.Nd half\(encycle cosine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn cospi "double x"
.Ft float
.Fn cospif "float x"
.Ft long double
.Fn cospil "long double x"
.Sh DESCRIPTION
The
.Fn cospi ,
.Fn cospif ,
and
.Fn cospil
functions compute the cosine of
.Fa "\(*p \(mu x" .
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn cospi ,
.Fn cospif ,
and
.Fn cospil
functions returns
.Fn cos "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of 
.Ar x ,
then the returned value is 1 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn cospi \*(Pm0
returns 1.
.It
.Fn cospi \*(Pmn/2
returns 0 for positive integers
.Ar n .
.It
.Fn cospi n
returns 1 for even integers
.Ar n .
.It
.Fn cospi n
returns \-1 for odd integers
.Ar n .
.It
.Fn cospi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn cospi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr sinpi 3 ,
.Xr tan 3 ,
.Xr tanpi 3
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 , 
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.






.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" All rights reserved.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" are met:
.\" 1. Redistributions of source code must retain the above copyright
.\"    notice, this list of conditions and the following disclaimer.
.\" 2. Redistributions in binary form must reproduce the above copyright
.\"    notice, this list of conditions and the following disclaimer in the
.\"    documentation and/or other materials provided with the distribution.
.\"
.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
.\" ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt SINPI 3
.Os
.Sh NAME
.Nm sinpi ,
.Nm sinpif ,
.Nm sinpil
.Nd half\(encycle sine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn sinpi "double x"
.Ft float
.Fn sinpif "float x"
.Ft long double
.Fn sinpil "long double x"
.Sh DESCRIPTION
The
.Fn sinpi ,
.Fn sinpif ,
and
.Fn sinpil
functions compute the sine of
.Fa "\(*p \(mu x" .
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn sinpi ,
.Fn sinpif ,
and
.Fn sinpil
functions returns
.Fn sin "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of 
.Ar x ,
then the returned value is \*(Pm0 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn sinpi \*(Pm0
returns \*(Pm0.
.It
.Fn sinpi \*(Pmn
returns \*(Pm0 for positive integers
.Ar n .
.It
.Fn sinpi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn sinpi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr cospi 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr tan 3 ,
.Xr tanpi 3
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 , 
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.





.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
.\" All rights reserved.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" are met:
.\" 1. Redistributions of source code must retain the above copyright
.\"    notice, this list of conditions and the following disclaimer.
.\" 2. Redistributions in binary form must reproduce the above copyright
.\"    notice, this list of conditions and the following disclaimer in the
.\"    documentation and/or other materials provided with the distribution.
.\"
.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
.\" ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
.\" $FreeBSD$
.\"
.Dd April 1, 2017
.Dt TANPI 3
.Os
.Sh NAME
.Nm tanpi ,
.Nm tanpif ,
.Nm tanpil
.Nd half\(encycle tangent functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft double
.Fn tanpi "double x"
.Ft float
.Fn tanpif "float x"
.Ft long double
.Fn tanpil "long double x"
.Sh DESCRIPTION
The
.Fn tanpi ,
.Fn tanpif ,
and
.Fn tanpil
functions compute the tangent of
.Fa "\(*p \(mu x" 
and measure angles in half-cycles.
.Sh RETURN VALUES
The
.Fn tanpi ,
.Fn tanpif ,
and
.Fn tanpil
functions returns
.Fn tan "\(*p \(mu x" .
If \*(Bax\*(Ba \*(Ge 2^(p - 1)
where p is the floating\(enpoint precision of 
.Ar x ,
then the returned value is \*(Pm0 and it has no significance.
.Sh SPECIAL VALUES
.Bl -tag
.It
.Fn tanpi \*(Pm0
returns \*(Pm0.
.It
.Fn tanpi \*(Pmn
returns \*(Pm0 for positive integers
.Ar n .
.It
.Fn tanpi \*(Pmn/2
returns \*(Na for n > 0 and raises an FE_INVALID exception.
.It
.Fn tanpi \*(Pm\(if
return an \*(Na and raises an FE_INVALID exception.
.It
.Fn tanpi \*(Na
return an \*(Na and raises an FE_INVALID exception.
.El
.Sh SEE ALSO
.Xr cos 3 ,
.Xr cospi 3 ,
.Xr fenv 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr sinpi 3 ,
.Xr tan 3 ,
.Sh AUTHORS
The half\(encycle trignometric functions were written by
.An Steven G. Kargl Aq Mt kargl@FreeBSD.org .
.Sh STANDARDS
These functions conform to
IEEE Std 754\(tm\(en2008 ,
\(dqIEEE Standard for Floating-Point Arithmetic\(dq
and to
ISO/IEC TS 18661-4 , 
\(dqInformation technology \(em Programming languages, their environments,
and system software interfaces \(em Floating\(enpoint extensions for
C\(dq \(em Part 4: Supplementary functions.






/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * Compute cospi(x) = cos(pi*x) via a polynomial approximation.  Coefficients 
 * were determined for the following function:
 *
 * f(x) = ((cos(pi*x) - 1) / (x*x)) = c1 + c2*x**2 + c3*x**4 + ...
 *
 * with an Remes algorithm.  The Remes algorithm used MPFR with a working
 * precision of 30 times that of the target precision (e.g., for 53-bit
 * precision MPFR computed f(x) with 1590-bit precision).  If |x| < 0x1p-6,
 * the polynomial can be summed using a Horner's method.  If |x| exceeds
 * 0x1p-6, then to accumulate the result with sufficient accuracy the c1 
 * coefficient and final multiplications needed to be handled specially.
 *
 * x2 = x * x
 * C(x2) = c1hi + c2*x2 + c3*x2**2 + ...
 *       = c1hi + (c2 + (c3 ...) * x2) * x2	(Horner's method)
 * cospi(x) = c0 + x2 * C(x2)			(Fused multiply add)
 *
 * On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the 
 * interval [0x1p-28,0.25] gave
 *
 *         a = 2.4885199032723904e-01
 *  cospi(a) = 7.0965241314283167e-01,	(0x3fe6b578 0xfa3f3aec)
 * cos(pi*a) = 7.0965241314283156e-01,	(0x3fe6b578 0xfa3f3aeb)
 *
 *          MAX ULP: 0.88903717
 *     Total tested: 134217727
 * 0.8 < ULP <= 0.9: 31605
 * 0.7 < ULP <= 0.8: 382675
 * 0.6 < ULP <= 0.7: 1367966
 */

static const double
c0   =  1.0000000000000000e+00,	/* 0x3ff00000 0x00000000 */
c1hi = -4.9348022005446790e+00,	/* 0xc013bd3c 0xc9be45de */
c1lo = -3.1132368291465357e-16,	/* 0xbcb66ee8 0x86887579 */
c2   =  4.0587121264167649e+00,	/* 0x40103c1f 0x081b5ac0 */
c3   = -1.3352627688538099e+00,	/* 0xbff55d3c 0x7e3cb243 */
c4   =  2.3533063028402096e-01,	/* 0x3fce1f50 0x68689166 */
c5   = -2.5806887965588832e-02,	/* 0xbf9a6d1e 0xef4b8280 */
c6   =  1.9294938875092037e-03,	/* 0x3f5f9ce2 0x49428066 */
c7   = -1.0370089699114971e-04;	/* 0xbf1b2f3f 0xd835d4fe */

static inline double
__kernel_cospi(double x)
{
	double c, chi, clo, x2, x2hi, x2lo;
	uint32_t hx, lx;

	x2 = x * x;
	c = (c2 + (c3 + (c4 + (c5 + (c6 + c7 * x2) * x2) * x2) * x2)
	    * x2) * x2 + c1hi;

	EXTRACT_WORDS(hx, lx, x2);
	if (hx < 0x3f900000)		/* |x| < 0x1p-6 */
		 return (c0 + c * x2);

	INSERT_WORDS(x2hi, hx, 0);
	x2lo = x2 - x2hi;

	EXTRACT_WORDS(hx, lx, c);
	INSERT_WORDS(chi, hx, 0);
	clo = c - chi + c1lo;

	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);

	return (c);
}




/*-
 * 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 const float
c0   =  1.00000000e+00f,	/* 0x3fc00000 */
c1hi = -4.93480206e+00f,	/* 0xc09de9e6 0xb418a08d */
c1lo = -1.42145112e-07f,	/* 0xc09de9e6 0xb418a08d */
c2   =  4.05871058e+00f,	/* 0x4081e0f5 0xb1fc37ea */
c3   = -1.33513784e+00f,	/* 0xbfaae5cc 0xb375e792 */
c4   =  2.32126340e-01f;	/* 0x3e6db287 0x31daa25f */

static inline float
__kernel_cospif(float x)
{
  	float c, chi, clo, x2, x2hi, x2lo;
	uint32_t ix;

	x2 = x * x;
	c = (c2 + (c3 + c4 * x2) * x2) * x2 + c1hi;

	GET_FLOAT_WORD(ix, x2);
	if (ix < 0x3c800000)		/* |x2| < 0x1p-6 */
		 return (c0 + c * x2);
	
	SET_FLOAT_WORD(x2hi, (ix >> 14) << 14);
	x2lo = x2 - x2hi;

	GET_FLOAT_WORD(ix, c);
	SET_FLOAT_WORD(chi, (ix >> 14) << 14);
	clo = c - chi + c1lo;

	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);

	return (c);
}




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

/*
 * Compute sinpi(x) = sin(pi*x) via a polynomial approximation.  Coefficients 
 * were determined for the following function:
 *
 * f(x) = sin(pi*x) / x = s0 + s1*x**2 + s2*x**4 + ...
 *
 * with an Remes algorithm.  The Remes algorithm used MPFR with a working
 * precision of 30 times that of the target precision (e.g., for 53-bit
 * precision MPFR computed f(x) with 1590-bit precision).  To accumulate the
 * result with sufficient accuracy the s0 coefficient and final multiplications
 * needed to be handled specially.
 *
 * x2 = x * x
 * S(x2) = s1 + s2*x2 + s3*x2**2 + ...
 *       = s1 + (s2 + (s3 ...) * x2) * x2		(Horner's method)
 * sinpi(x) = x * (s0 + x2 * S(x2))
 *
 * The final multiplications and additions in the above expression for
 * sinpi(x) are done in an extra precision arithmetic.
 *
 * On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the 
 * interval [0x1p-28,0.25] gave
 *
 *         a =  2.4959746748209000e-01
 *  sinpi(a) =  7.0621201359586339e-01	(0x3fe69949 0xefcdf824)
 * sin(pi*a) =  7.0621201359586350e-01	(0x3fe69949 0xefcdf825)
 *
 *          MAX ULP: 0.82844596
 *     Total tested: 134217727
 * 0.8 < ULP <= 0.9: 13
 * 0.7 < ULP <= 0.8: 9062
 * 0.6 < ULP <= 0.7: 181394
 */

static const double 
s0hi =  3.1415926218032837e+00,	/* 0x400921fb 0x50000000 */
s0md =  3.1786509547050787e-08,	/* 0x3e6110b4 0x611a5f14 */
s0lo =  1.6590549662519743e-24,	/* 0x3b000b9f 0x030b40fc */
s1   = -5.1677127800499703e+00,	/* 0xc014abbc 0xe625be53 */
s2   =  2.5501640398773415e+00,	/* 0x400466bc 0x6775aad9 */
s3   = -5.9926452932029806e-01,	/* 0xbfe32d2c 0xce62ac24 */
s4   =  8.2145886580065136e-02,	/* 0x3fb50783 0x485cc006 */
s5   = -7.3704298849218506e-03,	/* 0xbf7e3074 0xb502de6a */
s6   =  4.6628273194622205e-04,	/* 0x3f3e8eed 0x159b24d7 */
s7   = -2.1717412527405332e-05;	/* 0xbef6c5b9 0x3995dfbe */

static inline double
__kernel_sinpi(double x)
{
	double hi, lo, sm, x2, xhi, xlo;
	uint32_t hx, lx;

	EXTRACT_WORDS(hx, lx, x);
	INSERT_WORDS(xhi, hx, 0);
	xlo = x - xhi;

	x2 = x * x;
	sm = x2 * (x2 * (x2 * (x2 * (x2 * (x2 * s7 + s6) + s5) + s4) + s3)
	    + s2) + s1;
	lo = xlo * (x + xhi) * sm;
	hi = xhi * xhi * sm;
	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;

	return (lo + xhi * s0md + xhi * s0hi);
}




/*-
 * 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 const float
s0hi =  3.14062500e+00f,	/* 0x40490000 */
s0md =  9.67653585e-04f,	/* 0x3a7daa22 */
s0lo = -8.43987033e-12f,	/* 0xad1479cc */
s1   = -5.16771269e+00f,	/* 0xc0a55de7 */
s2   =  2.55016255e+00f,	/* 0x402335dd */
s3   = -5.99202096e-01f,	/* 0xbf19654f */
s4   =  8.10018554e-02f;	/* 0x3da5e44d */

static inline float
__kernel_sinpif(float x)
{
	float hi, lo, sm, x2, xhi, xlo;
	uint32_t ix;

	GET_FLOAT_WORD(ix, x);
	SET_FLOAT_WORD(xhi, (ix >> 14) << 14);
	xlo = x - xhi;

	x2 = x * x;
	sm = x2 * (x2 * (x2 * s4 + s3) + s2) + s1;
	lo = xlo * (x + xhi) * sm;
	hi = xhi * xhi * sm;
	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;

	return (lo + xhi * s0md + xhi * s0hi);
}

Lines 500-505 Link Here

(-)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 = 3 or 4.
 *
 * 2. For |x| < 1, argument reduction is not required and cospi(x) is 
 *    computed by polynomial approximations for sinpi(x) and cospi(x) in the
 *    the interval [0,0.25].  See k_sinpi.c and k_cospi.c for the details for
 *    these polynomial approximations.
 *
 * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
 *    |x| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder.  With
 *    the given domain, a simplified inline j0 = floor(x) is used.
 *
 *    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 polynomial approximation.
 *
 * 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"
#include "k_cospi.c"
#include "k_sinpi.c"

/*
 * To compute cospi(x) for 0 < x < 1, use trignometric identities to map
 * cospi(x) into the [0,0.25] quadrant of sinpi(r) or cospi(r).
 */
static inline double
__compute_cospi(double x)
{

	if (x < 0.5) {
		if (x < 0.25)
			return (__kernel_cospi(x));
		return (__kernel_sinpi(0.5 - x));
	} else {
		if (x < 0.75)
			return (-__kernel_sinpi(x - 0.5));
		return (-__kernel_cospi(1 - x));
	}
}

double
cospi(double x)
{
	const double huge = 1e+300;
	double ax;
	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 < 0x3e200000) {	/* |x| < 0x1p-29 */
			if (huge * x == 0)
				return (1);
			return (1);
		}
		if ((ix|lx) == 0x3fe00000)
			return (0);
		return (__compute_cospi(ax));
	}

	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;
		if (ax == 0.5)
			return (0);

		ax = ax == 0 ? 1 : __compute_cospi(ax);

		if (j0 > 30) x -= 0x1p30;
		j0 = (uint32_t)x;

		return (j0 & 1 ? -ax : ax);
	}

	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p52 is always an even integer, so return 1.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (1);
}

#if LDBL_MANT_DIG == 53
__weak_reference(cospi, cospil);
#endif




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_cospi.c for implementation details.
 */

#include "math.h"
#include "math_private.h"
#include "k_cospif.c"
#include "k_sinpif.c"

static inline float
__compute_cospif(float x)
{

	if (x <= 0.5f) {
		if (x <= 0.25f)
			return (__kernel_cospif(x));
		return (__kernel_sinpif(0.5f - x));
	} else {
		if (x <= 0.75f)
			return (-__kernel_sinpif(x - 0.5f));
		return (-__kernel_cospif(1 - x));
	}
}

float
cospif(float x)
{
	static const float huge = 1e30; 
	float ax;
	uint32_t ix, j0;

	GET_FLOAT_WORD(ix, x);
	ix = ix & 0x7fffffff;
	SET_FLOAT_WORD(ax, ix);

	if (ix < 0x3f800000) {			/* |x| < 1 */
		if (ix < 0x39000000) {		/* |x| < 0x1p-13 */
			if (huge * ax == 0)	/* Raise inexact iff != 0. */
				return (1);
		}
		return (__compute_cospif(ax));
	}

	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;
		if (ax == 0.5f)
			return (0);

		ax = ax == 0 ? 1 : __compute_cospif(ax);

		ix = (uint32_t)x;

		return (ix & 1 ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p23 is always an even integer, so return 1.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (1);
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/**
 * sinpi(x) computes sin(pi*x) without multiplication by pi (almost).  First,
 * note that sinpi(-x) = -sinpi(x), so the algorithm considers only |x| and
 * includes reflection symmetry by considering the sign of x on output.  The
 * method used depends on the magnitude of x.
 *
 * 1. For small |x|, sinpi(x) = pi * x where a sloppy threshold is used.  The
 *    threshold is |x| < 0x1pN with N = -(P/2+M).  P is the precision of the
 *    floating-point type and M = 3 or 4.  To achieve a max ULP less than
 *    0.67, 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 polynomial approximations for sinpi(x) and cospi(x) in the
 *    the interval [0,0.25].  See k_sinpi.c and k_cospi.c for the details for
 *    these polynomial approximations.
 *
 * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
 *    |x| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder.  With
 *    the given domain, a simplified inline j0 = floor(x) is used.
 *
 *    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 polynomial approximation.
 *
 * 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"
#include "k_cospi.c"
#include "k_sinpi.c"

static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */

/*
 * To compute sinpi(x) for 0 < x < 1, use trignometric identities to map
 * sinpi(x) into the [0,0.25] quadrant of sinpi(r) or cospi(r).
 */
static inline double
__compute_sinpi(double x)
{

	if (x < 0.5) {
		if (x <= 0.25)
			return (__kernel_sinpi(x));
		return (__kernel_cospi(0.5 - x));
	} else {
		if (x < 0.75)
			return (__kernel_cospi(x - 0.5));
		return (__kernel_sinpi(1 - x));
	}
}

double
sinpi(double x)
{
	double ax;
	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 < 0x3e300000) {	/* |x| < 0x1p-28 */
			if (x == 0)
				return (x);
			INSERT_WORDS(ax, hx, 0);
			x -= ax;
			return (pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		ax = __compute_sinpi(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	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;
		if (ax == 0) {
			ax = 0;
		} else {
			ax = __compute_sinpi(ax);
			if (j0 > 30) x -= 0x1p30;
			lx = (uint32_t)x;
			if (lx & 1) ax = -ax;
		}
		return ((hx & 0x80000000) ? -ax : ax);
	}

	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p52 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysign(0, x));
}

#if LDBL_MANT_DIG == 53
__weak_reference(sinpi, sinpil);
#endif




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_sinpi.c for implementation details.
 */

#include "math.h"
#include "math_private.h"
#include "k_cospif.c"
#include "k_sinpif.c"

static inline float
__compute_sinpif(float x)
{

	if (x <= 0.5f) {
		if (x <= 0.25f)
			return (__kernel_sinpif(x));
		return (__kernel_cospif(0.5f - x));
	} else {
		if (x <= 0.75f)
			return (__kernel_cospif(x - 0.5f));
		return (__kernel_sinpif(1 - x));
	}
}

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

float
sinpif(float x)
{
	float ax;
	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 < 0x38000000) {	/* |x| < 0x1p-15 */
			if (x == 0)
				return (x);
			SET_FLOAT_WORD(ax, (hx >> 13) << 13);
			x -= ax;
			return (pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		ax = __compute_sinpif(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	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;
		if (ax == 0) {
			ax = 0;
		} else {
			ix = (uint32_t)x;
			ax =  __compute_sinpif(ax);
			if (ix & 1) ax = -ax;
		}
		return ((hx & 0x80000000) ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p23 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysignf(0, x));
}




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/**
 * tanpi(x) computes tan(pi*x) without multiplication by pi (almost).  First,
 * note that tanpi(-x) = -tanpi(x), so the algorithm considers only |x| and
 * includes reflection symmetry by considering the sign of x on output.  The
 * method used depends on the magnitude of x.
 *
 * 1. For small |x|, tanpi(x) = pi * x where a sloppy threshold is used.  The
 *    threshold is |x| < 0x1pN with N = -(P/2+M).  P is the precision of the
 *    floating-point type and M = 3 or 4.  To achieve a max ULP less than
 *    0.59, 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 either a polynomial approximation or the definition of
 *    tan in terms of sin and cos.  For |x| in [0,0.25), a polynomial
 *    approximation is used.  For |x| in (0.25, 0.5], tanpi(x) is computed
 *    by tanpi(x) = cospi(0.5 - x) / sinpi(0.5 - x).
 *
 * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
 *    |x| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder.  With
 *    the given domain, a simplified inline j0 = floor(x) is used.
 *
 *    tanpi(x) = tan(pi*(j0+r))
 *                  tan(pi*j0) + tan(pi*r)
 *             = ----------------------------
 *                1 - tan(pi*j0) * tan(pi*r)
 *             = tanpi(r)
 *
 * 4. For |x| >= 0x1p(P-1), |x| is integral and tanpi(x) = copysign(0,x).
 *
 * 5. Special cases:
 *
 *    tanpif(+-0) = +-0
 *    tanpif(+-n) = +-0, for positive integers n.
 *    tanpif(+-n+1/4) = +-1, for positive integers n.
 *    tanpif(+-n+1/2) = NaN, for positive integers n.
 *    tanpif(+-inf) = nan.  Raises the "invalid" floating-point exception.
 *    tanpif(nan) = nan.  Raises the "invalid" floating-point exception.
 *
 * Compute tanpi(x) = tan(pi*x) via a polynomial approximation.  Coefficients 
 * were determined for the following function:
 *
 * f(x) = tan(pi*x) / x = t0 + t1*x**2 + t2*x**4 + ...
 *
 * with an Remes algorithm.  The Remes algorithm used MPFR with a working
 * precision of 30 times that of the target precision (e.g., for 53-bit
 * precision MPFR computed f(x) with 1590-bit precision).  To accumulate the
 * result with sufficient accuracy the t0, t1, and t2 coefficients and final
 * multiplications and additions needed to be handled specially.
 *
 * x2 = x * x
 * T(x2) = t3 + t4*x2 + t5*x2**2 + ...
 *       = t3 + (t4 + (t5 ...) * x2) * x2		(Horner's method)
 * tanpi(x) = x * (t0 + x2*(t1 + x2*(t2 + x2*T(x2))))
 *
 * The final multiplications and additions in the above expression for
 * tanpi(x) are done in an extra precision arithmetic.
 *
 * On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the 
 * interval [0x1p-28,0.25] gave
 *
 *         a =  2.4250162020325661e-01
 *  tanpi(a) =  9.5396227802027667e-01	(0x3fee86db 0xe636df37)
 * tan(pi*a) =  9.5396227802027678e-01	(0x3fee86db 0xe636df38)
 *
 *          MAX ULP: 1.37582114
 *     Total tested: 134217727
 * 1.0 < ULP       : 309144
 * 0.9 < ULP <= 1.0: 1047636
 * 0.8 < ULP <= 0.9: 2732819
 * 0.7 < ULP <= 0.8: 5491591
 * 0.6 < ULP <= 0.7: 8639239
 *
 * In the interval [0.25,0.48], limited testing finds
 *
 *         a =  4.1875240078199505e-01
 *  tanpi(a) =  3.8323217592887793e+00	(0x400ea898 0x4f7f27eb)
 * tan(pi*a) =  3.8323217592887784e+00	(0x400ea898 0x4f7f27e9)
 * 
 *          MAX ULP: 2.01127812
 *     Total tested: 134217727
 * 1.0 < ULP       : 5702704
 * 0.9 < ULP <= 1.0: 3602502
 * 0.8 < ULP <= 0.9: 5207602
 * 0.7 < ULP <= 0.8: 7173276
 * 0.6 < ULP <= 0.7: 9449882
 */

#include "math.h"
#include "math_private.h"
#include "k_cospi.c"
#include "k_sinpi.c"

static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */

static const double 
t00   =  3.1415926535897931e+00,	/* 0x400921fb 0x54442d18 */
t00lo =  1.2267625636978282e-16,	/* 0x3ca1adf4 0x7b8ee56b */
t01   =  1.0335425560099939e+01,	/* 0x4024abbc 0xe625be52 */
t01lo = -6.0289590484334045e-16,	/* 0xbcc5b8bb 0xb4f3d850 */
t02   =  4.0802624638040442e+01,	/* 0x404466bc 0x6775ac7c */
t02lo = -3.1181465177803440e-15,	/* 0xbcec15f4 0xd355f491 */
t03   =  1.6299995197351416e+02,	/* 0x40645fff 0x9b47fa0c */
t04   =  6.5190975669363422e+02,	/* 0x40845f47 0x2e8473cf */
t05   =  2.6075994007466215e+03,	/* 0x40a45f32 0xe4a797e0 */
t06   =  1.0430393632420044e+04,	/* 0x40c45f32 0x628c115e */
t07   =  4.1720374420664804e+04,	/* 0x40e45f0b 0xfb410bc9 */
t08   =  1.6695720713150888e+05,	/* 0x41046169 0xa8349085 */
t09   =  6.6429163117898861e+05,	/* 0x412445c7 0x4329e474 */
t10   =  2.7802499382503941e+06,	/* 0x4145362c 0xf81896c3 */
t11   =  7.9185101944778729e+06,	/* 0x415e34eb 0x8c725352 */
t12   =  9.3691291294289947e+07,	/* 0x41965676 0x6d2d5a58 */
t13   = -5.0585651903137898e+08,	/* 0xc1be26c2 0x07080874 */
t14   =  6.8772398519624548e+09,	/* 0x41f99ea5 0xa2bf6637 */
t15   = -3.3094464588263161e+10,	/* 0xc21ed255 0xd1310d7a */
t16   =  1.1673794006564838e+11;	/* 0x423b2e1f 0x9a61a5fc */

static inline double
__kernel_tanpi(double x)
{
	double hi, lo, t, x2, xhi, xlo;
	uint32_t hx, lx;

	if (x < 0.25) {
		x2 = x * x;
		t = x2 * (x2 * (x2 * (t15 + x2 * t16) + t14) + t13) + t12;
		t = x2 * (x2 * (x2 * (x2 * t + t11) + t10) + t09) + t08;
		t = x2 * (x2 * (x2 * (x2 * t + t07) + t06) + t05) + t04;

		EXTRACT_WORDS(hx, lx, x);
		INSERT_WORDS(xhi, hx, 0);
		xlo = x - xhi;
		xlo *= (xlo + xhi + xhi);
		xhi *= xhi;

		t = t02lo + x2 * (x2 * t + t03) + t02;
		t = t01lo + xlo * t + xhi * t + t01;
		t = t00lo + xlo * t + xhi * t + t00;
		return (t * x);
	} else if (x > 0.25) {
		x = 0.5 - x;
		t = __kernel_cospi(x) / __kernel_sinpi(x);
	} else
		t = 1;
	return (t);
}

static inline double
__compute_tanpi(double x)
{

	return (x < 0.5 ? __kernel_tanpi(x) : -__kernel_tanpi(1 - x));
}

double
tanpi(double x)
{
	double ax;
	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 < 0x3e200000) {	/* |x| < 0x1p-29 */
			if (x == 0)
				return (x);
			INSERT_WORDS(ax, hx, 0);
			x -= ax;
			return (pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		if (ix == 0x3fe00000 && !lx)
			return ((x - x) / (x - x));
		ax = __compute_tanpi(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	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;

		if (ax == 0.5)
			return ((x - x) / (x - x));

		ax = ax == 0 ? 0 : __compute_tanpi(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p53 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysign(0, x));
}

#if LDBL_MANT_DIG == 53
__weak_reference(tanpi, tanpil);
#endif




/*-
 * Copyright (c) 2017 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice unmodified, this list of conditions, and the following
 *    disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

/*
 * See ../src/s_tanpi.c for implementation details.
 */

#include "math.h"
#include "math_private.h"
#include "k_cospif.c"
#include "k_sinpif.c"

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

static const float
t0   =  3.14159274e+00f,	/* 0x40490fdb */
t0lo = -8.71231194e-08f,	/* 0xb3bb1871 */
t1   =  1.03354244e+01f,	/* 0x41255de6 */
t1lo =  3.83117083e-07f,	/* 0x34cdaf36 */
t2   =  4.08029366e+01f,	/* 0x42233635 */
t2lo = -1.28256147e-07f,	/* 0xb409b6c8 */
t3   =  1.62950455e+02f,	/* 0x4322f351 */
t3lo =  4.45962769e-06f,	/* 0x3695a3e9 */
t4   =  6.55820984e+02f,	/* 0x4423f48b */
t5   =  2.43586987e+03f,	/* 0x45183deb */
t6   =  1.47717539e+04f,	/* 0x4666cf04 */
t7   = -1.96741250e+04f,	/* 0xc699b440 */
t8   =  5.87398688e+05f;	/* 0x490f686b */

static inline float
__kernel_tanpif(float x)
{
	float hi, lo, t, x2, xhi, xlo;
	uint32_t ix;

	if (x < 0.25f) {		/* 0 <= |x| < 0.25 */
		x2 = x * x;
		t = t3lo + x2 * (t4 + x2 * (t5 + x2 * (t6 + x2 * (t7
		    + x2 * t8)))) + t3;

		GET_FLOAT_WORD(ix, x);
		SET_FLOAT_WORD(xhi, (ix >> 14) << 14);
		xlo = x - xhi;
		xlo *= (xlo + xhi + xhi);
		xhi *= xhi;

		t = t2lo + xlo * t + xhi * t + t2;
		t = t1lo + xlo * t + xhi * t + t1;
		t = t0lo + xlo * t + xhi * t + t0;
		t *= x;
	} else if (x > 0.25f) {		/* 0.25 < |x| < 0.5 */
		x = 0.5f - x;
		t = __kernel_cospif(x) / __kernel_sinpif(x);
	} else
		t = 1;
	return (t);
}

static inline float
__compute_tanpif(float x)
{

	return (x < 0.5f ? __kernel_tanpif(x) : -__kernel_tanpif(1 - x));
}

float
tanpif(float x)
{
	float ax;
	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 < 0x38800000) {	/* |x| < 0x1p-14 */
			if (x == 0)
				return (x);
			SET_FLOAT_WORD(ax, (hx >> 13) << 13);
			x -= ax;
			return (pilo * x + pihi * x + pilo * ax + pihi * ax);
		}
		if (ix == 0x3f000000)
			return ((x - x) / (x - x));
		ax = __compute_tanpif(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	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;

		if (ax == 0.5f)
			return ((x - x) / (x - x));

		ax = ax == 0 ? 0 : __compute_tanpif(ax);
		return ((hx & 0x80000000) ? -ax : ax);
	}

	/* x = +-inf or nan. */
	if (ix >= 0x7f800000)
		return (x - x);

	/*
	 * |x| >= 0x1p23 is always an integer, so return +-0.
	 * FIXME: should this raise FE_INEXACT or FE_INVALID.
	 */
	return (copysignf(0, x));
}

Return to bug 218514

Line 0 Link Here

(-)ld128/k_cospil.c (+79 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	* FIXME: This has not been compiled nor has it been tested for accuracy.
		31	* FIXME: This should use bit twiddling.
		32	*/
		33
		34	static const long double
		35	c0 = 1,
		36	c1hi = -4.93480220054467930941724549993807541e+00L,
		37	c1lo = -1.61159244006703024615160135944400084e-34L,
		38	c02 = 4.05871212641676821818501386202937963e+00L,
		39	c03 = -1.33526276885458949587530478285057448e+00L,
		40	c04 = 2.35330630358893204541879352772399571e-01L,
		41	c05 = -2.58068913900140600125982936922503890e-02L,
		42	c06 = 1.92957430940392304790328456072058283e-03L,
		43	c07 = -1.04638104924845707113757250204498056e-04L,
		44	c08 = 4.30306958703294680806457740580669010e-06L,
		45	c09 = -1.38789524622131344829253970500107607e-07L,
		46	c10 = 3.60473079732244660718578413933676910e-09L,
		47	c11 = -7.70070692297072534641720608808474790e-11L,
		48	c12 = 1.37684487309360112303144478033041762e-12L,
		49	c13 = -2.07957304470476734518573046614539573e-14L;
		50
		51	static inline long double
		52	__kernel_cospi(long double x)
		53	{
		54	union IEEEl2bits u;
		55	long double c, chi, clo, x2, x2hi, x2lo;
		56
		57	x2 = x * x;
		58
		59	c = c09 + (c10 + (c11 + (c12 + c13 * x2) * x2) * x2) * x2;
		60	c = c05 + (c06 + (c07 + (c08 + c * x2) * x2) * x2) * x2;
		61	c = (c02 + (c03 + (c04 + c * x2) * x2) * x2) * x2 + c1hi;
		62
		63	if (x2 < 0x1p-6) /* \|x\| < 0x1p-6 */
		64	return (c0 + c * x2);
65
66	u.e = x2;
67	u.bits.manl = 0;
68	x2hi = u.e;
69	x2lo = x2 - x2hi;
70
71	u.e = c;
72	u.bits.manl = 0;
73	chi = u.e;
74	clo = c - chi + c1lo;
75
76	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);
77
78	return (c);
79	}

Line 0 Link Here

(-)ld128/k_sinpil.c (+74 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	* FIXME: This has not been compiled nor has it been tested for accuracy.
		31	* FIXME: This should use bit twiddling.
		32	*
		33	* Compute sinpi(x) = sin(pi*x) in the interval [0,0.25]
		34	*/
		35
		36	static const long double
		37	s0lo = 3.14159265358979322702026593105983920e+00L,
		38	s0md = 1.14423774522196636802274144823600323e-17L,
		39	s0lo = -5.10459839933399352423791604618557491e-52L,
		40	s01 = -5.16771278004997002924605251118356578e+00L,
		41	s02 = 2.55016403987734544385617758369521912e+00L,
		42	s03 = -5.99264529320792076887739383518724498e-01L,
		43	s04 = 8.21458866111282287988023605285125252e-02L,
		44	s05 = -7.37043094571435077725854643888422319e-03L,
		45	s06 = 4.66302805767612564382677987452304496e-04L,
		46	s07 = -2.19153534478302140522509092972929589e-05L,
		47	s08 = 7.95205400147494466063630770553231865e-07L,
		48	s09 = -2.29484289960219849488354905723109270e-08L,
		49	s10 = 5.39266447760508581521426149354268607e-10L,
		50	s11 = -1.05182947998557665557194979309998767e-11L,
		51	s12 = 1.72036427603492093074818577722479034e-13L;
		52
		53	static inline long double
		54	__kernel_sinpi(long double x)
		55	{
		56	union IEEEl2bits u;
		57	long double hi, lo, sm, x2, xhi, xlo;
		58
		59	/* Split x into high and low parts. */
		60	u.e = x;
		61	u.bits.manl = 0;
		62	xhi = u.e;
		63	xlo = x - xhi;
		64
65	x2 = x * x;
66	sm = x2 * (x2 * (x2 * (x2 * s12 + s11) + s10) + s09) + s08;
67	sm = x2 * (x2 * (x2 * (x2 * sm + s07) + s06) + s05) + s04;
68	sm = x2 * (x2 * (x2 * sm + s03) + s02) + s01;
69	lo = xlo * (x + xhi) * sm;
70	hi = xhi * xhi * sm;
71	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;
72
73	return (lo + xhi * s0md + xhi * s0hi);
74	}

Line 0 Link Here

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

Line 0 Link Here

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

Line 0 Link Here

(-)ld128/s_tanpil.c (+98 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 pihi = 3.14159265358979322702026593105983920e+00L;
		40	static const long double pilo = 1.14423774522196636802434264184180742e-17L;
		41
		42	static inline long double
		43	__k_tanpi(long double x)
		44	{
		45
		46	return (sinpil(x) / cospil(x));
		47	}
		48
		49	static inline long double
		50	__compute_tanpi(long double x)
		51	{
		52
		53	return (x < 0.5 ? __k_tanpil(x) : -__k_tanpil(1 - x));
		54	}
		55
		56	long double
		57	tanpil(long double x)
		58	{
		59	long double ax, xf;
		60	uint32_t ix;
		61
		62	ax = fabsl(ax);
		63
		64	if (ax < 1) { /* \|x\| < 1 */
65	if (ax < 0x1p-60) { /* \|x\| < 0x1p-60 */
66	if (x == 0)
67	return (x);
68	return (pilo * x + pihi * x);
69	}
70	if (ax == 0.5)
71	return ((x - x) / (x - x));
72	ax = __compute_tanpil(ax);
73	return (copysignl(ax, x));
74	}
75
76	if (ix < 0x1p112) { /* 1 <= \|x\| < 0x1p52 */
77
78	xf = floorl(ax);
79
80	ax -= xf;
81
82	if (ax == 0.5)
83	return ((x - x) / (x - x));
84
85	ax = ax == 0 ? 0 : __compute_tanpil(ax);
86	return (copysignl(ax, x));
87	}
88
89	/* x = +-inf or nan. */
90	if (isinf(x) \|\| isnan(x))
91	return (x - x);
92
93	/*
94	* \|x\| >= 0x1p53 is always an integer, so return +-0.
95	* FIXME: should this raise FE_INEXACT or FE_INVALID.
96	*/
97	return (copysignl(0, x));
98	}

Line 0 Link Here

(-)ld80/k_cospil.c (+75 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 const union IEEEl2bits
		32	c1hiu = LD80C(0x9de9e64df22ef2d2, 2, -4.93480220054467930927e+00L),
		33	c1lou = LD80C(0xadc4baf0dd791c3e, -63, -1.47187625629472133798e-19L),
		34	c2u = LD80C(0x81e0f840dad61d9b, 2, 4.05871212641676821836e+00L);
		35	#define c1hi (c1hiu.e)
		36	#define c1lo (c1lou.e)
		37	#define c2 (c2u.e)
		38
		39	static const double
		40	c3 = -1.3352627688545895e+00, /* 0xbff55d3c, 0x7e3cbffa */
		41	c4 = 2.3533063035889312e-01, /* 0x3fce1f50, 0x6891bab8 */
		42	c5 = -2.5806891390008115e-02, /* 0xbf9a6d1f, 0x2a2043da */
		43	c6 = 1.9295743091289794e-03, /* 0x3f5f9d38, 0xa362e3d1 */
		44	c7 = -1.0463809745687055e-04, /* 0xbf1b6e24, 0xd372e145 */
		45	c8 = 4.3029513401606985e-06, /* 0x3ed20c42, 0x4210907c */
		46	c9 = -1.3777927680211844e-07; /* 0xbe827e0f, 0x55d52bbb */
		47
		48	static const double c0 = 1;
		49
		50	static inline long double
		51	__kernel_cospil(long double x)
		52	{
		53	long double c, chi, clo, x2, x2hi, x2lo;
		54	uint64_t lx;
		55	uint16_t ix;
		56
		57	x2 = x * x;
		58	c = (c2 + (c3 + (c4 + (c5 + (c6 + (c7 + (c8 + c9 * x2) * x2)
		59	* x2) * x2) * x2) * x2) * x2) * x2 + c1hi;
		60
		61	EXTRACT_LDBL80_WORDS(ix, lx, x2);
		62	if (ix < 0x3ff9) /* \|x\| < 0x1p-6 */
		63	return (c0 + c * x2);
		64
65	INSERT_LDBL80_WORDS(x2hi, ix, (lx >> 32) << 32);
66	x2lo = x2 - x2hi;
67
68	EXTRACT_LDBL80_WORDS(ix, lx, c);
69	INSERT_LDBL80_WORDS(chi, ix, (lx >> 40) << 40);
70	clo = c - chi + c1lo;
71
72	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);
73
74	return (c);
75	}

Line 0 Link Here

(-)ld80/k_sinpil.c (+71 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 const union IEEEl2bits
		32	s0hiu = LD80C(0xc90fdaa200000000, 1, 3.14159265346825122833e+00L),
		33	s0lou = LD80C(0x85a308d313167b73, -33, 1.21542010130121322891e-10L),
		34	s0mdu = LD80C(0xbe32e3e32bca9d07, -98, -4.68876217131824210625e-30L),
		35	s1u = LD80C(0xa55de7312df295f5, 2, -5.16771278004997002909e+00L),
		36	s2u = LD80C(0xa335e33bad570e85, 1, 2.55016403987734544098e+00L);
		37
		38	#define s0hi (s0hiu.e)
		39	#define s0md (s0mdu.e)
		40	#define s0lo (s0lou.e)
		41	#define s1 (s1u.e)
		42	#define s2 (s2u.e)
		43
		44	static const double
		45	s3 = -5.9926452932079166e-01, /* 0xbfe32d2c 0xce62bd82 */
		46	s4 = 8.2145886611090388e-02, /* 0x3fb50783 0x487edcdb */
		47	s5 = -7.3704309439650085e-03, /* 0xbf7e3074 0xfdc9c0c9 */
		48	s6 = 4.6630275825043123e-04, /* 0x3f3e8f43 0x18c2a7a0 */
		49	s7 = -2.1914601021339888e-05, /* 0xbef6faa7 0xea41adee */
		50	s8 = 7.8877623841310396e-07; /* 0x3eaa7789 0x4aacabad */
		51
		52	static inline long double
		53	__kernel_sinpil(long double x)
		54	{
		55	long double hi, lo, sm, x2, xhi, xlo;
		56	uint64_t lx;
		57	uint16_t ix;
		58
		59	EXTRACT_LDBL80_WORDS(ix, lx, x);
		60	INSERT_LDBL80_WORDS(xhi, ix, (lx >> 32) << 32);
		61	xlo = x - xhi;
		62
		63	x2 = x * x;
		64	sm = x2 * (x2 * (x2 * (x2 * (x2 * (x2 * (x2 * s8 + s7) + s6)
65	+ s5) + s4) + s3) + s2) + s1;
66	lo = xlo * (x + xhi) * sm;
67	hi = xhi * xhi * sm;
68	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;
69
70	return (lo + xhi * s0md + xhi * s0hi);
71	}

Line 0 Link Here

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

Line 0 Link Here

(-)ld80/s_sinpil.c (+120 lines)
		1	/*-
		2	* Copyright (c) 2017 Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	/*
		28	* See ../src/s_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	#include "k_cospil.c"
		40	#include "k_sinpil.c"
		41
		42	static inline long double
		43	__compute_sinpil(long double x)
		44	{
		45
		46	if (x <= 0.5) {
		47	if (x <= 0.25)
		48	return (__kernel_sinpil(x));
		49	return (__kernel_cospil(0.5 - x));
		50	} else {
		51	if (x <= 0.75)
		52	return (__kernel_cospil(x - 0.5));
		53	return (__kernel_sinpil(1 - x));
		54	}
		55	}
		56
		57	static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
		58	static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
		59
		60	long double
		61	sinpil(long double x)
		62	{
		63	long double ax;
		64	uint32_t j0;
65	uint64_t lx;
66	uint16_t hx, ix;
67
68	EXTRACT_LDBL80_WORDS(hx, lx, x);
69	ix = hx & 0x7fff;
70	INSERT_LDBL80_WORDS(ax, ix, lx);
71
72	ENTERI();
73
74	if (ix < 0x3fff) { /* \|x\| < 1 */
75	if (ix < 0x3fdc) { /* \|x\| < 0x1p-35 */
76	if (x == 0)
77	RETURNI(x);
78	INSERT_LDBL80_WORDS(ax, hx, (lx >> 32) << 32);
79	x -= ax;
80	RETURNI(pilo * x + pihi * x + pilo * ax + pihi * ax);
81	}
82	ax = __compute_sinpil(ax);
83	RETURNI((hx & 0x8000) ? -ax : ax);
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	if (ax == 0) {
100	ax = 0;
101	} else {
102	ax = __compute_sinpil(ax);
103	if (j0 > 40) x -= 0x1p40;
104	if (j0 > 30) x -= 0x1p30;
105	j0 = (uint32_t)x;
106	if (j0 & 1) ax = -ax;
107	}
108	RETURNI((hx & 0x8000) ? -ax : ax);
109	}
110
111	/* x = +-inf or nan. */
112	if (ix >= 0x7fff)
113	RETURNI(x - x);
114
115	/*
116	* \|x\| >= 0x1p63 is always an integer, so return +-0.
117	* FIXME: should this raise FE_INEXACT or FE_INVALID.
118	*/
119	RETURNI(copysignl(0, x));
120	}

Line 0 Link Here

(-)ld80/s_tanpil.c (+180 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	#include "k_cospil.c"
		40	#include "k_sinpil.c"
		41
		42	static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
		43	static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
		44
		45	static const union IEEEl2bits
		46	t00hiu = LD80C(0xc90fdaa22168c235, 1, 3.14159265358979323851e+00L),
		47	t00lou = LD80C(0xecdfa9e9520ede79, -65, -5.01599541779735313611e-20L),
		48	t01hiu = LD80C(0xa55de7312df295f5, 3, 1.03354255600999400582e+01L),
		49	t01lou = LD80C(0x8d51cc1ca8412de2, -62, 2.39404580332694307957e-19L),
		50	t02hiu = LD80C(0xa335e33bad570ec7, 5, 4.08026246380375272847e+01L),
		51	t02lou = LD80C(0xe72db1e59703a7ba, -60, -1.56652642070969805771e-18L),
		52	t03u = LD80C(0xa2fffcda474ac1d3, 7, 1.62999951975255277944e+02L),
		53	t04u = LD80C(0xa2fa3971d7988829, 9, 6.51909756145994949750e+02L),
		54	t05u = LD80C(0xa2f99792aec80f33, 11, 2.60759950512193695427e+03L),
		55	t06u = LD80C(0xa2f985abafb74d99, 13, 1.04303805377441052764e+04L),
		56	t07u = LD80C(0xa2f983819a4712aa, 15, 4.17215136963294523085e+04L);
		57	#define t00 (t00hiu.e)
		58	#define t00lo (t00lou.e)
		59	#define t01 (t01hiu.e)
		60	#define t01lo (t01lou.e)
		61	#define t02 (t02hiu.e)
		62	#define t02lo (t02lou.e)
		63	#define t03 (t03u.e)
		64	#define t04 (t04u.e)
65	#define t05 (t05u.e)
66	#define t06 (t06u.e)
67	#define t07 (t07u.e)
68
69	static const double
70	t08 = 1.6688612339742415e+05, /* 0x41045f30, 0xfcb7c9e9 */
71	t09 = 6.6753884815733030e+05, /* 0x41245f25, 0xb241ad77 */
72	t10 = 2.6704954502371075e+06, /* 0x41445fcf, 0xb9a15e9a */
73	t11 = 1.0666012381827524e+07, /* 0x41645803, 0x8c37ee5b */
74	t12 = 4.3253268125143178e+07, /* 0x41849ff0, 0xa1004b11 */
75	t13 = 1.5587209092892912e+08, /* 0x41a294d6, 0xb5db9c99 */
76	t14 = 1.0162479548145952e+09, /* 0x41ce495b, 0x496844a8 */
77	t15 = -2.9808831901354914e+09, /* 0xc1e63595, 0x5ec455f2 */
78	t16 = 8.5833091951381027e+10, /* 0x4233fc0d, 0x0b6f618b */
79	t17 = -6.8395913998321326e+11, /* 0xc263e7e4, 0x87d1e6d3 */
80	t18 = 5.1641660627562871e+12, /* 0x4292c981, 0x228a9126 */
81	t19 = -2.1258706455671141e+13, /* 0xc2b355ad, 0xa58c7724 */
82	t20 = 5.2041706779955102e+13; /* 0x42c7aa73, 0xb91a998d */
83
84	static inline long double
85	__kernel_tanpil(long double x)
86	{
87	long double t, x2, xhi, xlo;
88	uint64_t lx;
89	uint16_t ix;
90
91	EXTRACT_LDBL80_WORDS(ix, lx, x);
92	INSERT_LDBL80_WORDS(xhi, ix, (lx >> 32) << 32);
93	xlo = x - xhi;
94	xlo *= (xlo + xhi + xhi);
95	xhi *= xhi;
96	x2 = x * x;
97
98	if (x < 0.25) {
99	t = x2 * (x2 * (x2 * (x2 * t20 + t19) + t18) + t17) + t16;
100	t = x2 * (x2 * (x2 * (x2 * t + t15) + t14) + t13) + t12;
101	t = x2 * (x2 * (x2 * (x2 * t + t11) + t10) + t09) + t08;
102	t = x2 * (x2 * (x2 * (x2 * t + t07) + t06) + t05) + t04;
103	t = t02lo + x2 * (x2 * t + t03) + t02;
104	t = t01lo + xlo * t + xhi * t + t01;
105	t = t00lo + xlo * t + xhi * t + t00;
106	return (t * x);
107	} else if (x > 0.25) {
108	x = 0.5 - x;
109	t = __kernel_cospil(x) / __kernel_sinpil(x);
110	} else
111	t = 1;
112	return (t);
113	}
114
115	static inline long double
116	__compute_tanpil(long double x)
117	{
118
119	return (x < 0.5 ? __kernel_tanpil(x) : -__kernel_tanpil(1 - x));
120	}
121
122	long double
123	tanpil(long double x)
124	{
125	long double ax;
126	uint32_t j0;
127	uint64_t lx;
128	uint16_t hx, ix;
129
130	EXTRACT_LDBL80_WORDS(hx, lx, x);
131	ix = hx & 0x7fff;
132	INSERT_LDBL80_WORDS(ax, ix, lx);
133
134	ENTERI();
135
136	if (ix < 0x3fff) { /* \|x\| < 1 */
137	if (ix < 0x3fdd) { /* \|x\| < 0x1p-34 */
138	if (x == 0)
139	RETURNI(x);
140	INSERT_LDBL80_WORDS(ax, ix, (lx >> 32) << 32);
141	x -= ax;
142	RETURNI(pilo * x + pihi * x + pilo * ax + pihi * ax);
143	}
144	if (ix == 0x3ffe && lx == 0x8000000000000000ull)
145	RETURNI((x - x) / (x - x));
146	ax = __compute_tanpil(ax);
147	RETURNI((hx & 0x8000) ? -ax : ax);
148	}
149
150	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
151	/* Determine integer part of ax. */
152	j0 = ix - 0x3fff + 1;
153	if (j0 < 32) {
154	lx = (lx >> 32) << 32;
155	lx &= ~(((lx << 32)-1) >> j0);
156	} else {
157	uint64_t m = (uint64_t)-1 >> (j0 + 1);
158	if (lx & m) lx &= ~m;
159	}
160	INSERT_LDBL80_WORDS(x, ix, lx);
161
162	ax -= x;
163
164	if (ax == 0.5)
165	RETURNI((x - x) / (x - x));
166
167	ax = ax == 0 ? 0 : __compute_tanpil(ax);
168	RETURNI((hx & 0x8000) ? -ax : ax);
169	}
170
171	/* x = +-inf or nan. */
172	if (ix >= 0x7fff)
173	RETURNI(x - x);
174
175	/*
176	* \|x\| >= 0x1p53 is always an integer, so return +-0.
177	* FIXME: should this raise FE_INEXACT or FE_INVALID.
178	*/
179	RETURNI(copysignl(0, x));
180	}

Line 0 Link Here

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

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

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

(-)src/k_cospi.c (+94 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	* Compute cospi(x) = cos(pi*x) via a polynomial approximation. Coefficients
		29	* were determined for the following function:
		30	*
		31	* f(x) = ((cos(pix) - 1) / (xx)) = c1 + c2x2 + c3x**4 + ...
		32	*
		33	* with an Remes algorithm. The Remes algorithm used MPFR with a working
		34	* precision of 30 times that of the target precision (e.g., for 53-bit
		35	* precision MPFR computed f(x) with 1590-bit precision). If \|x\| < 0x1p-6,
		36	* the polynomial can be summed using a Horner's method. If \|x\| exceeds
		37	* 0x1p-6, then to accumulate the result with sufficient accuracy the c1
		38	* coefficient and final multiplications needed to be handled specially.
		39	*
		40	* x2 = x * x
		41	* C(x2) = c1hi + c2x2 + c3x2**2 + ...
		42	* = c1hi + (c2 + (c3 ...) * x2) * x2 (Horner's method)
		43	* cospi(x) = c0 + x2 * C(x2) (Fused multiply add)
		44	*
		45	* On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the
		46	* interval [0x1p-28,0.25] gave
		47	*
		48	* a = 2.4885199032723904e-01
		49	* cospi(a) = 7.0965241314283167e-01, (0x3fe6b578 0xfa3f3aec)
		50	* cos(pi*a) = 7.0965241314283156e-01, (0x3fe6b578 0xfa3f3aeb)
		51	*
		52	* MAX ULP: 0.88903717
		53	* Total tested: 134217727
		54	* 0.8 < ULP <= 0.9: 31605
		55	* 0.7 < ULP <= 0.8: 382675
		56	* 0.6 < ULP <= 0.7: 1367966
		57	*/
		58
		59	static const double
		60	c0 = 1.0000000000000000e+00, /* 0x3ff00000 0x00000000 */
		61	c1hi = -4.9348022005446790e+00, /* 0xc013bd3c 0xc9be45de */
		62	c1lo = -3.1132368291465357e-16, /* 0xbcb66ee8 0x86887579 */
		63	c2 = 4.0587121264167649e+00, /* 0x40103c1f 0x081b5ac0 */
		64	c3 = -1.3352627688538099e+00, /* 0xbff55d3c 0x7e3cb243 */
65	c4 = 2.3533063028402096e-01, /* 0x3fce1f50 0x68689166 */
66	c5 = -2.5806887965588832e-02, /* 0xbf9a6d1e 0xef4b8280 */
67	c6 = 1.9294938875092037e-03, /* 0x3f5f9ce2 0x49428066 */
68	c7 = -1.0370089699114971e-04; /* 0xbf1b2f3f 0xd835d4fe */
69
70	static inline double
71	__kernel_cospi(double x)
72	{
73	double c, chi, clo, x2, x2hi, x2lo;
74	uint32_t hx, lx;
75
76	x2 = x * x;
77	c = (c2 + (c3 + (c4 + (c5 + (c6 + c7 * x2) * x2) * x2) * x2)
78	* x2) * x2 + c1hi;
79
80	EXTRACT_WORDS(hx, lx, x2);
81	if (hx < 0x3f900000) /* \|x\| < 0x1p-6 */
82	return (c0 + c * x2);
83
84	INSERT_WORDS(x2hi, hx, 0);
85	x2lo = x2 - x2hi;
86
87	EXTRACT_WORDS(hx, lx, c);
88	INSERT_WORDS(chi, hx, 0);
89	clo = c - chi + c1lo;
90
91	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);
92
93	return (c);
94	}

Line 0 Link Here

(-)src/k_cospif.c (+62 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 const float
	32	c0 = 1.00000000e+00f, /* 0x3fc00000 */
	33	c1hi = -4.93480206e+00f, /* 0xc09de9e6 0xb418a08d */
	34	c1lo = -1.42145112e-07f, /* 0xc09de9e6 0xb418a08d */
	35	c2 = 4.05871058e+00f, /* 0x4081e0f5 0xb1fc37ea */
	36	c3 = -1.33513784e+00f, /* 0xbfaae5cc 0xb375e792 */
	37	c4 = 2.32126340e-01f; /* 0x3e6db287 0x31daa25f */
	38
	39	static inline float
	40	__kernel_cospif(float x)
	41	{
	42	float c, chi, clo, x2, x2hi, x2lo;
	43	uint32_t ix;
	44
	45	x2 = x * x;
	46	c = (c2 + (c3 + c4 * x2) * x2) * x2 + c1hi;
	47
	48	GET_FLOAT_WORD(ix, x2);
	49	if (ix < 0x3c800000) /* \|x2\| < 0x1p-6 */
	50	return (c0 + c * x2);
	51
	52	SET_FLOAT_WORD(x2hi, (ix >> 14) << 14);
	53	x2lo = x2 - x2hi;
	54
	55	GET_FLOAT_WORD(ix, c);
	56	SET_FLOAT_WORD(chi, (ix >> 14) << 14);
	57	clo = c - chi + c1lo;
	58
	59	c = clo * x2lo + chi * x2lo + clo * x2hi + (chi * x2hi + c0);
	60
	61	return (c);
	62	}

Line 0 Link Here

(-)src/k_sinpi.c (+91 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	* Compute sinpi(x) = sin(pi*x) via a polynomial approximation. Coefficients
		29	* were determined for the following function:
		30	*
		31	* f(x) = sin(pix) / x = s0 + s1x*2 + s2x**4 + ...
		32	*
		33	* with an Remes algorithm. The Remes algorithm used MPFR with a working
		34	* precision of 30 times that of the target precision (e.g., for 53-bit
		35	* precision MPFR computed f(x) with 1590-bit precision). To accumulate the
		36	* result with sufficient accuracy the s0 coefficient and final multiplications
		37	* needed to be handled specially.
		38	*
		39	* x2 = x * x
		40	* S(x2) = s1 + s2x2 + s3x2**2 + ...
		41	* = s1 + (s2 + (s3 ...) * x2) * x2 (Horner's method)
		42	* sinpi(x) = x * (s0 + x2 * S(x2))
		43	*
		44	* The final multiplications and additions in the above expression for
		45	* sinpi(x) are done in an extra precision arithmetic.
		46	*
		47	* On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the
		48	* interval [0x1p-28,0.25] gave
		49	*
		50	* a = 2.4959746748209000e-01
		51	* sinpi(a) = 7.0621201359586339e-01 (0x3fe69949 0xefcdf824)
		52	* sin(pi*a) = 7.0621201359586350e-01 (0x3fe69949 0xefcdf825)
		53	*
		54	* MAX ULP: 0.82844596
		55	* Total tested: 134217727
		56	* 0.8 < ULP <= 0.9: 13
		57	* 0.7 < ULP <= 0.8: 9062
		58	* 0.6 < ULP <= 0.7: 181394
		59	*/
		60
		61	static const double
		62	s0hi = 3.1415926218032837e+00, /* 0x400921fb 0x50000000 */
		63	s0md = 3.1786509547050787e-08, /* 0x3e6110b4 0x611a5f14 */
		64	s0lo = 1.6590549662519743e-24, /* 0x3b000b9f 0x030b40fc */
65	s1 = -5.1677127800499703e+00, /* 0xc014abbc 0xe625be53 */
66	s2 = 2.5501640398773415e+00, /* 0x400466bc 0x6775aad9 */
67	s3 = -5.9926452932029806e-01, /* 0xbfe32d2c 0xce62ac24 */
68	s4 = 8.2145886580065136e-02, /* 0x3fb50783 0x485cc006 */
69	s5 = -7.3704298849218506e-03, /* 0xbf7e3074 0xb502de6a */
70	s6 = 4.6628273194622205e-04, /* 0x3f3e8eed 0x159b24d7 */
71	s7 = -2.1717412527405332e-05; /* 0xbef6c5b9 0x3995dfbe */
72
73	static inline double
74	__kernel_sinpi(double x)
75	{
76	double hi, lo, sm, x2, xhi, xlo;
77	uint32_t hx, lx;
78
79	EXTRACT_WORDS(hx, lx, x);
80	INSERT_WORDS(xhi, hx, 0);
81	xlo = x - xhi;
82
83	x2 = x * x;
84	sm = x2 * (x2 * (x2 * (x2 * (x2 * (x2 * s7 + s6) + s5) + s4) + s3)
85	+ s2) + s1;
86	lo = xlo * (x + xhi) * sm;
87	hi = xhi * xhi * sm;
88	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;
89
90	return (lo + xhi * s0md + xhi * s0hi);
91	}

Line 0 Link Here

(-)src/k_sinpif.c (+57 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 const float
	32	s0hi = 3.14062500e+00f, /* 0x40490000 */
	33	s0md = 9.67653585e-04f, /* 0x3a7daa22 */
	34	s0lo = -8.43987033e-12f, /* 0xad1479cc */
	35	s1 = -5.16771269e+00f, /* 0xc0a55de7 */
	36	s2 = 2.55016255e+00f, /* 0x402335dd */
	37	s3 = -5.99202096e-01f, /* 0xbf19654f */
	38	s4 = 8.10018554e-02f; /* 0x3da5e44d */
	39
	40	static inline float
	41	__kernel_sinpif(float x)
	42	{
	43	float hi, lo, sm, x2, xhi, xlo;
	44	uint32_t ix;
	45
	46	GET_FLOAT_WORD(ix, x);
	47	SET_FLOAT_WORD(xhi, (ix >> 14) << 14);
	48	xlo = x - xhi;
	49
	50	x2 = x * x;
	51	sm = x2 * (x2 * (x2 * s4 + s3) + s2) + s1;
	52	lo = xlo * (x + xhi) * sm;
	53	hi = xhi * xhi * sm;
	54	lo = x * (lo + s0lo) + xlo * (s0hi + s0md) + x * hi;
	55
	56	return (lo + xhi * s0md + xhi * s0hi);
	57	}

Line 0 Link Here

(-)src/s_cospi.c (+144 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 = 3 or 4.
		35	*
		36	* 2. For \|x\| < 1, argument reduction is not required and cospi(x) is
		37	* computed by polynomial approximations for sinpi(x) and cospi(x) in the
		38	* the interval [0,0.25]. See k_sinpi.c and k_cospi.c for the details for
		39	* these polynomial approximations.
		40	*
		41	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
		42	* \|x\| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder. With
		43	* the given domain, a simplified inline j0 = floor(x) is used.
		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 polynomial approximation.
		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	#include "k_cospi.c"
66	#include "k_sinpi.c"
67
68	/*
69	* To compute cospi(x) for 0 < x < 1, use trignometric identities to map
70	* cospi(x) into the [0,0.25] quadrant of sinpi(r) or cospi(r).
71	*/
72	static inline double
73	__compute_cospi(double x)
74	{
75
76	if (x < 0.5) {
77	if (x < 0.25)
78	return (__kernel_cospi(x));
79	return (__kernel_sinpi(0.5 - x));
80	} else {
81	if (x < 0.75)
82	return (-__kernel_sinpi(x - 0.5));
83	return (-__kernel_cospi(1 - x));
84	}
85	}
86
87	double
88	cospi(double x)
89	{
90	const double huge = 1e+300;
91	double ax;
92	uint32_t hx, ix, j0, lx;
93
94	EXTRACT_WORDS(hx, lx, x);
95	ix = hx & 0x7fffffff;
96	INSERT_WORDS(ax, ix, lx);
97
98	if (ix < 0x3ff00000) { /* \|x\| < 1 */
99	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
100	if (huge * x == 0)
101	return (1);
102	return (1);
103	}
104	if ((ix\|lx) == 0x3fe00000)
105	return (0);
106	return (__compute_cospi(ax));
107	}
108
109	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
110	/* Determine integer part of ax. */
111	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
112	if (j0 < 20) {
113	ix &= ~(0x000fffff >> j0);
114	lx = 0;
115	} else {
116	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
117	}
118	INSERT_WORDS(x,ix,lx);
119
120	ax -= x;
121	if (ax == 0.5)
122	return (0);
123
124	ax = ax == 0 ? 1 : __compute_cospi(ax);
125
126	if (j0 > 30) x -= 0x1p30;
127	j0 = (uint32_t)x;
128
129	return (j0 & 1 ? -ax : ax);
130	}
131
132	if (ix >= 0x7f800000)
133	return (x - x);
134
135	/*
136	* \|x\| >= 0x1p52 is always an even integer, so return 1.
137	* FIXME: should this raise FE_INEXACT or FE_INVALID.
138	*/
139	return (1);
140	}
141
142	#if LDBL_MANT_DIG == 53
143	__weak_reference(cospi, cospil);
144	#endif

Line 0 Link Here

(-)src/s_cospif.c (+96 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	#include "math.h"
		32	#include "math_private.h"
		33	#include "k_cospif.c"
		34	#include "k_sinpif.c"
		35
		36	static inline float
		37	__compute_cospif(float x)
		38	{
		39
		40	if (x <= 0.5f) {
		41	if (x <= 0.25f)
		42	return (__kernel_cospif(x));
		43	return (__kernel_sinpif(0.5f - x));
		44	} else {
		45	if (x <= 0.75f)
		46	return (-__kernel_sinpif(x - 0.5f));
		47	return (-__kernel_cospif(1 - x));
		48	}
		49	}
		50
		51	float
		52	cospif(float x)
		53	{
		54	static const float huge = 1e30;
		55	float ax;
		56	uint32_t ix, j0;
		57
		58	GET_FLOAT_WORD(ix, x);
		59	ix = ix & 0x7fffffff;
		60	SET_FLOAT_WORD(ax, ix);
		61
		62	if (ix < 0x3f800000) { /* \|x\| < 1 */
		63	if (ix < 0x39000000) { /* \|x\| < 0x1p-13 */
		64	if (huge * ax == 0) /* Raise inexact iff != 0. */
65	return (1);
66	}
67	return (__compute_cospif(ax));
68	}
69
70	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
71	/* Determine integer part of ax. */
72	j0 = ((ix >> 23) & 0xff) - 0x7f;
73	ix &= ~(0x007fffff >> j0);
74	SET_FLOAT_WORD(x, ix);
75
76	ax -= x;
77	if (ax == 0.5f)
78	return (0);
79
80	ax = ax == 0 ? 1 : __compute_cospif(ax);
81
82	ix = (uint32_t)x;
83
84	return (ix & 1 ? -ax : ax);
85	}
86
87	/* x = +-inf or nan. */
88	if (ix >= 0x7f800000)
89	return (x - x);
90
91	/*
92	* \|x\| >= 0x1p23 is always an even integer, so return 1.
93	* FIXME: should this raise FE_INEXACT or FE_INVALID.
94	*/
95	return (1);
96	}

Line 0 Link Here

(-)src/s_sinpi.c (+151 lines)
		1	/*-
		2	* Copyright (c) 2017 Steven G. Kargl
		3	* All rights reserved.
		4	*
		5	* Redistribution and use in source and binary forms, with or without
		6	* modification, are permitted provided that the following conditions
		7	* are met:
		8	* 1. Redistributions of source code must retain the above copyright
		9	* notice unmodified, this list of conditions, and the following
		10	* disclaimer.
		11	* 2. Redistributions in binary form must reproduce the above copyright
		12	* notice, this list of conditions and the following disclaimer in the
		13	* documentation and/or other materials provided with the distribution.
		14	*
		15	* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
		16	* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
		17	* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
		18	* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
		19	* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
		20	* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
		21	* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
		22	* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
		23	* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
		24	* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
		25	*/
		26
		27	/**
		28	* 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 = 3 or 4. To achieve a max ULP less than
		36	* 0.67, pi is decomposed into high and low parts with the high part
		37	* containing a number of trailing zero bits. x is also split into high
		38	* and low parts.
		39	*
		40	* 2. For \|x\| < 1, argument reduction is not required and sinpi(x) is
		41	* computed by polynomial approximations for sinpi(x) and cospi(x) in the
		42	* the interval [0,0.25]. See k_sinpi.c and k_cospi.c for the details for
		43	* these polynomial approximations.
		44	*
		45	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
		46	* \|x\| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder. With
		47	* the given domain, a simplified inline j0 = floor(x) is used.
		48	*
		49	* sinpi(x) = sin(pi*(j0+r))
		50	* = sin(pij0) cos(pir) + cos(pij0) * sin(pi*r)
		51	* = cos(pij0) sin(pi*r)
		52	* = +-sinpi(r)
		53	*
		54	* If j0 is even, then cos(pij0) = 1. If j0 is odd, then cos(pij0) = -1.
		55	* sinpi(r) is then computed via polynomial approximation.
		56	*
		57	* 4. For \|x\| >= 0x1p(P-1), \|x\| is integral and sinpi(x) = copysign(0,x).
		58	*
		59	* 5. Special cases:
		60	*
		61	* sinpi(+-0) = +-0
		62	* sinpi(+-n) = +-0, for positive integers n.
		63	* sinpi(+-inf) = nan. Raises the "invalid" floating-point exception.
		64	* sinpi(nan) = nan. Raises the "invalid" floating-point exception.
65	*/
66
67	#include "math.h"
68	#include "math_private.h"
69	#include "k_cospi.c"
70	#include "k_sinpi.c"
71
72	static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
73	static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
74
75	/*
76	* To compute sinpi(x) for 0 < x < 1, use trignometric identities to map
77	* sinpi(x) into the [0,0.25] quadrant of sinpi(r) or cospi(r).
78	*/
79	static inline double
80	__compute_sinpi(double x)
81	{
82
83	if (x < 0.5) {
84	if (x <= 0.25)
85	return (__kernel_sinpi(x));
86	return (__kernel_cospi(0.5 - x));
87	} else {
88	if (x < 0.75)
89	return (__kernel_cospi(x - 0.5));
90	return (__kernel_sinpi(1 - x));
91	}
92	}
93
94	double
95	sinpi(double x)
96	{
97	double ax;
98	uint32_t hx, ix, j0, lx;
99
100	EXTRACT_WORDS(hx, lx, x);
101	ix = hx & 0x7fffffff;
102	INSERT_WORDS(ax, ix, lx);
103
104	if (ix < 0x3ff00000) { /* \|x\| < 1 */
105	if (ix < 0x3e300000) { /* \|x\| < 0x1p-28 */
106	if (x == 0)
107	return (x);
108	INSERT_WORDS(ax, hx, 0);
109	x -= ax;
110	return (pilo * x + pihi * x + pilo * ax + pihi * ax);
111	}
112	ax = __compute_sinpi(ax);
113	return ((hx & 0x80000000) ? -ax : ax);
114	}
115
116	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
117	/* Determine integer part of ax. */
118	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
119	if (j0 < 20) {
120	ix &= ~(0x000fffff >> j0);
121	lx = 0;
122	} else {
123	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
124	}
125	INSERT_WORDS(x,ix,lx);
126
127	ax -= x;
128	if (ax == 0) {
129	ax = 0;
130	} else {
131	ax = __compute_sinpi(ax);
132	if (j0 > 30) x -= 0x1p30;
133	lx = (uint32_t)x;
134	if (lx & 1) ax = -ax;
135	}
136	return ((hx & 0x80000000) ? -ax : ax);
137	}
138
139	if (ix >= 0x7f800000)
140	return (x - x);
141
142	/*
143	* \|x\| >= 0x1p52 is always an integer, so return +-0.
144	* FIXME: should this raise FE_INEXACT or FE_INVALID.
145	*/
146	return (copysign(0, x));
147	}
148
149	#if LDBL_MANT_DIG == 53
150	__weak_reference(sinpi, sinpil);
151	#endif

Line 0 Link Here

(-)src/s_sinpif.c (+102 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	#include "math.h"
		32	#include "math_private.h"
		33	#include "k_cospif.c"
		34	#include "k_sinpif.c"
		35
		36	static inline float
		37	__compute_sinpif(float x)
		38	{
		39
		40	if (x <= 0.5f) {
		41	if (x <= 0.25f)
		42	return (__kernel_sinpif(x));
		43	return (__kernel_cospif(0.5f - x));
		44	} else {
		45	if (x <= 0.75f)
		46	return (__kernel_cospif(x - 0.5f));
		47	return (__kernel_sinpif(1 - x));
		48	}
		49	}
		50
		51	static const float pihi = 3.14160156e+00f; /* 0x40491000 */
		52	static const float pilo = -8.90890988e-06f; /* 0xb715777a */
		53
		54	float
		55	sinpif(float x)
		56	{
		57	float ax;
		58	uint32_t hx, ix, j0;
		59
		60	GET_FLOAT_WORD(hx, x);
		61	ix = hx & 0x7fffffff;
		62	SET_FLOAT_WORD(ax, ix);
		63
		64	if (ix < 0x3f800000) { /* \|x\| < 1 */
65	if (ix < 0x38000000) { /* \|x\| < 0x1p-15 */
66	if (x == 0)
67	return (x);
68	SET_FLOAT_WORD(ax, (hx >> 13) << 13);
69	x -= ax;
70	return (pilo * x + pihi * x + pilo * ax + pihi * ax);
71	}
72	ax = __compute_sinpif(ax);
73	return ((hx & 0x80000000) ? -ax : ax);
74	}
75
76	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
77	/* Determine integer part of ax. */
78	j0 = ((ix >> 23) & 0xff) - 0x7f;
79	ix &= ~(0x007fffff >> j0);
80	SET_FLOAT_WORD(x, ix);
81
82	ax -= x;
83	if (ax == 0) {
84	ax = 0;
85	} else {
86	ix = (uint32_t)x;
87	ax = __compute_sinpif(ax);
88	if (ix & 1) ax = -ax;
89	}
90	return ((hx & 0x80000000) ? -ax : ax);
91	}
92
93	/* x = +-inf or nan. */
94	if (ix >= 0x7f800000)
95	return (x - x);
96
97	/*
98	* \|x\| >= 0x1p23 is always an integer, so return +-0.
99	* FIXME: should this raise FE_INEXACT or FE_INVALID.
100	*/
101	return (copysignf(0, x));
102	}

Line 0 Link Here

(-)src/s_tanpi.c (+240 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 = 3 or 4. To achieve a max ULP less than
		36	* 0.59, pi is decomposed into high and low parts with the high part
		37	* containing a number of trailing zero bits. x is also split into high
		38	* and low parts.
		39	*
		40	* 2. For \|x\| < 1, argument reduction is not required and tanpi(x) is
		41	* computed by either a polynomial approximation or the definition of
		42	* tan in terms of sin and cos. For \|x\| in [0,0.25), a polynomial
		43	* approximation is used. For \|x\| in (0.25, 0.5], tanpi(x) is computed
		44	* by tanpi(x) = cospi(0.5 - x) / sinpi(0.5 - x).
		45	*
		46	* 3. For 1 <= \|x\| < 0x1p(P-1), argument reduction is required where
		47	* \|x\| = j0 + r with j0 an integer and 0 <= r < 1 is a remainder. With
		48	* the given domain, a simplified inline j0 = floor(x) is used.
		49	*
		50	* tanpi(x) = tan(pi*(j0+r))
		51	* tan(pij0) + tan(pir)
		52	* = ----------------------------
		53	* 1 - tan(pij0) tan(pi*r)
		54	* = tanpi(r)
		55	*
		56	* 4. For \|x\| >= 0x1p(P-1), \|x\| is integral and tanpi(x) = copysign(0,x).
		57	*
		58	* 5. Special cases:
		59	*
		60	* tanpif(+-0) = +-0
		61	* tanpif(+-n) = +-0, for positive integers n.
		62	* tanpif(+-n+1/4) = +-1, for positive integers n.
		63	* tanpif(+-n+1/2) = NaN, for positive integers n.
		64	* tanpif(+-inf) = nan. Raises the "invalid" floating-point exception.
65	* tanpif(nan) = nan. Raises the "invalid" floating-point exception.
66	*
67	* Compute tanpi(x) = tan(pi*x) via a polynomial approximation. Coefficients
68	* were determined for the following function:
69	*
70	* f(x) = tan(pix) / x = t0 + t1x*2 + t2x**4 + ...
71	*
72	* with an Remes algorithm. The Remes algorithm used MPFR with a working
73	* precision of 30 times that of the target precision (e.g., for 53-bit
74	* precision MPFR computed f(x) with 1590-bit precision). To accumulate the
75	* result with sufficient accuracy the t0, t1, and t2 coefficients and final
76	* multiplications and additions needed to be handled specially.
77	*
78	* x2 = x * x
79	* T(x2) = t3 + t4x2 + t5x2**2 + ...
80	* = t3 + (t4 + (t5 ...) * x2) * x2 (Horner's method)
81	* tanpi(x) = x * (t0 + x2(t1 + x2(t2 + x2*T(x2))))
82	*
83	* The final multiplications and additions in the above expression for
84	* tanpi(x) are done in an extra precision arithmetic.
85	*
86	* On an Intel(R) Core(TM)2 Duo CPU T7250, limited testing in the
87	* interval [0x1p-28,0.25] gave
88	*
89	* a = 2.4250162020325661e-01
90	* tanpi(a) = 9.5396227802027667e-01 (0x3fee86db 0xe636df37)
91	* tan(pi*a) = 9.5396227802027678e-01 (0x3fee86db 0xe636df38)
92	*
93	* MAX ULP: 1.37582114
94	* Total tested: 134217727
95	* 1.0 < ULP : 309144
96	* 0.9 < ULP <= 1.0: 1047636
97	* 0.8 < ULP <= 0.9: 2732819
98	* 0.7 < ULP <= 0.8: 5491591
99	* 0.6 < ULP <= 0.7: 8639239
100	*
101	* In the interval [0.25,0.48], limited testing finds
102	*
103	* a = 4.1875240078199505e-01
104	* tanpi(a) = 3.8323217592887793e+00 (0x400ea898 0x4f7f27eb)
105	* tan(pi*a) = 3.8323217592887784e+00 (0x400ea898 0x4f7f27e9)
106	*
107	* MAX ULP: 2.01127812
108	* Total tested: 134217727
109	* 1.0 < ULP : 5702704
110	* 0.9 < ULP <= 1.0: 3602502
111	* 0.8 < ULP <= 0.9: 5207602
112	* 0.7 < ULP <= 0.8: 7173276
113	* 0.6 < ULP <= 0.7: 9449882
114	*/
115
116	#include "math.h"
117	#include "math_private.h"
118	#include "k_cospi.c"
119	#include "k_sinpi.c"
120
121	static const double pihi = 3.1415926814079285e+00; /* 0x400921fb 0x58000000 */
122	static const double pilo =-2.7818135228334233e-08; /* 0xbe5dde97 0x3dcb3b3a */
123
124	static const double
125	t00 = 3.1415926535897931e+00, /* 0x400921fb 0x54442d18 */
126	t00lo = 1.2267625636978282e-16, /* 0x3ca1adf4 0x7b8ee56b */
127	t01 = 1.0335425560099939e+01, /* 0x4024abbc 0xe625be52 */
128	t01lo = -6.0289590484334045e-16, /* 0xbcc5b8bb 0xb4f3d850 */
129	t02 = 4.0802624638040442e+01, /* 0x404466bc 0x6775ac7c */
130	t02lo = -3.1181465177803440e-15, /* 0xbcec15f4 0xd355f491 */
131	t03 = 1.6299995197351416e+02, /* 0x40645fff 0x9b47fa0c */
132	t04 = 6.5190975669363422e+02, /* 0x40845f47 0x2e8473cf */
133	t05 = 2.6075994007466215e+03, /* 0x40a45f32 0xe4a797e0 */
134	t06 = 1.0430393632420044e+04, /* 0x40c45f32 0x628c115e */
135	t07 = 4.1720374420664804e+04, /* 0x40e45f0b 0xfb410bc9 */
136	t08 = 1.6695720713150888e+05, /* 0x41046169 0xa8349085 */
137	t09 = 6.6429163117898861e+05, /* 0x412445c7 0x4329e474 */
138	t10 = 2.7802499382503941e+06, /* 0x4145362c 0xf81896c3 */
139	t11 = 7.9185101944778729e+06, /* 0x415e34eb 0x8c725352 */
140	t12 = 9.3691291294289947e+07, /* 0x41965676 0x6d2d5a58 */
141	t13 = -5.0585651903137898e+08, /* 0xc1be26c2 0x07080874 */
142	t14 = 6.8772398519624548e+09, /* 0x41f99ea5 0xa2bf6637 */
143	t15 = -3.3094464588263161e+10, /* 0xc21ed255 0xd1310d7a */
144	t16 = 1.1673794006564838e+11; /* 0x423b2e1f 0x9a61a5fc */
145
146	static inline double
147	__kernel_tanpi(double x)
148	{
149	double hi, lo, t, x2, xhi, xlo;
150	uint32_t hx, lx;
151
152	if (x < 0.25) {
153	x2 = x * x;
154	t = x2 * (x2 * (x2 * (t15 + x2 * t16) + t14) + t13) + t12;
155	t = x2 * (x2 * (x2 * (x2 * t + t11) + t10) + t09) + t08;
156	t = x2 * (x2 * (x2 * (x2 * t + t07) + t06) + t05) + t04;
157
158	EXTRACT_WORDS(hx, lx, x);
159	INSERT_WORDS(xhi, hx, 0);
160	xlo = x - xhi;
161	xlo *= (xlo + xhi + xhi);
162	xhi *= xhi;
163
164	t = t02lo + x2 * (x2 * t + t03) + t02;
165	t = t01lo + xlo * t + xhi * t + t01;
166	t = t00lo + xlo * t + xhi * t + t00;
167	return (t * x);
168	} else if (x > 0.25) {
169	x = 0.5 - x;
170	t = __kernel_cospi(x) / __kernel_sinpi(x);
171	} else
172	t = 1;
173	return (t);
174	}
175
176	static inline double
177	__compute_tanpi(double x)
178	{
179
180	return (x < 0.5 ? __kernel_tanpi(x) : -__kernel_tanpi(1 - x));
181	}
182
183	double
184	tanpi(double x)
185	{
186	double ax;
187	uint32_t hx, ix, j0, lx;
188
189	EXTRACT_WORDS(hx, lx, x);
190	ix = hx & 0x7fffffff;
191	INSERT_WORDS(ax, ix, lx);
192
193	if (ix < 0x3ff00000) { /* \|x\| < 1 */
194	if (ix < 0x3e200000) { /* \|x\| < 0x1p-29 */
195	if (x == 0)
196	return (x);
197	INSERT_WORDS(ax, hx, 0);
198	x -= ax;
199	return (pilo * x + pihi * x + pilo * ax + pihi * ax);
200	}
201	if (ix == 0x3fe00000 && !lx)
202	return ((x - x) / (x - x));
203	ax = __compute_tanpi(ax);
204	return ((hx & 0x80000000) ? -ax : ax);
205	}
206
207	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
208	/* Determine integer part of ax. */
209	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
210	if (j0 < 20) {
211	ix &= ~(0x000fffff >> j0);
212	lx = 0;
213	} else {
214	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
215	}
216	INSERT_WORDS(x,ix,lx);
217
218	ax -= x;
219
220	if (ax == 0.5)
221	return ((x - x) / (x - x));
222
223	ax = ax == 0 ? 0 : __compute_tanpi(ax);
224	return ((hx & 0x80000000) ? -ax : ax);
225	}
226
227	/* x = +-inf or nan. */
228	if (ix >= 0x7f800000)
229	return (x - x);
230
231	/*
232	* \|x\| >= 0x1p53 is always an integer, so return +-0.
233	* FIXME: should this raise FE_INEXACT or FE_INVALID.
234	*/
235	return (copysign(0, x));
236	}
237
238	#if LDBL_MANT_DIG == 53
239	__weak_reference(tanpi, tanpil);
240	#endif

Line 0 Link Here

(-)src/s_tanpif.c (+138 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	#include "math.h"
		32	#include "math_private.h"
		33	#include "k_cospif.c"
		34	#include "k_sinpif.c"
		35
		36	static const float pihi = 3.14160156e+00f; /* 0x40491000 */
		37	static const float pilo = -8.90890988e-06f; /* 0xb715777a */
		38
		39	static const float
		40	t0 = 3.14159274e+00f, /* 0x40490fdb */
		41	t0lo = -8.71231194e-08f, /* 0xb3bb1871 */
		42	t1 = 1.03354244e+01f, /* 0x41255de6 */
		43	t1lo = 3.83117083e-07f, /* 0x34cdaf36 */
		44	t2 = 4.08029366e+01f, /* 0x42233635 */
		45	t2lo = -1.28256147e-07f, /* 0xb409b6c8 */
		46	t3 = 1.62950455e+02f, /* 0x4322f351 */
		47	t3lo = 4.45962769e-06f, /* 0x3695a3e9 */
		48	t4 = 6.55820984e+02f, /* 0x4423f48b */
		49	t5 = 2.43586987e+03f, /* 0x45183deb */
		50	t6 = 1.47717539e+04f, /* 0x4666cf04 */
		51	t7 = -1.96741250e+04f, /* 0xc699b440 */
		52	t8 = 5.87398688e+05f; /* 0x490f686b */
		53
		54	static inline float
		55	__kernel_tanpif(float x)
		56	{
		57	float hi, lo, t, x2, xhi, xlo;
		58	uint32_t ix;
		59
		60	if (x < 0.25f) { /* 0 <= \|x\| < 0.25 */
		61	x2 = x * x;
		62	t = t3lo + x2 * (t4 + x2 * (t5 + x2 * (t6 + x2 * (t7
		63	+ x2 * t8)))) + t3;
		64
65	GET_FLOAT_WORD(ix, x);
66	SET_FLOAT_WORD(xhi, (ix >> 14) << 14);
67	xlo = x - xhi;
68	xlo *= (xlo + xhi + xhi);
69	xhi *= xhi;
70
71	t = t2lo + xlo * t + xhi * t + t2;
72	t = t1lo + xlo * t + xhi * t + t1;
73	t = t0lo + xlo * t + xhi * t + t0;
74	t *= x;
75	} else if (x > 0.25f) { /* 0.25 < \|x\| < 0.5 */
76	x = 0.5f - x;
77	t = __kernel_cospif(x) / __kernel_sinpif(x);
78	} else
79	t = 1;
80	return (t);
81	}
82
83	static inline float
84	__compute_tanpif(float x)
85	{
86
87	return (x < 0.5f ? __kernel_tanpif(x) : -__kernel_tanpif(1 - x));
88	}
89
90	float
91	tanpif(float x)
92	{
93	float ax;
94	uint32_t hx, ix, j0;
95
96	GET_FLOAT_WORD(hx, x);
97	ix = hx & 0x7fffffff;
98	SET_FLOAT_WORD(ax, ix);
99
100	if (ix < 0x3f800000) { /* \|x\| < 1 */
101	if (ix < 0x38800000) { /* \|x\| < 0x1p-14 */
102	if (x == 0)
103	return (x);
104	SET_FLOAT_WORD(ax, (hx >> 13) << 13);
105	x -= ax;
106	return (pilo * x + pihi * x + pilo * ax + pihi * ax);
107	}
108	if (ix == 0x3f000000)
109	return ((x - x) / (x - x));
110	ax = __compute_tanpif(ax);
111	return ((hx & 0x80000000) ? -ax : ax);
112	}
113
114	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
115	/* Determine integer part of ax. */
116	j0 = ((ix >> 23) & 0xff) - 0x7f;
117	ix &= ~(0x007fffff >> j0);
118	SET_FLOAT_WORD(x, ix);
119
120	ax -= x;
121
122	if (ax == 0.5f)
123	return ((x - x) / (x - x));
124
125	ax = ax == 0 ? 0 : __compute_tanpif(ax);
126	return ((hx & 0x80000000) ? -ax : ax);
127	}
128
129	/* x = +-inf or nan. */
130	if (ix >= 0x7f800000)
131	return (x - x);
132
133	/*
134	* \|x\| >= 0x1p23 is always an integer, so return +-0.
135	* FIXME: should this raise FE_INEXACT or FE_INVALID.
136	*/
137	return (copysignf(0, x));
138	}