Attachment #181395 for bug #218300




	s_nexttowardf.c s_remquo.c s_remquof.c \
	s_rint.c s_rintf.c s_round.c s_roundf.c \
	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
	s_signgam.c s_significand.c s_significandf.c s_sin.c \
	s_sincos.c s_sincosf.c s_sinf.c \
	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
	w_cabs.c w_cabsf.c w_drem.c w_dremf.c


	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
	s_scalbnl.c s_sinl.c s_sincosl.c \
	s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
.endif

# C99 complex functions

	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
	nextafter.3 remainder.3 rint.3 \
	round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
	sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
	complex.3

MLINKS+=acos.3 acosf.3 acos.3 acosl.3

MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
MLINKS+=sin.3 sinf.3 sin.3 sinl.3
MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
	sqrt.3 sqrtl.3

Lines 294-297 Link Here

(-)Symbol.map (+3 lines)
294	casinl;	294	casinl;
295	catanl;	295	catanl;
296	catanhl;	296	catanhl;
		297	sincos;
		298	sincosf;
		299	sincosl;
297	};	300	};




.\" Copyright (c) 2011 Steven G. Kargl.
.\"
.\" 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.
.\"
.\"
.Dd March 12, 2011
.Dt SINCOS 3
.Os
.Sh NAME
.Nm sincos ,
.Nm sincosf ,
.Nm sincosl
.Nd sine and cosine functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In math.h
.Ft void
.Fn sincos "double x" "double *s" "double *c"
.Ft void
.Fn sincosf "float x" "float *s" "float *c"
.Ft void
.Fn sincosl "long double x" "long double *s" "long double *c"
.Sh DESCRIPTION
The
.Fn sincos ,
.Fn sincosf ,
and
.Fn sincosl
functions compute the sine and cosine of
.Fa x .
Using these functions allows argument reduction to occur only
once instead of twice with individual invocations of 
.Fn sin
and 
.Fn cos .
Like 
.Fn sin
and 
.Fn cos ,
a large magnitude argument may yield a result with little
or no significance.
.Sh RETURN VALUES
Upon returning from 
.Fn sincos ,
.Fn sincosf ,
and
.Fn sincosl ,
the memory pointed to by 
.Ar "*s" 
and
.Ar "*c" 
are assigned the values of sine and cosine, respectively.
.Sh SEE ALSO
.Xr cos 3 ,
.Xr sin 3 ,
.Sh HISTORY
These functions were added to 
.Fx 9.0
to aid in writing various complex function contained in 
.St -isoC-99 .


Line 0 Link Here

(-)src/k_sincos.h (+49 lines)
	1	/*
	2	* ====================================================
	3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	4	*
	5	* Developed at SunSoft, a Sun Microsystems, Inc. business.
	6	* Permission to use, copy, modify, and distribute this
	7	* software is freely granted, provided that this notice
	8	* is preserved.
	9	* ====================================================
	10	*
	11	* k_sin.c and k_cos.c merged by Steven G. Kargl.
	12	*/
	13
	14	static const double
	15	S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
	16	S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
	17	S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
	18	S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
	19	S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
	20	S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
	21
	22	static const double
	23	C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
	24	C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
	25	C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
	26	C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
	27	C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
	28	C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
	29
	30	static inline void
	31	__kernel_sincos(double x, double y, int iy, double sn, double cs)
	32	{
	33	double hz, r, v, w, z;
	34
	35	z = x * x;
	36	w = z * z;
	37	r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
	38	v = z * x;
	39
	40	if (iy == 0)
	41	sn = x + v (S1 + z * r);
	42	else
	43	sn = x - ((z (y / 2 - v * r) - y) - v * S1);
	44
	45	r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
	46	hz = z / 2;
	47	w = 1 - hz;
	48	cs = w + (((1 - w) - hz) + (z r - x * y));
	49	}

Line 0 Link Here

(-)src/k_sincosf.h (+40 lines)
	1	/*
	2	* ====================================================
	3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
	4	*
	5	* Developed at SunPro, a Sun Microsystems, Inc. business.
	6	* Permission to use, copy, modify, and distribute this
	7	* software is freely granted, provided that this notice
	8	* is preserved.
	9	* ====================================================
	10	*
	11	* k_sinf.c and k_cosf.c merged by Steven G. Kargl.
	12	*/
	13
	14	/* \|sin(x)/x - s(x)\| < 2*-37.5 (~[-4.89e-12, 4.824e-12]). /
	15	static const double
	16	S1 = -0x15555554cbac77.0p-55, /* -0.166666666416265235595 */
	17	S2 = 0x111110896efbb2.0p-59, /* 0.0083333293858894631756 */
	18	S3 = -0x1a00f9e2cae774.0p-65, /* -0.000198393348360966317347 */
	19	S4 = 0x16cd878c3b46a7.0p-71; /* 0.0000027183114939898219064 */
	20
	21	/* \|cos(x) - c(x)\| < 2*-34.1 (~[-5.37e-11, 5.295e-11]). /
	22	static const double
	23	C0 = -0x1ffffffd0c5e81.0p-54, /* -0.499999997251031003120 */
	24	C1 = 0x155553e1053a42.0p-57, /* 0.0416666233237390631894 */
	25	C2 = -0x16c087e80f1e27.0p-62, /* -0.00138867637746099294692 */
	26	C3 = 0x199342e0ee5069.0p-68; /* 0.0000243904487962774090654 */
	27
	28	static inline void
	29	__kernel_sincosdf(double x, float sn, float cs)
	30	{
	31	double r, s, w, z;
	32
	33	z = x * x;
	34	w = z * z;
	35	r = S3 + z * S4;
	36	s = z * x;
	37	sn = (x + s (S1 + z * S2)) + s * w * r;
	38	r = C2 + z * C3;
	39	cs = ((1 + z C0) + w * C1) + (w * z) * r;
	40	}

Line 0 Link Here

(-)src/k_sincosl.h (+131 lines)
		1	/*
		2	* ====================================================
		3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		4	* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
		5	*
		6	* Developed at SunSoft, a Sun Microsystems, Inc. business.
		7	* Permission to use, copy, modify, and distribute this
		8	* software is freely granted, provided that this notice
		9	* is preserved.
		10	* ====================================================
		11	*
		12	* k_sinl.c and k_cosl.c merged by Steven G. Kargl
		13	*/
		14
		15	#if LDBL_MANT_DIG == 64 /* ld80 version of k_sincosl.c. */
		16
		17	#if defined(__amd64__) \|\| defined(__i386__)
		18	/* Long double constants are slow on these arches, and broken on i386. */
		19	static const volatile double
		20	C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */
		21	C1lo = 2.2598839032744733e-18, /* 0x14d80000000000.0p-111 */
		22	S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */
		23	S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */
		24	#define S1 ((long double)S1hi + S1lo)
		25	#define C1 ((long double)C1hi + C1lo)
		26	#else
		27	static const long double
		28	C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
		29	S1 = -0.166666666666666666671L, /* -0xaaaaaaaaaaaaaaab.0p-66 */
		30	#endif
		31
		32	static const double
		33	C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */
		34	C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */
		35	C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */
		36	C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */
		37	C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */
		38	C7 = 4.7383039476436467e-14, /* 0x1aac9d9af5c43e.0p-97 */
		39	S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
		40	S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
		41	S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
		42	S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
		43	S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
		44	S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
		45	S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
		46
		47	static inline void
		48	__kernel_sincosl(long double x, long double y, int iy, long double *sn,
		49	long double *cs)
		50	{
		51	long double hz, r, v, w, z;
		52
		53	z = x * x;
		54	v = z * x;
		55	/*
		56	* XXX Replace Horner scheme with an algorithm suitable for CPUs
		57	* with more complex pipelines.
		58	*/
		59	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8)))));
		60
		61	if (iy == 0)
		62	sn = x + v (S1 + z * r);
		63	else
		64	sn = x - ((z (y / 2 - v * r) - y) - v * S1);
65
66	hz = z / 2;
67	w = 1 - hz;
68	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
69	z * C7))))));
70	cs = w + (((1 - w) - hz) + (z r - x * y));
71	}
72
73	#elif LDBL_MANT_DIG == 113 /* ld128 version of k_sincosl.c. */
74
75	static const long double
76	C1 = 0.04166666666666666666666666666666658424671L,
77	C2 = -0.001388888888888888888888888888863490893732L,
78	C3 = 0.00002480158730158730158730158600795304914210L,
79	C4 = -0.2755731922398589065255474947078934284324e-6L,
80	C5 = 0.2087675698786809897659225313136400793948e-8L,
81	C6 = -0.1147074559772972315817149986812031204775e-10L,
82	C7 = 0.4779477332386808976875457937252120293400e-13L,
83	S1 = -0.16666666666666666666666666666666666606732416116558L,
84	S2 = 0.0083333333333333333333333333333331135404851288270047L,
85	S3 = -0.00019841269841269841269841269839935785325638310428717L,
86	S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
87	S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
88	S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
89	S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
90	S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
91
92	static const double
93	C8 = -0.1561920696721507929516718307820958119868e-15,
94	C9 = 0.4110317413744594971475941557607804508039e-18,
95	C10 = -0.8896592467191938803288521958313920156409e-21,
96	C11 = 0.1601061435794535138244346256065192782581e-23,
97	S9 = -0.82206352458348947812512122163446202498005154296863e-17,
98	S10 = 0.19572940011906109418080609928334380560135358385256e-19,
99	S11 = -0.38680813379701966970673724299207480965452616911420e-22,
100	S12 = 0.64038150078671872796678569586315881020659912139412e-25;
101
102	static inline void
103	__kernel_sincosl(long double x, long double y, int iy, long double *sn,
104	long double *cs)
105	{
106	long double hz, r, v, w, z;
107
108	z = x * x;
109	v = z * x;
110	/*
111	* XXX Replace Horner scheme with an algorithm suitable for CPUs
112	* with more complex pipelines.
113	*/
114	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
115	z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
116
117	if (iy == 0)
118	sn = x + v (S1 + z * r);
119	else
120	cs = x - ((z (y / 2 - v * r) - y) - v * S1);
121
122	hz = z / 2;
123	w = 1 - hz;
124	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
125	z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
126
127	cs = w + (((1 - w) - hz) + (z r - x * y));
128	}
129	#else
130	#error "Unsupported long double format"
131	#endif

Lines 500-505 Link Here

(-)src/math.h (+3 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	void sincos(double, double , double );
		504	void sincosf(float, float , float );
		505	void sincosl(long double, long double , long double );
503	#endif	506	#endif
504		507
505	__END_DECLS	508	__END_DECLS

Lines 306-314 Link Here

(-)src/math_private.h (+6 lines)
306	fpsetprec(__oprec); \	306	fpsetprec(__oprec); \
307	RETURNF(__retval); \	307	RETURNF(__retval); \
308	} while (0)	308	} while (0)
		309	#define RETURNV() do { \
		310	if (__oprec != FP_PE) \
		311	fpsetprec(__oprec); \
		312	return; \
		313	} while (0)
309	#else	314	#else
310	#define ENTERI(x)	315	#define ENTERI(x)
311	#define RETURNI(x) RETURNF(x)	316	#define RETURNI(x) RETURNF(x)
		317	#define RETURNV() return
312	#endif	318	#endif
313		319
314	/* Default return statement if hack_t() is not used. /	320	/* Default return statement if hack_t() is not used. /




/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 * 
 * s_sin.c and s_cos.c merged by Steven G. Kargl.  Descriptions of the
 * algorithms are contained in the original files.
 */

#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

#include <float.h>

#include "math.h"
#define INLINE_REM_PIO2
#include "math_private.h"
#include "e_rem_pio2.c"
#include "k_sincos.h"

void
sincos(double x, double *sn, double *cs)
{
	double y[2];
	int32_t n, ix;

	/* High word of x. */
	GET_HIGH_WORD(ix, x);

	/* |x| ~< pi/4 */
	ix &= 0x7fffffff;
	if (ix <= 0x3fe921fb) {
		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
			if ((int)x == 0) {	/* Generate inexact. */
				*sn = x;
				*cs = 1;
				return;
			}
		}
		__kernel_sincos(x, 0, 0, sn, cs);
		return;
	}

	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
	if (ix >= 0x7ff00000) {
		*sn = x - x;
		*cs = x - x;
		return;
	}

	/* Argument reduction. */
	n = __ieee754_rem_pio2(x, y);

	switch(n & 3) {
	case 0:
		__kernel_sincos(y[0], y[1], 1, sn, cs);
		break;
	case 1:
		__kernel_sincos(y[0], y[1], 1, cs, sn);
		*cs = -*cs;
		break;
	case 2:
		__kernel_sincos(y[0], y[1], 1, sn, cs);
		*sn = -*sn;
		*cs = -*cs;
		break;
	default:
		__kernel_sincos(y[0], y[1], 1, cs, sn);
		*sn = -*sn;
	}
}

#if (LDBL_MANT_DIG == 53)
__weak_reference(sincos, sincosl);
#endif




/* s_sincosf.c -- float version of s_sincos.c.
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 * Optimized by Bruce D. Evans.
 * Merged s_sinf.c and s_cosf.c by Steven G. Kargl.
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

#include <float.h>

#include "math.h"
#define INLINE_REM_PIO2F
#include "math_private.h"
#include "e_rem_pio2f.c"
#include "k_sincosf.h"

/* Small multiples of pi/2 rounded to double precision. */
static const double
p1pio2 = 1*M_PI_2,			/* 0x3FF921FB, 0x54442D18 */
p2pio2 = 2*M_PI_2,			/* 0x400921FB, 0x54442D18 */
p3pio2 = 3*M_PI_2,			/* 0x4012D97C, 0x7F3321D2 */
p4pio2 = 4*M_PI_2;			/* 0x401921FB, 0x54442D18 */

void
sincosf(float x, float *sn, float *cs)
{
	float c, s;
	double y;
	int32_t n, hx, ix;

	GET_FLOAT_WORD(hx, x);
	ix = hx & 0x7fffffff;

	if (ix <= 0x3f490fda) {		/* |x| ~<= pi/4 */
		if (ix < 0x39800000) {	/* |x| < 2**-12 */
			if ((int)x == 0) {
				*sn = x;	/* x with inexact if x != 0 */
				*cs = 1;
				return;
			}
		}
		__kernel_sincosdf(x, sn, cs);
		return;
	}

	if (ix <= 0x407b53d1) {		/* |x| ~<= 5*pi/4 */
		if (ix <= 0x4016cbe3) {	/* |x| ~<= 3pi/4 */
			if (hx > 0) {
				__kernel_sincosdf(x - p1pio2, cs, sn);
				*cs = -*cs;
			} else {
				__kernel_sincosdf(x + p1pio2, cs, sn);
				*sn = -*sn;
			}
		} else {
			if (hx > 0)
				__kernel_sincosdf(x - p2pio2, sn, cs);
			else
				__kernel_sincosdf(x + p2pio2, sn, cs);
			*sn = -*sn;
			*cs = -*cs;
		}
		return;
	}

	if (ix <= 0x40e231d5) {		/* |x| ~<= 9*pi/4 */
		if (ix <= 0x40afeddf) {	/* |x| ~<= 7*pi/4 */
			if (hx > 0) {
				__kernel_sincosdf(x - p3pio2, cs, sn);
				*sn = -*sn;
			} else {
				__kernel_sincosdf(x + p3pio2, cs, sn);
				*cs = -*cs;
			}
		} else {
			if (hx > 0)
				__kernel_sincosdf(x - p4pio2, sn, cs);
			else
				__kernel_sincosdf(x + p4pio2, sn, cs);
		}
		return;
	}

	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
	if (ix >= 0x7f800000) {
		*sn = x - x;
		*cs = x - x;
		return;
	}

	/* Argument reduction. */
	n = __ieee754_rem_pio2f(x, &y);
	__kernel_sincosdf(y, &s, &c);

	switch(n & 3) {
	case 0:
		*sn = s;
		*cs = c;
		break;
	case 1:
		*sn = c;
		*cs = -s;
		break;
	case 2:
		*sn = -s;
		*cs = -c;
		break;
	default:
		*sn = -c;
		*cs = s;
	}
}






/*-
 * Copyright (c) 2007, 2010-2013 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.
 *
 * s_sinl.c and s_cosl.c merged by Steven G. Kargl.
 */

#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

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

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

#if LDBL_MANT_DIG == 64
#include "../ld80/e_rem_pio2l.h"
#elif LDBL_MANT_DIG == 113
#include "../ld128/e_rem_pio2l.h"
#else
#error "Unsupported long double format"
#endif

void
sincosl(long double x, long double *sn, long double *cs)
{
	union IEEEl2bits z;
	int e0, sgn;
	long double y[2];

	z.e = x;
	sgn = z.bits.sign;
	z.bits.sign = 0;

	ENTERI();

	/* Optimize the case where x is already within range. */
	if (z.e < M_PI_4) {
		/*
		 * If x = +-0 or x is a subnormal number, then sin(x) = x and
		 * cos(x) = 1.
		 */
		if (z.bits.exp == 0) {
			*sn = x;
			*cs = 1;
		} else
			__kernel_sincosl(x, 0, 0, sn, cs);
		RETURNV();
	}

	/* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
	if (z.bits.exp == 32767) {
		*sn = x - x;
		*cs = x - x;
		RETURNV();
	}

	/* Range reduction. */
	e0 = __ieee754_rem_pio2l(x, y);

	switch (e0 & 3) {
	case 0:
		__kernel_sincosl(y[0], y[1], 1, sn, cs);
		break;
	case 1:
		__kernel_sincosl(y[0], y[1], 1, cs, sn);
		*cs = -*cs;
		break;
	case 2:
		__kernel_sincosl(y[0], y[1], 1, sn, cs);
		*sn = -*sn;
		*cs = -*cs;
		break;
	default:
		__kernel_sincosl(y[0], y[1], 1, cs, sn);
		*sn = -*sn;
	}

	RETURNV();
}

Return to bug 218300

Lines 73-79 Link Here

(-)Makefile (-3 / +7 lines)
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_sincos.c s_sincosf.c s_sinf.c \
77	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \	78	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
78	w_cabs.c w_cabsf.c w_drem.c w_dremf.c	79	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
79		80
Lines 104-110 Link Here
104	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \	105	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 \	106	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 \	107	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	108	s_scalbnl.c s_sinl.c s_sincosl.c \
		109	s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
108	.endif	110	.endif
109		111
110	# C99 complex functions	112	# C99 complex functions
Lines 137-143 Link Here
137	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \	139	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 \	140	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
139	nextafter.3 remainder.3 rint.3 \	141	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 \	142	round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
		143	sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
141	complex.3	144	complex.3
142		145
143	MLINKS+=acos.3 acosf.3 acos.3 acosl.3	146	MLINKS+=acos.3 acosf.3 acos.3 acosl.3
Lines 215-220 Link Here
215	MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3	218	MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
216	MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3	219	MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
217	MLINKS+=sin.3 sinf.3 sin.3 sinl.3	220	MLINKS+=sin.3 sinf.3 sin.3 sinl.3
		221	MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
218	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3	222	MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
219	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \	223	MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
220	sqrt.3 sqrtl.3	224	sqrt.3 sqrtl.3

Line 0 Link Here

(-)man/sincos.3 (+81 lines)
		1	.\" Copyright (c) 2011 Steven G. Kargl.
		2	.\"
		3	.\" Redistribution and use in source and binary forms, with or without
		4	.\" modification, are permitted provided that the following conditions
		5	.\" are met:
		6	.\" 1. Redistributions of source code must retain the above copyright
		7	.\" notice, this list of conditions and the following disclaimer.
		8	.\" 2. Redistributions in binary form must reproduce the above copyright
		9	.\" notice, this list of conditions and the following disclaimer in the
		10	.\" documentation and/or other materials provided with the distribution.
		11	.\"
		12	.\" THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
		13	.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		14	.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		15	.\" ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
		16	.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		17	.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		18	.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		19	.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		20	.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		21	.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		22	.\" SUCH DAMAGE.
		23	.\"
		24	.\"
		25	.Dd March 12, 2011
		26	.Dt SINCOS 3
		27	.Os
		28	.Sh NAME
		29	.Nm sincos ,
		30	.Nm sincosf ,
		31	.Nm sincosl
		32	.Nd sine and cosine functions
		33	.Sh LIBRARY
		34	.Lb libm
		35	.Sh SYNOPSIS
		36	.In math.h
		37	.Ft void
		38	.Fn sincos "double x" "double s" "double c"
		39	.Ft void
		40	.Fn sincosf "float x" "float s" "float c"
		41	.Ft void
		42	.Fn sincosl "long double x" "long double s" "long double c"
		43	.Sh DESCRIPTION
		44	The
		45	.Fn sincos ,
		46	.Fn sincosf ,
		47	and
		48	.Fn sincosl
		49	functions compute the sine and cosine of
		50	.Fa x .
		51	Using these functions allows argument reduction to occur only
		52	once instead of twice with individual invocations of
		53	.Fn sin
		54	and
		55	.Fn cos .
		56	Like
		57	.Fn sin
		58	and
		59	.Fn cos ,
		60	a large magnitude argument may yield a result with little
		61	or no significance.
		62	.Sh RETURN VALUES
		63	Upon returning from
		64	.Fn sincos ,
65	.Fn sincosf ,
66	and
67	.Fn sincosl ,
68	the memory pointed to by
69	.Ar "*s"
70	and
71	.Ar "*c"
72	are assigned the values of sine and cosine, respectively.
73	.Sh SEE ALSO
74	.Xr cos 3 ,
75	.Xr sin 3 ,
76	.Sh HISTORY
77	These functions were added to
78	.Fx 9.0
79	to aid in writing various complex function contained in
80	.St -isoC-99 .
81

Line 0 Link Here

(-)src/s_sincos.c (+80 lines)
		1	/*
		2	* ====================================================
		3	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		4	*
		5	* Developed at SunPro, a Sun Microsystems, Inc. business.
		6	* Permission to use, copy, modify, and distribute this
		7	* software is freely granted, provided that this notice
		8	* is preserved.
		9	* ====================================================
		10	*
		11	* s_sin.c and s_cos.c merged by Steven G. Kargl. Descriptions of the
		12	* algorithms are contained in the original files.
		13	*/
		14
		15	#include <sys/cdefs.h>
		16	__FBSDID("$FreeBSD$");
		17
		18	#include <float.h>
		19
		20	#include "math.h"
		21	#define INLINE_REM_PIO2
		22	#include "math_private.h"
		23	#include "e_rem_pio2.c"
		24	#include "k_sincos.h"
		25
		26	void
		27	sincos(double x, double sn, double cs)
		28	{
		29	double y[2];
		30	int32_t n, ix;
		31
		32	/* High word of x. */
		33	GET_HIGH_WORD(ix, x);
		34
		35	/* \|x\| ~< pi/4 */
		36	ix &= 0x7fffffff;
		37	if (ix <= 0x3fe921fb) {
		38	if (ix < 0x3e400000) { /* \|x\| < 2*-27 /
		39	if ((int)x == 0) { /* Generate inexact. */
		40	*sn = x;
		41	*cs = 1;
		42	return;
		43	}
		44	}
		45	__kernel_sincos(x, 0, 0, sn, cs);
		46	return;
		47	}
		48
		49	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
		50	if (ix >= 0x7ff00000) {
		51	*sn = x - x;
		52	*cs = x - x;
		53	return;
		54	}
		55
		56	/* Argument reduction. */
		57	n = __ieee754_rem_pio2(x, y);
		58
		59	switch(n & 3) {
		60	case 0:
		61	__kernel_sincos(y[0], y[1], 1, sn, cs);
		62	break;
		63	case 1:
		64	__kernel_sincos(y[0], y[1], 1, cs, sn);
65	cs = -cs;
66	break;
67	case 2:
68	__kernel_sincos(y[0], y[1], 1, sn, cs);
69	sn = -sn;
70	cs = -cs;
71	break;
72	default:
73	__kernel_sincos(y[0], y[1], 1, cs, sn);
74	sn = -sn;
75	}
76	}
77
78	#if (LDBL_MANT_DIG == 53)
79	__weak_reference(sincos, sincosl);
80	#endif

Line 0 Link Here

(-)src/s_sincosf.c (+126 lines)
		1	/* s_sincosf.c -- float version of s_sincos.c.
		2	* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
		3	* Optimized by Bruce D. Evans.
		4	* Merged s_sinf.c and s_cosf.c by Steven G. Kargl.
		5	*/
		6
		7	/*
		8	* ====================================================
		9	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		10	*
		11	* Developed at SunPro, a Sun Microsystems, Inc. business.
		12	* Permission to use, copy, modify, and distribute this
		13	* software is freely granted, provided that this notice
		14	* is preserved.
		15	* ====================================================
		16	*/
		17
		18	#include <sys/cdefs.h>
		19	__FBSDID("$FreeBSD$");
		20
		21	#include <float.h>
		22
		23	#include "math.h"
		24	#define INLINE_REM_PIO2F
		25	#include "math_private.h"
		26	#include "e_rem_pio2f.c"
		27	#include "k_sincosf.h"
		28
		29	/* Small multiples of pi/2 rounded to double precision. */
		30	static const double
		31	p1pio2 = 1M_PI_2, / 0x3FF921FB, 0x54442D18 */
		32	p2pio2 = 2M_PI_2, / 0x400921FB, 0x54442D18 */
		33	p3pio2 = 3M_PI_2, / 0x4012D97C, 0x7F3321D2 */
		34	p4pio2 = 4M_PI_2; / 0x401921FB, 0x54442D18 */
		35
		36	void
		37	sincosf(float x, float sn, float cs)
		38	{
		39	float c, s;
		40	double y;
		41	int32_t n, hx, ix;
		42
		43	GET_FLOAT_WORD(hx, x);
		44	ix = hx & 0x7fffffff;
		45
		46	if (ix <= 0x3f490fda) { /* \|x\| ~<= pi/4 */
		47	if (ix < 0x39800000) { /* \|x\| < 2*-12 /
		48	if ((int)x == 0) {
		49	sn = x; / x with inexact if x != 0 */
		50	*cs = 1;
		51	return;
		52	}
		53	}
		54	__kernel_sincosdf(x, sn, cs);
		55	return;
		56	}
		57
		58	if (ix <= 0x407b53d1) { /* \|x\| ~<= 5pi/4 /
		59	if (ix <= 0x4016cbe3) { /* \|x\| ~<= 3pi/4 */
		60	if (hx > 0) {
		61	__kernel_sincosdf(x - p1pio2, cs, sn);
		62	cs = -cs;
		63	} else {
		64	__kernel_sincosdf(x + p1pio2, cs, sn);
65	sn = -sn;
66	}
67	} else {
68	if (hx > 0)
69	__kernel_sincosdf(x - p2pio2, sn, cs);
70	else
71	__kernel_sincosdf(x + p2pio2, sn, cs);
72	sn = -sn;
73	cs = -cs;
74	}
75	return;
76	}
77
78	if (ix <= 0x40e231d5) { /* \|x\| ~<= 9pi/4 /
79	if (ix <= 0x40afeddf) { /* \|x\| ~<= 7pi/4 /
80	if (hx > 0) {
81	__kernel_sincosdf(x - p3pio2, cs, sn);
82	sn = -sn;
83	} else {
84	__kernel_sincosdf(x + p3pio2, cs, sn);
85	cs = -cs;
86	}
87	} else {
88	if (hx > 0)
89	__kernel_sincosdf(x - p4pio2, sn, cs);
90	else
91	__kernel_sincosdf(x + p4pio2, sn, cs);
92	}
93	return;
94	}
95
96	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
97	if (ix >= 0x7f800000) {
98	*sn = x - x;
99	*cs = x - x;
100	return;
101	}
102
103	/* Argument reduction. */
104	n = __ieee754_rem_pio2f(x, &y);
105	__kernel_sincosdf(y, &s, &c);
106
107	switch(n & 3) {
108	case 0:
109	*sn = s;
110	*cs = c;
111	break;
112	case 1:
113	*sn = c;
114	*cs = -s;
115	break;
116	case 2:
117	*sn = -s;
118	*cs = -c;
119	break;
120	default:
121	*sn = -c;
122	*cs = s;
123	}
124	}
125
126

Line 0 Link Here

(-)src/s_sincosl.c (+105 lines)
		1	/*-
		2	* Copyright (c) 2007, 2010-2013 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	* s_sinl.c and s_cosl.c merged by Steven G. Kargl.
		27	*/
		28
		29	#include <sys/cdefs.h>
		30	__FBSDID("$FreeBSD$");
		31
		32	#include <float.h>
		33	#ifdef __i386__
		34	#include <ieeefp.h>
		35	#endif
		36
		37	#include "math.h"
		38	#include "math_private.h"
		39	#include "k_sincosl.h"
		40
		41	#if LDBL_MANT_DIG == 64
		42	#include "../ld80/e_rem_pio2l.h"
		43	#elif LDBL_MANT_DIG == 113
		44	#include "../ld128/e_rem_pio2l.h"
		45	#else
		46	#error "Unsupported long double format"
		47	#endif
		48
		49	void
		50	sincosl(long double x, long double sn, long double cs)
		51	{
		52	union IEEEl2bits z;
		53	int e0, sgn;
		54	long double y[2];
		55
		56	z.e = x;
		57	sgn = z.bits.sign;
		58	z.bits.sign = 0;
		59
		60	ENTERI();
		61
		62	/* Optimize the case where x is already within range. */
		63	if (z.e < M_PI_4) {
		64	/*
65	* If x = +-0 or x is a subnormal number, then sin(x) = x and
66	* cos(x) = 1.
67	*/
68	if (z.bits.exp == 0) {
69	*sn = x;
70	*cs = 1;
71	} else
72	__kernel_sincosl(x, 0, 0, sn, cs);
73	RETURNV();
74	}
75
76	/* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
77	if (z.bits.exp == 32767) {
78	*sn = x - x;
79	*cs = x - x;
80	RETURNV();
81	}
82
83	/* Range reduction. */
84	e0 = __ieee754_rem_pio2l(x, y);
85
86	switch (e0 & 3) {
87	case 0:
88	__kernel_sincosl(y[0], y[1], 1, sn, cs);
89	break;
90	case 1:
91	__kernel_sincosl(y[0], y[1], 1, cs, sn);
92	cs = -cs;
93	break;
94	case 2:
95	__kernel_sincosl(y[0], y[1], 1, sn, cs);
96	sn = -sn;
97	cs = -cs;
98	break;
99	default:
100	__kernel_sincosl(y[0], y[1], 1, cs, sn);
101	sn = -sn;
102	}
103
104	RETURNV();
105	}