Attachment 181395 Details for Bug 218300 – patch

[patch] patch

sincos.diff (text/plain), 21.52 KB, created by sgk on 2017-04-01 23:39:07 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: sgk

Created: 2017-04-01 23:39:07 UTC

Size: 21.52 KB

patch

obsolete

>Index: Makefile
>===================================================================
>--- Makefile	(revision 316379)
>+++ Makefile	(working copy)
>@@ -73,7 +73,8 @@
> 	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_sinf.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
> 
>@@ -104,7 +105,8 @@
> 	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_tanhl.c s_tanl.c s_truncl.c w_cabsl.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
>@@ -137,7 +139,8 @@
> 	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 sinh.3 sqrt.3 tan.3 tanh.3 trunc.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
>@@ -215,6 +218,7 @@
> 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
>Index: Symbol.map
>===================================================================
>--- Symbol.map	(revision 316379)
>+++ Symbol.map	(working copy)
>@@ -294,4 +294,7 @@
> 	casinl;
> 	catanl;
> 	catanhl;
>+	sincos;
>+	sincosf;
>+	sincosl;
> };
>Index: man/sincos.3
>===================================================================
>--- man/sincos.3	(nonexistent)
>+++ man/sincos.3	(working copy)
>@@ -0,0 +1,81 @@
>+.\" 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 .
>+
>
>Index: src/k_sincos.h
>===================================================================
>--- src/k_sincos.h	(nonexistent)
>+++ src/k_sincos.h	(working copy)
>@@ -0,0 +1,49 @@
>+/*
>+ * ====================================================
>+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
>+ *
>+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
>+ * Permission to use, copy, modify, and distribute this
>+ * software is freely granted, provided that this notice 
>+ * is preserved.
>+ * ====================================================
>+ * 
>+ * k_sin.c and k_cos.c merged by Steven G. Kargl.
>+ */
>+
>+static const double
>+S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
>+S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
>+S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
>+S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
>+S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
>+S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
>+
>+static const double
>+C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
>+C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
>+C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
>+C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
>+C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
>+C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
>+
>+static inline void
>+__kernel_sincos(double x, double y, int iy, double *sn, double *cs)
>+{
>+	double hz, r, v, w, z;
>+
>+	z = x * x;
>+	w = z * z;
>+	r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
>+	v = z * x;
>+
>+	if (iy == 0)
>+		*sn = x + v * (S1 + z * r);
>+	else
>+		*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
>+
>+	r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
>+	hz = z / 2;
>+	w = 1 - hz;
>+	*cs = w + (((1 - w) - hz) + (z * r - x * y));
>+}
>
>Index: src/k_sincosf.h
>===================================================================
>--- src/k_sincosf.h	(nonexistent)
>+++ src/k_sincosf.h	(working copy)
>@@ -0,0 +1,40 @@
>+/*
>+ * ====================================================
>+ * 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.
>+ * ====================================================
>+ *
>+ * k_sinf.c and k_cosf.c merged by Steven G. Kargl.
>+ */
>+
>+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
>+static const double
>+S1 = -0x15555554cbac77.0p-55,	/* -0.166666666416265235595 */
>+S2 =  0x111110896efbb2.0p-59,	/*  0.0083333293858894631756 */
>+S3 = -0x1a00f9e2cae774.0p-65,	/* -0.000198393348360966317347 */
>+S4 =  0x16cd878c3b46a7.0p-71;	/*  0.0000027183114939898219064 */
>+
>+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
>+static const double
>+C0  = -0x1ffffffd0c5e81.0p-54,	/* -0.499999997251031003120 */
>+C1  =  0x155553e1053a42.0p-57,	/*  0.0416666233237390631894 */
>+C2  = -0x16c087e80f1e27.0p-62,	/* -0.00138867637746099294692 */
>+C3  =  0x199342e0ee5069.0p-68;	/*  0.0000243904487962774090654 */
>+
>+static inline void
>+__kernel_sincosdf(double x, float *sn, float *cs)
>+{
>+	double r, s, w, z;
>+
>+	z = x * x;
>+	w = z * z;
>+	r = S3 + z * S4;
>+	s = z * x;
>+	*sn = (x + s * (S1 + z * S2)) + s * w * r;
>+	r = C2 + z * C3;
>+	*cs = ((1 + z * C0) + w * C1) + (w * z) * r;
>+}
>
>Index: src/k_sincosl.h
>===================================================================
>--- src/k_sincosl.h	(nonexistent)
>+++ src/k_sincosl.h	(working copy)
>@@ -0,0 +1,131 @@
>+/*
>+ * ====================================================
>+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
>+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
>+ *
>+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
>+ * Permission to use, copy, modify, and distribute this
>+ * software is freely granted, provided that this notice 
>+ * is preserved.
>+ * ====================================================
>+ *
>+ * k_sinl.c and k_cosl.c merged by Steven G. Kargl
>+ */
>+
>+#if LDBL_MANT_DIG == 64		/* ld80 version of k_sincosl.c. */
>+
>+#if defined(__amd64__) || defined(__i386__)
>+/* Long double constants are slow on these arches, and broken on i386. */
>+static const volatile double
>+C1hi = 0.041666666666666664,		/*  0x15555555555555.0p-57 */
>+C1lo = 2.2598839032744733e-18,		/*  0x14d80000000000.0p-111 */
>+S1hi = -0.16666666666666666,		/* -0x15555555555555.0p-55 */
>+S1lo = -9.2563760475949941e-18;		/* -0x15580000000000.0p-109 */
>+#define	S1	((long double)S1hi + S1lo)
>+#define	C1	((long double)C1hi + C1lo)
>+#else
>+static const long double
>+C1 =  0.0416666666666666666136L;	/*  0xaaaaaaaaaaaaaa9b.0p-68 */
>+S1 = -0.166666666666666666671L,		/* -0xaaaaaaaaaaaaaaab.0p-66 */
>+#endif
>+
>+static const double
>+C2 = -0.0013888888888888874,		/* -0x16c16c16c16c10.0p-62 */
>+C3 =  0.000024801587301571716,		/*  0x1a01a01a018e22.0p-68 */
>+C4 = -0.00000027557319215507120,	/* -0x127e4fb7602f22.0p-74 */
>+C5 =  0.0000000020876754400407278,	/*  0x11eed8caaeccf1.0p-81 */
>+C6 = -1.1470297442401303e-11,		/* -0x19393412bd1529.0p-89 */
>+C7 =  4.7383039476436467e-14,		/*  0x1aac9d9af5c43e.0p-97 */
>+S2 =  0.0083333333333333332,		/*  0x11111111111111.0p-59 */
>+S3 = -0.00019841269841269427,		/* -0x1a01a01a019f81.0p-65 */
>+S4 =  0.0000027557319223597490,		/*  0x171de3a55560f7.0p-71 */
>+S5 = -0.000000025052108218074604,	/* -0x1ae64564f16cad.0p-78 */
>+S6 =  1.6059006598854211e-10,		/*  0x161242b90243b5.0p-85 */
>+S7 = -7.6429779983024564e-13,		/* -0x1ae42ebd1b2e00.0p-93 */
>+S8 =  2.6174587166648325e-15;		/*  0x179372ea0b3f64.0p-101 */
>+
>+static inline void
>+__kernel_sincosl(long double x, long double y, int iy, long double *sn,
>+    long double *cs)
>+{
>+	long double hz, r, v, w, z;
>+
>+	z = x * x;
>+	v = z * x;
>+	/*
>+	 * XXX Replace Horner scheme with an algorithm suitable for CPUs
>+	 * with more complex pipelines.
>+	 */
>+	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8)))));
>+
>+	if (iy == 0)
>+		*sn = x + v * (S1 + z * r);
>+	else
>+		*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
>+
>+	hz = z / 2;
>+	w = 1 - hz;
>+	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
>+	    z * C7))))));
>+	*cs = w + (((1 - w) - hz) + (z * r - x * y));
>+}
>+
>+#elif LDBL_MANT_DIG == 113	/* ld128 version of k_sincosl.c. */
>+
>+static const long double
>+C1 =  0.04166666666666666666666666666666658424671L,
>+C2 = -0.001388888888888888888888888888863490893732L,
>+C3 =  0.00002480158730158730158730158600795304914210L,
>+C4 = -0.2755731922398589065255474947078934284324e-6L,
>+C5 =  0.2087675698786809897659225313136400793948e-8L,
>+C6 = -0.1147074559772972315817149986812031204775e-10L,
>+C7 =  0.4779477332386808976875457937252120293400e-13L,
>+S1 = -0.16666666666666666666666666666666666606732416116558L,
>+S2 =  0.0083333333333333333333333333333331135404851288270047L,
>+S3 = -0.00019841269841269841269841269839935785325638310428717L,
>+S4 =  0.27557319223985890652557316053039946268333231205686e-5L,
>+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
>+S6 =  0.16059043836821614596571832194524392581082444805729e-9L,
>+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
>+S8 =  0.28114572543451292625024967174638477283187397621303e-14L;
>+
>+static const double
>+C8  = -0.1561920696721507929516718307820958119868e-15,
>+C9  =  0.4110317413744594971475941557607804508039e-18,
>+C10 = -0.8896592467191938803288521958313920156409e-21,
>+C11 =  0.1601061435794535138244346256065192782581e-23,
>+S9  = -0.82206352458348947812512122163446202498005154296863e-17,
>+S10 =  0.19572940011906109418080609928334380560135358385256e-19,
>+S11 = -0.38680813379701966970673724299207480965452616911420e-22,
>+S12 =  0.64038150078671872796678569586315881020659912139412e-25;
>+
>+static inline void
>+__kernel_sincosl(long double x, long double y, int iy, long double *sn, 
>+    long double *cs)
>+{
>+	long double hz, r, v, w, z;
>+
>+	z = x * x;
>+	v = z * x;
>+	/*
>+	 * XXX Replace Horner scheme with an algorithm suitable for CPUs
>+	 * with more complex pipelines.
>+	 */
>+	r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
>+	    z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
>+
>+	if (iy == 0)
>+		*sn = x + v * (S1 + z * r);
>+	else
>+		*cs = x - ((z * (y / 2 - v * r) - y) - v * S1);
>+
>+	hz = z / 2;
>+	w = 1 - hz;
>+	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + 
>+	    z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
>+
>+	*cs =  w + (((1 - w) - hz) + (z * r - x * y));
>+}
>+#else
>+#error "Unsupported long double format"
>+#endif
>
>Index: src/math.h
>===================================================================
>--- src/math.h	(revision 316379)
>+++ src/math.h	(working copy)
>@@ -500,6 +500,9 @@
> 
> #if __BSD_VISIBLE
> long double	lgammal_r(long double, int *);
>+void		sincos(double, double *, double *);
>+void		sincosf(float, float *, float *);
>+void		sincosl(long double, long double *, long double *);
> #endif
> 
> __END_DECLS
>Index: src/math_private.h
>===================================================================
>--- src/math_private.h	(revision 316379)
>+++ src/math_private.h	(working copy)
>@@ -306,9 +306,15 @@
> 		fpsetprec(__oprec);		\
> 	RETURNF(__retval);			\
> } while (0)
>+#define	RETURNV() do {				\
>+	if (__oprec != FP_PE)			\
>+		fpsetprec(__oprec);		\
>+	return;			\
>+} while (0)
> #else
> #define	ENTERI(x)
> #define	RETURNI(x)	RETURNF(x)
>+#define	RETURNV()	return
> #endif
> 
> /* Default return statement if hack*_t() is not used. */
>Index: src/s_sincos.c
>===================================================================
>--- src/s_sincos.c	(nonexistent)
>+++ src/s_sincos.c	(working copy)
>@@ -0,0 +1,80 @@
>+/*
>+ * ====================================================
>+ * 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
>
>Index: src/s_sincosf.c
>===================================================================
>--- src/s_sincosf.c	(nonexistent)
>+++ src/s_sincosf.c	(working copy)
>@@ -0,0 +1,126 @@
>+/* 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;
>+	}
>+}
>+
>+
>
>Index: src/s_sincosl.c
>===================================================================
>--- src/s_sincosl.c	(nonexistent)
>+++ src/s_sincosl.c	(working copy)
>@@ -0,0 +1,105 @@
>+/*-
>+ * 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();
>+}

Actions: View | Diff

Attachments on bug 218300: 181395 | 182677