Module Name: src
Committed By: christos
Date: Sat Aug 27 08:31:59 UTC 2022
Modified Files:
src/distrib/sets/lists/comp: mi
src/distrib/sets/lists/debug: mi
src/distrib/sets/lists/tests: mi
src/include: math.h
src/lib/libm: Makefile
src/lib/libm/noieee_src: n_sincos.c
src/lib/libm/src: math_private.h namespace.h
src/tests/lib/libm: Makefile
Added Files:
src/lib/libm/ld128: e_rem_pio2l.h
src/lib/libm/ld80: e_rem_pio2l.h
src/lib/libm/man: sincos.3
src/lib/libm/src: e_rem_pio2f.h e_rem_pio2l.h k_sincos.h k_sincosf.h
k_sincosl.h s_sincos.c s_sincosf.c s_sincosl.c
src/tests/lib/libm: t_sincos.c
Log Message:
Add sincos{,f,l} from FreeBSD
To generate a diff of this commit:
cvs rdiff -u -r1.2417 -r1.2418 src/distrib/sets/lists/comp/mi
cvs rdiff -u -r1.387 -r1.388 src/distrib/sets/lists/debug/mi
cvs rdiff -u -r1.1219 -r1.1220 src/distrib/sets/lists/tests/mi
cvs rdiff -u -r1.66 -r1.67 src/include/math.h
cvs rdiff -u -r1.216 -r1.217 src/lib/libm/Makefile
cvs rdiff -u -r0 -r1.1 src/lib/libm/ld128/e_rem_pio2l.h
cvs rdiff -u -r0 -r1.1 src/lib/libm/ld80/e_rem_pio2l.h
cvs rdiff -u -r0 -r1.1 src/lib/libm/man/sincos.3
cvs rdiff -u -r1.7 -r1.8 src/lib/libm/noieee_src/n_sincos.c
cvs rdiff -u -r0 -r1.1 src/lib/libm/src/e_rem_pio2f.h \
src/lib/libm/src/e_rem_pio2l.h src/lib/libm/src/k_sincos.h \
src/lib/libm/src/k_sincosf.h src/lib/libm/src/k_sincosl.h \
src/lib/libm/src/s_sincos.c src/lib/libm/src/s_sincosf.c \
src/lib/libm/src/s_sincosl.c
cvs rdiff -u -r1.25 -r1.26 src/lib/libm/src/math_private.h
cvs rdiff -u -r1.15 -r1.16 src/lib/libm/src/namespace.h
cvs rdiff -u -r1.47 -r1.48 src/tests/lib/libm/Makefile
cvs rdiff -u -r0 -r1.1 src/tests/lib/libm/t_sincos.c
Please note that diffs are not public domain; they are subject to the
copyright notices on the relevant files.
Modified files:
Index: src/distrib/sets/lists/comp/mi
diff -u src/distrib/sets/lists/comp/mi:1.2417 src/distrib/sets/lists/comp/mi:1.2418
--- src/distrib/sets/lists/comp/mi:1.2417 Fri Jul 22 11:43:36 2022
+++ src/distrib/sets/lists/comp/mi Sat Aug 27 04:31:58 2022
@@ -1,4 +1,4 @@
-# $NetBSD: mi,v 1.2417 2022/07/22 15:43:36 wiz Exp $
+# $NetBSD: mi,v 1.2418 2022/08/27 08:31:58 christos Exp $
#
# Note: don't delete entries from here - mark them as "obsolete" instead.
./etc/mtree/set.comp comp-sys-root
@@ -10280,6 +10280,7 @@
./usr/share/man/cat3/simpleq_next.0 comp-obsolete obsolete
./usr/share/man/cat3/simpleq_remove_head.0 comp-obsolete obsolete
./usr/share/man/cat3/sin.0 comp-c-catman .cat
+./usr/share/man/cat3/sincos.0 comp-c-catman .cat
./usr/share/man/cat3/sinf.0 comp-c-catman .cat
./usr/share/man/cat3/sinh.0 comp-c-catman .cat
./usr/share/man/cat3/sinhf.0 comp-c-catman .cat
@@ -18562,6 +18563,7 @@
./usr/share/man/html3/sigsetops.html comp-c-htmlman html
./usr/share/man/html3/sigvec.html comp-c-htmlman html
./usr/share/man/html3/sin.html comp-c-htmlman html
+./usr/share/man/html3/sincos.html comp-c-htmlman html
./usr/share/man/html3/sinf.html comp-c-htmlman html
./usr/share/man/html3/sinh.html comp-c-htmlman html
./usr/share/man/html3/sinhf.html comp-c-htmlman html
@@ -26878,6 +26880,7 @@
./usr/share/man/man3/simpleq_next.3 comp-obsolete obsolete
./usr/share/man/man3/simpleq_remove_head.3 comp-obsolete obsolete
./usr/share/man/man3/sin.3 comp-c-man .man
+./usr/share/man/man3/sincos.3 comp-c-man .man
./usr/share/man/man3/sinf.3 comp-c-man .man
./usr/share/man/man3/sinh.3 comp-c-man .man
./usr/share/man/man3/sinhf.3 comp-c-man .man
Index: src/distrib/sets/lists/debug/mi
diff -u src/distrib/sets/lists/debug/mi:1.387 src/distrib/sets/lists/debug/mi:1.388
--- src/distrib/sets/lists/debug/mi:1.387 Mon Jun 6 06:56:27 2022
+++ src/distrib/sets/lists/debug/mi Sat Aug 27 04:31:58 2022
@@ -1,4 +1,4 @@
-# $NetBSD: mi,v 1.387 2022/06/06 10:56:27 nia Exp $
+# $NetBSD: mi,v 1.388 2022/08/27 08:31:58 christos Exp $
./etc/mtree/set.debug comp-sys-root
./usr/lib comp-sys-usr compatdir
./usr/lib/i18n/libBIG5_g.a comp-c-debuglib debuglib,compatfile
@@ -2311,6 +2311,7 @@
./usr/libdata/debug/usr/tests/lib/libm/t_round.debug tests-lib-debug debug,atf,compattestfile
./usr/libdata/debug/usr/tests/lib/libm/t_scalbn.debug tests-lib-debug debug,atf,compattestfile
./usr/libdata/debug/usr/tests/lib/libm/t_sin.debug tests-lib-debug debug,atf,compattestfile
+./usr/libdata/debug/usr/tests/lib/libm/t_sincos.debug tests-lib-debug debug,atf,compattestfile
./usr/libdata/debug/usr/tests/lib/libm/t_sinh.debug tests-lib-debug debug,atf,compattestfile
./usr/libdata/debug/usr/tests/lib/libm/t_sqrt.debug tests-lib-debug debug,atf,compattestfile
./usr/libdata/debug/usr/tests/lib/libm/t_tan.debug tests-lib-debug debug,atf,compattestfile
Index: src/distrib/sets/lists/tests/mi
diff -u src/distrib/sets/lists/tests/mi:1.1219 src/distrib/sets/lists/tests/mi:1.1220
--- src/distrib/sets/lists/tests/mi:1.1219 Fri Aug 12 06:49:17 2022
+++ src/distrib/sets/lists/tests/mi Sat Aug 27 04:31:58 2022
@@ -1,4 +1,4 @@
-# $NetBSD: mi,v 1.1219 2022/08/12 10:49:17 riastradh Exp $
+# $NetBSD: mi,v 1.1220 2022/08/27 08:31:58 christos Exp $
#
# Note: don't delete entries from here - mark them as "obsolete" instead.
#
@@ -3841,6 +3841,7 @@
./usr/tests/lib/libm/t_round tests-lib-tests compattestfile,atf
./usr/tests/lib/libm/t_scalbn tests-lib-tests compattestfile,atf
./usr/tests/lib/libm/t_sin tests-lib-tests compattestfile,atf
+./usr/tests/lib/libm/t_sincos tests-lib-tests compattestfile,atf
./usr/tests/lib/libm/t_sinh tests-lib-tests compattestfile,atf
./usr/tests/lib/libm/t_sqrt tests-lib-tests compattestfile,atf
./usr/tests/lib/libm/t_tan tests-lib-tests compattestfile,atf
Index: src/include/math.h
diff -u src/include/math.h:1.66 src/include/math.h:1.67
--- src/include/math.h:1.66 Sat Feb 22 17:47:35 2020
+++ src/include/math.h Sat Aug 27 04:31:59 2022
@@ -1,4 +1,4 @@
-/* $NetBSD: math.h,v 1.66 2020/02/22 22:47:35 joerg Exp $ */
+/* $NetBSD: math.h,v 1.67 2022/08/27 08:31:59 christos Exp $ */
/*
* ====================================================
@@ -553,6 +553,10 @@ float significandf(float);
* float versions of BSD math library entry points
*/
float dremf(float, float);
+
+void sincos(double, double *, double *);
+void sincosf(float, float *, float *);
+void sincosl(long double, long double *, long double *);
#endif /* _NETBSD_SOURCE */
#if defined(_NETBSD_SOURCE) || defined(_REENTRANT)
Index: src/lib/libm/Makefile
diff -u src/lib/libm/Makefile:1.216 src/lib/libm/Makefile:1.217
--- src/lib/libm/Makefile:1.216 Thu Jun 23 12:42:50 2022
+++ src/lib/libm/Makefile Sat Aug 27 04:31:58 2022
@@ -1,4 +1,4 @@
-# $NetBSD: Makefile,v 1.216 2022/06/23 16:42:50 martin Exp $
+# $NetBSD: Makefile,v 1.217 2022/08/27 08:31:58 christos Exp $
#
# @(#)Makefile 5.1beta 93/09/24
#
@@ -282,7 +282,8 @@ COMMON_SRCS+= b_exp.c b_log.c b_tgamma.c
s_matherr.c s_modff.c s_modfl.c s_nearbyint.c s_nextafter.c s_nextafterl.c \
s_nextafterf.c s_nexttowardf.c s_remquo.c s_remquof.c s_rint.c s_rintf.c \
s_round.c s_roundf.c s_roundl.c s_scalbn.c \
- s_scalbnf.c s_scalbnl.c s_signgam.c s_significand.c s_significandf.c s_sin.c \
+ s_scalbnf.c s_scalbnl.c s_signgam.c s_significand.c s_significandf.c \
+ s_sincos.c s_sincosf.c s_sincosl.c s_sin.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 s_truncl.c \
w_acos.c w_acosf.c w_acosh.c w_acoshf.c w_asin.c w_asinf.c w_atan2.c \
@@ -354,7 +355,7 @@ MAN+= acos.3 acosh.3 asin.3 asinh.3 atan
ieee_test.3 ilogb.3 isinff.3 j0.3 ldexp.3 lgamma.3 log.3 lrint.3 \
math.3 modf.3 nextafter.3 pow.3 \
remainder.3 rint.3 round.3 \
- scalbn.3 sin.3 sinh.3 sqrt.3 \
+ scalbn.3 sincos.3 sin.3 sinh.3 sqrt.3 \
tan.3 tanh.3 trunc.3 fmax.3 fdim.3
# fenv.h interface
Index: src/lib/libm/noieee_src/n_sincos.c
diff -u src/lib/libm/noieee_src/n_sincos.c:1.7 src/lib/libm/noieee_src/n_sincos.c:1.8
--- src/lib/libm/noieee_src/n_sincos.c:1.7 Fri Oct 10 16:58:09 2014
+++ src/lib/libm/noieee_src/n_sincos.c Sat Aug 27 04:31:59 2022
@@ -1,4 +1,4 @@
-/* $NetBSD: n_sincos.c,v 1.7 2014/10/10 20:58:09 martin Exp $ */
+/* $NetBSD: n_sincos.c,v 1.8 2022/08/27 08:31:59 christos Exp $ */
/*
* Copyright (c) 1987, 1993
* The Regents of the University of California. All rights reserved.
@@ -41,6 +41,7 @@ static char sccsid[] = "@(#)sincos.c 8.1
#ifdef __weak_alias
__weak_alias(_sinl, sin);
__weak_alias(_cosl, cos);
+__weak_alias(_sincosl, sincos);
#endif
double
@@ -113,3 +114,17 @@ cosf(float x)
{
return cos(x);
}
+
+void
+sincos(double x, double *s, double *c)
+{
+ *s = sin(x);
+ *c = cos(x);
+}
+
+void
+sincosf(float x, float *s, float *c)
+{
+ *s = sinf(x);
+ *c = cosf(x);
+}
Index: src/lib/libm/src/math_private.h
diff -u src/lib/libm/src/math_private.h:1.25 src/lib/libm/src/math_private.h:1.26
--- src/lib/libm/src/math_private.h:1.25 Wed Aug 24 09:51:19 2022
+++ src/lib/libm/src/math_private.h Sat Aug 27 04:31:59 2022
@@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
- * $NetBSD: math_private.h,v 1.25 2022/08/24 13:51:19 christos Exp $
+ * $NetBSD: math_private.h,v 1.26 2022/08/27 08:31:59 christos Exp $
*/
#ifndef _MATH_PRIVATE_H_
@@ -194,6 +194,39 @@ do { \
} while (0)
#endif
+/* Support switching the mode to FP_PE if necessary. */
+#if defined(__i386__) && !defined(NO_FPSETPREC)
+#define ENTERI() ENTERIT(long double)
+#define ENTERIT(returntype) \
+ returntype __retval; \
+ fp_prec_t __oprec; \
+ \
+ if ((__oprec = fpgetprec()) != FP_PE) \
+ fpsetprec(FP_PE)
+#define RETURNI(x) do { \
+ __retval = (x); \
+ if (__oprec != FP_PE) \
+ fpsetprec(__oprec); \
+ RETURNF(__retval); \
+} while (0)
+#define ENTERV() \
+ fp_prec_t __oprec; \
+ \
+ if ((__oprec = fpgetprec()) != FP_PE) \
+ fpsetprec(FP_PE)
+#define RETURNV() do { \
+ if (__oprec != FP_PE) \
+ fpsetprec(__oprec); \
+ return; \
+} while (0)
+#else
+#define ENTERI()
+#define ENTERIT(x)
+#define RETURNI(x) RETURNF(x)
+#define ENTERV()
+#define RETURNV() return
+#endif
+
#ifdef _COMPLEX_H
/*
Index: src/lib/libm/src/namespace.h
diff -u src/lib/libm/src/namespace.h:1.15 src/lib/libm/src/namespace.h:1.16
--- src/lib/libm/src/namespace.h:1.15 Sat Oct 26 13:57:20 2019
+++ src/lib/libm/src/namespace.h Sat Aug 27 04:31:59 2022
@@ -1,4 +1,4 @@
-/* $NetBSD: namespace.h,v 1.15 2019/10/26 17:57:20 christos Exp $ */
+/* $NetBSD: namespace.h,v 1.16 2022/08/27 08:31:59 christos Exp $ */
#define atan2 _atan2
#define atan2f _atan2f
@@ -22,6 +22,9 @@
#define finite _finite
#define finitef _finitef
#endif /* notyet */
+#define sincos _sincos
+#define sincosf _sincosf
+#define sincosl _sincosl
#define sinh _sinh
#define sinhf _sinhf
#define sinhl _sinhl
Index: src/tests/lib/libm/Makefile
diff -u src/tests/lib/libm/Makefile:1.47 src/tests/lib/libm/Makefile:1.48
--- src/tests/lib/libm/Makefile:1.47 Sun Jun 21 02:58:16 2020
+++ src/tests/lib/libm/Makefile Sat Aug 27 04:31:58 2022
@@ -1,4 +1,4 @@
-# $NetBSD: Makefile,v 1.47 2020/06/21 06:58:16 lukem Exp $
+# $NetBSD: Makefile,v 1.48 2022/08/27 08:31:58 christos Exp $
.include <bsd.own.mk>
@@ -40,6 +40,7 @@ TESTS_C+= t_precision
TESTS_C+= t_round
TESTS_C+= t_scalbn
TESTS_C+= t_sin
+TESTS_C+= t_sincos
TESTS_C+= t_sinh
TESTS_C+= t_sqrt
TESTS_C+= t_tan
Added files:
Index: src/lib/libm/ld128/e_rem_pio2l.h
diff -u /dev/null src/lib/libm/ld128/e_rem_pio2l.h:1.1
--- /dev/null Sat Aug 27 04:31:59 2022
+++ src/lib/libm/ld128/e_rem_pio2l.h Sat Aug 27 04:31:58 2022
@@ -0,0 +1,139 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/ld128/e_rem_pio2l.h 336545 2018-07-20 12:42:24Z bde $");
+#endif
+
+/* ld128 version of __ieee754_rem_pio2l(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+#include <machine/ieee.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * XXX need to verify that nonzero integer multiples of pi/2 within the
+ * range get no closer to a long double than 2**-140, or that
+ * ilogb(x) + ilogb(min_delta) < 45 - -140.
+ */
+/*
+ * invpio2: 113 bits of 2/pi
+ * pio2_1: first 68 bits of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ * pio2_2: second 68 bits of pi/2
+ * pio2_2t: pi/2 - (pio2_1+pio2_2)
+ * pio2_3: third 68 bits of pi/2
+ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 = 1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */
+
+static const long double
+invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
+pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */
+pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */
+pio2_2 = 2.0222662487959507323994779168837751e-21L, /* 0x13198a2e03707344a400000000000.0p-181 */
+pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df531d89cd91.0p-254 */
+pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */
+pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
+
+static inline __always_inline int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+ union ieee_ext_u u,u1;
+ long double z,w,t,r,fn;
+ double tx[5],ty[3];
+ int64_t n;
+ int e0,ex,i,j,nx;
+ int16_t expsign, expsign1;
+
+ u.extu_ld = x;
+ ex = u.extu_exp;
+ expsign = u.extu_exp | (u.extu_sign << 15);
+ if (ex < BIAS + 45 || (ex == BIAS + 45 &&
+ u.extu_frach < 0x921fb54442d1LL)) {
+ /* |x| ~< 2^45*(pi/2), medium size */
+ /* TODO: use only double precision for fn, as in expl(). */
+ fn = rnintl(x * invpio2);
+ n = i64rint(fn);
+ r = x-fn*pio2_1;
+ w = fn*pio2_1t; /* 1st round good to 180 bit */
+ {
+ union ieee_ext_u u2;
+ int ex1;
+ j = ex;
+ y[0] = r-w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_exp;
+ i = j-ex1;
+ if(i>51) { /* 2nd iteration needed, good to 248 */
+ t = r;
+ w = fn*pio2_2;
+ r = t-w;
+ w = fn*pio2_2t-((t-r)-w);
+ y[0] = r-w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_exp;
+ i = j-ex1;
+ if(i>119) { /* 3rd iteration need, 316 bits acc */
+ t = r; /* will cover all possible cases */
+ w = fn*pio2_3;
+ r = t-w;
+ w = fn*pio2_3t-((t-r)-w);
+ y[0] = r-w;
+ }
+ }
+ }
+ y[1] = (r-y[0])-w;
+ return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(ex==0x7fff) { /* x is inf or NaN */
+ y[0]=y[1]=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(x)-23) */
+ u1.extu_ld = x;
+ e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
+ expsign1 = ex - e0;
+ u1.extu_exp = expsign1;
+ u1.extu_sign = expsign1 >> 15;
+ z = u1.extu_ld;
+ for(i=0;i<4;i++) {
+ tx[i] = (double)((int32_t)(z));
+ z = (z-tx[i])*two24;
+ }
+ tx[4] = z;
+ nx = 5;
+ while(tx[nx-1]==zero) nx--; /* skip zero term */
+ n = __kernel_rem_pio2(tx,ty,e0,nx,3);
+ t = (long double)ty[2] + ty[1];
+ r = t + ty[0];
+ w = ty[0] - (r - t);
+ if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+ y[0] = r; y[1] = w; return n;
+}
Index: src/lib/libm/ld80/e_rem_pio2l.h
diff -u /dev/null src/lib/libm/ld80/e_rem_pio2l.h:1.1
--- /dev/null Sat Aug 27 04:31:59 2022
+++ src/lib/libm/ld80/e_rem_pio2l.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,147 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/ld80/e_rem_pio2l.h 336545 2018-07-20 12:42:24Z bde $");
+#endif
+
+/* ld80 version of __ieee754_rem_pio2l(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+#include <machine/ieee.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * invpio2: 64 bits of 2/pi
+ * pio2_1: first 39 bits of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ * pio2_2: second 39 bits of pi/2
+ * pio2_2t: pi/2 - (pio2_1+pio2_2)
+ * pio2_3: third 39 bits of pi/2
+ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+pio2_1 = 1.57079632679597125389e+00, /* 0x3FF921FB, 0x54444000 */
+pio2_2 = -1.07463465549783099519e-12, /* -0x12e7b967674000.0p-92 */
+pio2_3 = 6.36831716351370313614e-25; /* 0x18a2e037074000.0p-133 */
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+invpio2hi = 6.3661977236758138e-01, /* 0x145f306dc9c883.0p-53 */
+invpio2lo = -3.9356538861223811e-17, /* -0x16b00000000000.0p-107 */
+pio2_1thi = -1.0746346554971943e-12, /* -0x12e7b9676733af.0p-92 */
+pio2_1tlo = 8.8451028997905949e-29, /* 0x1c080000000000.0p-146 */
+pio2_2thi = 6.3683171635109499e-25, /* 0x18a2e03707344a.0p-133 */
+pio2_2tlo = 2.3183081793789774e-41, /* 0x10280000000000.0p-187 */
+pio2_3thi = -2.7529965190440717e-37, /* -0x176b7ed8fbbacc.0p-174 */
+pio2_3tlo = -4.2006647512740502e-54; /* -0x19c00000000000.0p-230 */
+#define invpio2 ((long double)invpio2hi + invpio2lo)
+#define pio2_1t ((long double)pio2_1thi + pio2_1tlo)
+#define pio2_2t ((long double)pio2_2thi + pio2_2tlo)
+#define pio2_3t ((long double)pio2_3thi + pio2_3tlo)
+#else
+static const long double
+invpio2 = 6.36619772367581343076e-01L, /* 0xa2f9836e4e44152a.0p-64 */
+pio2_1t = -1.07463465549719416346e-12L, /* -0x973dcb3b399d747f.0p-103 */
+pio2_2t = 6.36831716351095013979e-25L, /* 0xc51701b839a25205.0p-144 */
+pio2_3t = -2.75299651904407171810e-37L; /* -0xbb5bf6c7ddd660ce.0p-185 */
+#endif
+
+static inline __always_inline int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+ union ieee_ext_u u, u1;
+ long double z,w,t,r,fn;
+ double tx[3],ty[2];
+ int e0,ex,i,j,nx,n;
+ int16_t expsign, expsign1;
+
+ u.extu_ld = x;
+ ex = u.extu_exp;
+ expsign = u.extu_exp | (u.extu_sign << EXT_EXPBITS);
+ if (ex < BIAS + 25 || (ex == BIAS + 25 && u.extu_frach < 0xc90fdaa2U)) {
+ /* |x| ~< 2^25*(pi/2), medium size */
+ fn = rnintl(x*invpio2);
+ n = irint(fn);
+ r = x-fn*pio2_1;
+ w = fn*pio2_1t; /* 1st round good to 102 bit */
+ {
+ union ieee_ext_u u2;
+ int ex1;
+ j = ex;
+ y[0] = r-w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_exp;
+ i = j-ex1;
+ if(i>22) { /* 2nd iteration needed, good to 141 */
+ t = r;
+ w = fn*pio2_2;
+ r = t-w;
+ w = fn*pio2_2t-((t-r)-w);
+ y[0] = r-w;
+ u2.extu_ld = y[0];
+ ex1 = u2.extu_exp;
+ i = j-ex1;
+ if(i>61) { /* 3rd iteration need, 180 bits acc */
+ t = r; /* will cover all possible cases */
+ w = fn*pio2_3;
+ r = t-w;
+ w = fn*pio2_3t-((t-r)-w);
+ y[0] = r-w;
+ }
+ }
+ }
+ y[1] = (r-y[0])-w;
+ return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(ex==0x7fff) { /* x is inf or NaN */
+ y[0]=y[1]=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(x)-23) */
+ u1.extu_ld = x;
+ e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
+ expsign1 = ex - e0;
+ u1.extu_exp = expsign1;
+ u1.extu_sign = expsign1 >> EXT_EXPBITS;
+ z = u1.extu_ld;
+ for(i=0;i<2;i++) {
+ tx[i] = (double)((int32_t)(z));
+ z = (z-tx[i])*two24;
+ }
+ tx[2] = z;
+ nx = 3;
+ while(tx[nx-1]==zero) nx--; /* skip zero term */
+ n = __kernel_rem_pio2(tx,ty,e0,nx,2);
+ r = (long double)ty[0] + ty[1];
+ w = ty[1] - (r - ty[0]);
+ if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+ y[0] = r; y[1] = w; return n;
+}
Index: src/lib/libm/man/sincos.3
diff -u /dev/null src/lib/libm/man/sincos.3:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/man/sincos.3 Sat Aug 27 04:31:59 2022
@@ -0,0 +1,85 @@
+.\" 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.
+.\"
+.\" $FreeBSD: head/lib/msun/man/sincos.3 366583 2020-10-09 19:12:44Z gbe $
+.\" $NetBSD: sincos.3,v 1.1 2022/08/27 08:31:59 christos Exp $
+.\"
+.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
+and
+.Nx 10.0
+to aid in writing various complex function contained in
+.St -isoC-99 .
+
Index: src/lib/libm/src/e_rem_pio2f.h
diff -u /dev/null src/lib/libm/src/e_rem_pio2f.h:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/e_rem_pio2f.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,80 @@
+/* e_rem_pio2f.c -- float version of e_rem_pio2.c
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
+ * Debugged and optimized by Bruce D. Evans.
+ */
+
+/*
+ * ====================================================
+ * 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>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/e_rem_pio2f.c 336545 2018-07-20 12:42:24Z bde $");
+#endif
+
+/* __ieee754_rem_pio2fd(x,y)
+ *
+ * return the remainder of x rem pi/2 in *y
+ * use double precision for everything except passing x
+ * use __kernel_rem_pio2() for large x
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+/*
+ * invpio2: 53 bits of 2/pi
+ * pio2_1: first 25 bits of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ */
+
+static const double
+invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
+pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
+pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
+
+#ifdef INLINE_REM_PIO2F
+static __inline __always_inline
+#endif
+int
+__ieee754_rem_pio2fd(float x, double *y)
+{
+ double w,r,fn;
+ double tx[1],ty[1];
+ float z;
+ int32_t e0,n,ix,hx;
+
+ GET_FLOAT_WORD(hx,x);
+ ix = hx&0x7fffffff;
+ /* 33+53 bit pi is good enough for medium size */
+ if(ix<0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */
+ fn = rnint((float_t)x*invpio2);
+ n = irint(fn);
+ r = x-fn*pio2_1;
+ w = fn*pio2_1t;
+ *y = r-w;
+ return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(ix>=0x7f800000) { /* x is inf or NaN */
+ *y=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(|x|)-23) */
+ e0 = (ix>>23)-150; /* e0 = ilogb(|x|)-23; */
+ SET_FLOAT_WORD(z, ix - ((int32_t)(e0<<23)));
+ tx[0] = z;
+ n = __kernel_rem_pio2(tx,ty,e0,1,0);
+ if(hx<0) {*y = -ty[0]; return -n;}
+ *y = ty[0]; return n;
+}
Index: src/lib/libm/src/e_rem_pio2l.h
diff -u /dev/null src/lib/libm/src/e_rem_pio2l.h:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/e_rem_pio2l.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,135 @@
+/* From: @(#)e_rem_pio2.c 1.4 95/01/18 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD: head/lib/msun/ld128/e_rem_pio2l.h 336545 2018-07-20 12:42:24Z bde $");
+
+/* ld128 version of __ieee754_rem_pio2l(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+
+#define BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * XXX need to verify that nonzero integer multiples of pi/2 within the
+ * range get no closer to a long double than 2**-140, or that
+ * ilogb(x) + ilogb(min_delta) < 45 - -140.
+ */
+/*
+ * invpio2: 113 bits of 2/pi
+ * pio2_1: first 68 bits of pi/2
+ * pio2_1t: pi/2 - pio2_1
+ * pio2_2: second 68 bits of pi/2
+ * pio2_2t: pi/2 - (pio2_1+pio2_2)
+ * pio2_3: third 68 bits of pi/2
+ * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 = 1.67772160000000000000e+07; /* 0x41700000, 0x00000000 */
+
+static const long double
+invpio2 = 6.3661977236758134307553505349005747e-01L, /* 0x145f306dc9c882a53f84eafa3ea6a.0p-113 */
+pio2_1 = 1.5707963267948966192292994253909555e+00L, /* 0x1921fb54442d18469800000000000.0p-112 */
+pio2_1t = 2.0222662487959507323996846200947577e-21L, /* 0x13198a2e03707344a4093822299f3.0p-181 */
+pio2_2 = 2.0222662487959507323994779168837751e-21L, /* 0x13198a2e03707344a400000000000.0p-181 */
+pio2_2t = 2.0670321098263988236496903051604844e-43L, /* 0x127044533e63a0105df531d89cd91.0p-254 */
+pio2_3 = 2.0670321098263988236499468110329591e-43L, /* 0x127044533e63a0105e00000000000.0p-254 */
+pio2_3t = -2.5650587247459238361625433492959285e-65L; /* -0x159c4ec64ddaeb5f78671cbfb2210.0p-327 */
+
+static inline __always_inline int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+ union IEEEl2bits u,u1;
+ long double z,w,t,r,fn;
+ double tx[5],ty[3];
+ int64_t n;
+ int e0,ex,i,j,nx;
+ int16_t expsign;
+
+ u.e = x;
+ expsign = u.xbits.expsign;
+ ex = expsign & 0x7fff;
+ if (ex < BIAS + 45 || ex == BIAS + 45 &&
+ u.bits.manh < 0x921fb54442d1LL) {
+ /* |x| ~< 2^45*(pi/2), medium size */
+ /* TODO: use only double precision for fn, as in expl(). */
+ fn = rnintl(x * invpio2);
+ n = i64rint(fn);
+ r = x-fn*pio2_1;
+ w = fn*pio2_1t; /* 1st round good to 180 bit */
+ {
+ union IEEEl2bits u2;
+ int ex1;
+ j = ex;
+ y[0] = r-w;
+ u2.e = y[0];
+ ex1 = u2.xbits.expsign & 0x7fff;
+ i = j-ex1;
+ if(i>51) { /* 2nd iteration needed, good to 248 */
+ t = r;
+ w = fn*pio2_2;
+ r = t-w;
+ w = fn*pio2_2t-((t-r)-w);
+ y[0] = r-w;
+ u2.e = y[0];
+ ex1 = u2.xbits.expsign & 0x7fff;
+ i = j-ex1;
+ if(i>119) { /* 3rd iteration need, 316 bits acc */
+ t = r; /* will cover all possible cases */
+ w = fn*pio2_3;
+ r = t-w;
+ w = fn*pio2_3t-((t-r)-w);
+ y[0] = r-w;
+ }
+ }
+ }
+ y[1] = (r-y[0])-w;
+ return n;
+ }
+ /*
+ * all other (large) arguments
+ */
+ if(ex==0x7fff) { /* x is inf or NaN */
+ y[0]=y[1]=x-x; return 0;
+ }
+ /* set z = scalbn(|x|,ilogb(x)-23) */
+ u1.e = x;
+ e0 = ex - BIAS - 23; /* e0 = ilogb(|x|)-23; */
+ u1.xbits.expsign = ex - e0;
+ z = u1.e;
+ for(i=0;i<4;i++) {
+ tx[i] = (double)((int32_t)(z));
+ z = (z-tx[i])*two24;
+ }
+ tx[4] = z;
+ nx = 5;
+ while(tx[nx-1]==zero) nx--; /* skip zero term */
+ n = __kernel_rem_pio2(tx,ty,e0,nx,3);
+ t = (long double)ty[2] + ty[1];
+ r = t + ty[0];
+ w = ty[0] - (r - t);
+ if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+ y[0] = r; y[1] = w; return n;
+}
Index: src/lib/libm/src/k_sincos.h
diff -u /dev/null src/lib/libm/src/k_sincos.h:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/k_sincos.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,57 @@
+/*-
+ * ====================================================
+ * 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.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/k_sincos.h 319047 2017-05-28 06:13:38Z mmel $");
+#endif
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: k_sincos.h,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+
+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/lib/libm/src/k_sincosf.h
diff -u /dev/null src/lib/libm/src/k_sincosf.h:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/k_sincosf.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,48 @@
+/*-
+ * ====================================================
+ * 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.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/k_sincosf.h 319047 2017-05-28 06:13:38Z mmel $");
+#endif
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: k_sincosf.h,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+
+/* |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/lib/libm/src/k_sincosl.h
diff -u /dev/null src/lib/libm/src/k_sincosl.h:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/k_sincosl.h Sat Aug 27 04:31:59 2022
@@ -0,0 +1,139 @@
+/*-
+ * ====================================================
+ * 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
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/k_sincosl.h 354520 2019-11-07 23:57:48Z lwhsu $");
+#endif
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: k_sincosl.h,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+
+#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
+ *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 + 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/lib/libm/src/s_sincos.c
diff -u /dev/null src/lib/libm/src/s_sincos.c:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/s_sincos.c Sat Aug 27 04:31:59 2022
@@ -0,0 +1,90 @@
+/*-
+ * ====================================================
+ * 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>
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: s_sincos.c,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/s_sincos.c 319047 2017-05-28 06:13:38Z mmel $");
+#endif
+
+#include "namespace.h"
+#include <float.h>
+
+#include "math.h"
+// #define INLINE_REM_PIO2
+#include "math_private.h"
+// #include "e_rem_pio2.c"
+#include "k_sincos.h"
+
+#ifdef __weak_alias
+__weak_alias(sincos,_sincos)
+#endif
+
+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/lib/libm/src/s_sincosf.c
diff -u /dev/null src/lib/libm/src/s_sincosf.c:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/s_sincosf.c Sat Aug 27 04:31:59 2022
@@ -0,0 +1,136 @@
+/*-
+ * ====================================================
+ * 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_sincosf.c -- float version of s_sincos.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
+ * Optimized by Bruce D. Evans.
+ * Merged s_sinf.c and s_cosf.c by Steven G. Kargl.
+ */
+
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/s_sincosf.c 319047 2017-05-28 06:13:38Z mmel $");
+#endif
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: s_sincosf.c,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+
+#include "namespace.h"
+#include <float.h>
+
+#include "math.h"
+#define INLINE_REM_PIO2F
+#include "math_private.h"
+#include "e_rem_pio2f.h"
+#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 */
+
+#ifdef __weak_alias
+__weak_alias(sincosf,_sincosf)
+#endif
+
+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_pio2fd(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/lib/libm/src/s_sincosl.c
diff -u /dev/null src/lib/libm/src/s_sincosl.c:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/lib/libm/src/s_sincosl.c Sat Aug 27 04:31:59 2022
@@ -0,0 +1,116 @@
+/*-
+ * 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>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/s_sincosl.c 319047 2017-05-28 06:13:38Z mmel $");
+#endif
+#if defined(LIBM_SCCS) && !defined(lint)
+__RCSID("$NetBSD: s_sincosl.c,v 1.1 2022/08/27 08:31:59 christos Exp $");
+#endif
+
+#include "namespace.h"
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "math.h"
+#include "math_private.h"
+#include "k_sincosl.h"
+
+#ifdef __HAVE_LONG_DOUBLE
+#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
+
+#ifdef __weak_alias
+__weak_alias(sincosl,_sincosl)
+#endif
+
+void
+sincosl(long double x, long double *sn, long double *cs)
+{
+ union ieee_ext_u z;
+ int e0;
+ long double y[2];
+
+ z.extu_ld = x;
+ z.extu_sign = 0;
+
+ ENTERV();
+
+ /* Optimize the case where x is already within range. */
+ if (z.extu_ld < M_PI_4) {
+ /*
+ * If x = +-0 or x is a subnormal number, then sin(x) = x and
+ * cos(x) = 1.
+ */
+ if (z.extu_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.extu_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();
+}
+#endif
Index: src/tests/lib/libm/t_sincos.c
diff -u /dev/null src/tests/lib/libm/t_sincos.c:1.1
--- /dev/null Sat Aug 27 04:32:00 2022
+++ src/tests/lib/libm/t_sincos.c Sat Aug 27 04:31:58 2022
@@ -0,0 +1,478 @@
+/* $NetBSD: t_sincos.c,v 1.1 2022/08/27 08:31:58 christos Exp $ */
+
+/*-
+ * Copyright (c) 2011, 2022 The NetBSD Foundation, Inc.
+ * All rights reserved.
+ *
+ * This code is derived from software contributed to The NetBSD Foundation
+ * by Jukka Ruohonen and Christos Zoulas
+ *
+ * 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 NETBSD FOUNDATION, INC. 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 FOUNDATION 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.
+ */
+
+#include <assert.h>
+#include <atf-c.h>
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+static const struct {
+ int angle;
+ double x;
+ double y;
+ float fy;
+} sin_angles[] = {
+// { -360, -6.283185307179586, 2.4492935982947064e-16, -1.7484555e-07 },
+ { -180, -3.141592653589793, -1.2246467991473532e-16, 8.7422777e-08 },
+ { -135, -2.356194490192345, -0.7071067811865476, 999 },
+// { -90, -1.570796326794897, -1.0000000000000000, 999 },
+ { -45, -0.785398163397448, -0.7071067811865472, 999 },
+ { 0, 0.000000000000000, 0.0000000000000000, 999 },
+ { 30, 0.5235987755982989, 0.5000000000000000, 999 },
+ { 45, 0.785398163397448, 0.7071067811865472, 999 },
+// { 60, 1.047197551196598, 0.8660254037844388, 999 },
+ { 90, 1.570796326794897, 1.0000000000000000, 999 },
+// { 120, 2.094395102393195, 0.8660254037844389, 999 },
+ { 135, 2.356194490192345, 0.7071067811865476, 999 },
+ { 150, 2.617993877991494, 0.5000000000000003, 999 },
+ { 180, 3.141592653589793, 1.2246467991473532e-16, -8.7422777e-08 },
+ { 270, 4.712388980384690, -1.0000000000000000, 999 },
+ { 360, 6.283185307179586, -2.4492935982947064e-16, 1.7484555e-07 },
+};
+
+static const struct {
+ int angle;
+ double x;
+ double y;
+ float fy;
+} cos_angles[] = {
+ { -180, -3.141592653589793, -1.0000000000000000, 999 },
+ { -135, -2.356194490192345, -0.7071067811865476, 999 },
+// { -90, -1.5707963267948966, 6.123233995736766e-17, -4.3711388e-08 },
+// { -90, -1.5707963267948968, -1.6081226496766366e-16, -4.3711388e-08 },
+ { -45, -0.785398163397448, 0.7071067811865478, 999 },
+ { 0, 0.000000000000000, 1.0000000000000000, 999 },
+ { 30, 0.5235987755982989, 0.8660254037844386, 999 },
+ { 45, 0.785398163397448, 0.7071067811865478, 999 },
+// { 60, 1.0471975511965976, 0.5000000000000001, 999 },
+// { 60, 1.0471975511965979, 0.4999999999999999, 999 },
+ { 90, 1.570796326794897, -3.8285686989269494e-16, -4.3711388e-08 },
+// { 120, 2.0943951023931953, -0.4999999999999998, 999 },
+// { 120, 2.0943951023931957, -0.5000000000000002, 999 },
+ { 135, 2.356194490192345, -0.7071067811865476, 999 },
+ { 150, 2.617993877991494, -0.8660254037844386, 999 },
+ { 180, 3.141592653589793, -1.0000000000000000, 999 },
+ { 270, 4.712388980384690, -1.8369701987210297e-16, 1.1924881e-08 },
+ { 360, 6.283185307179586, 1.0000000000000000, 999 },
+};
+
+#ifdef __HAVE_LONG_DOUBLE
+/*
+ * sincosl(3)
+ */
+ATF_TC(sincosl_angles);
+ATF_TC_HEAD(sincosl_angles, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test some selected angles");
+}
+
+ATF_TC_BODY(sincosl_angles, tc)
+{
+ /*
+ * XXX The given data is for double, so take that
+ * into account and expect less precise results..
+ */
+ const long double eps = DBL_EPSILON;
+ size_t i;
+
+ ATF_CHECK(__arraycount(sin_angles) == __arraycount(cos_angles));
+
+ for (i = 0; i < __arraycount(sin_angles); i++) {
+ ATF_CHECK_MSG(sin_angles[i].angle == cos_angles[i].angle,
+ "%zu %d %d", i, sin_angles[i].angle, cos_angles[i].angle);
+ int deg = sin_angles[i].angle;
+ ATF_CHECK_MSG(sin_angles[i].x == cos_angles[i].x,
+ "%zu %g %g", i, sin_angles[i].x, cos_angles[i].x);
+ long double theta = sin_angles[i].x;
+ long double sin_theta = sin_angles[i].y;
+ long double cos_theta = cos_angles[i].y;
+ long double s, c;
+
+ sincosl(theta, &s, &c);
+
+ if (fabsl((s - sin_theta)/sin_theta) > eps) {
+ atf_tc_fail_nonfatal("sin(%d deg = %.17Lg) = %.17Lg"
+ " != %.17Lg",
+ deg, theta, s, sin_theta);
+ }
+ if (fabsl((c - cos_theta)/cos_theta) > eps) {
+ atf_tc_fail_nonfatal("cos(%d deg = %.17Lg) = %.17Lg"
+ " != %.17Lg",
+ deg, theta, c, cos_theta);
+ }
+ }
+}
+
+ATF_TC(sincosl_nan);
+ATF_TC_HEAD(sincosl_nan, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosl(NaN) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosl_nan, tc)
+{
+ const long double x = 0.0L / 0.0L;
+ long double s, c;
+
+ sincosl(x, &s, &c);
+ ATF_CHECK(isnan(x) && isnan(s) && isnan(c));
+}
+
+ATF_TC(sincosl_inf_neg);
+ATF_TC_HEAD(sincosl_inf_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosl(-Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosl_inf_neg, tc)
+{
+ const long double x = -1.0L / 0.0L;
+ long double s, c;
+
+ sincosl(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+}
+
+ATF_TC(sincosl_inf_pos);
+ATF_TC_HEAD(sincosl_inf_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosl(+Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosl_inf_pos, tc)
+{
+ const long double x = 1.0L / 0.0L;
+ long double s, c;
+
+ sincosl(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+}
+
+
+ATF_TC(sincosl_zero_neg);
+ATF_TC_HEAD(sincosl_zero_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosl(-0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincosl_zero_neg, tc)
+{
+ const long double x = -0.0L;
+ long double s, c;
+
+ sincosl(x, &s, &c);
+ ATF_CHECK(s == 0.0 && c == 1.0);
+}
+
+ATF_TC(sincosl_zero_pos);
+ATF_TC_HEAD(sincosl_zero_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosl(+0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincosl_zero_pos, tc)
+{
+ const long double x = 0.0L;
+ long double s, c;
+
+ sincosl(x, &s, &c);
+ ATF_CHECK(s == 0.0 && c == 1.0);
+}
+#endif
+
+/*
+ * sincos(3)
+ */
+ATF_TC(sincos_angles);
+ATF_TC_HEAD(sincos_angles, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test some selected angles");
+}
+
+ATF_TC_BODY(sincos_angles, tc)
+{
+ const double eps = DBL_EPSILON;
+ size_t i;
+
+ for (i = 0; i < __arraycount(sin_angles); i++) {
+ ATF_CHECK_MSG(sin_angles[i].angle == cos_angles[i].angle,
+ "%zu %d %d", i, sin_angles[i].angle, cos_angles[i].angle);
+ int deg = sin_angles[i].angle;
+ ATF_CHECK_MSG(sin_angles[i].x == cos_angles[i].x,
+ "%zu %g %g", i, sin_angles[i].x, cos_angles[i].x);
+ double theta = sin_angles[i].x;
+ double sin_theta = sin_angles[i].y;
+ double cos_theta = cos_angles[i].y;
+ double s, c;
+
+ sincos(theta, &s, &c);
+
+ if (fabs((s - sin_theta)/sin_theta) > eps) {
+ atf_tc_fail_nonfatal("sin(%d deg = %.17g) = %.17g"
+ " != %.17g",
+ deg, theta, s, sin_theta);
+ }
+ if (fabs((c - cos_theta)/cos_theta) > eps) {
+ atf_tc_fail_nonfatal("cos(%d deg = %.17g) = %.17g"
+ " != %.17g",
+ deg, theta, c, cos_theta);
+ }
+ }
+}
+
+ATF_TC(sincos_nan);
+ATF_TC_HEAD(sincos_nan, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincos(NaN) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincos_nan, tc)
+{
+ const double x = 0.0L / 0.0L;
+ double s, c;
+
+ sincos(x, &s, &c);
+ ATF_CHECK(isnan(x) && isnan(s) && isnan(c));
+}
+
+ATF_TC(sincos_inf_neg);
+ATF_TC_HEAD(sincos_inf_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincos(-Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincos_inf_neg, tc)
+{
+ const double x = -1.0L / 0.0L;
+ double s, c;
+
+ sincos(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+}
+
+ATF_TC(sincos_inf_pos);
+ATF_TC_HEAD(sincos_inf_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincos(+Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincos_inf_pos, tc)
+{
+ const double x = 1.0L / 0.0L;
+ double s, c;
+
+ sincos(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+}
+
+
+ATF_TC(sincos_zero_neg);
+ATF_TC_HEAD(sincos_zero_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincos(-0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincos_zero_neg, tc)
+{
+ const double x = -0.0L;
+ double s, c;
+
+ sincos(x, &s, &c);
+ ATF_CHECK(s == 0 && c == 1.0);
+}
+
+ATF_TC(sincos_zero_pos);
+ATF_TC_HEAD(sincos_zero_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cos(+0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincos_zero_pos, tc)
+{
+ const double x = 0.0L;
+ double s, c;
+
+ sincos(x, &s, &c);
+ ATF_CHECK(s == 0 && c == 1.0);
+}
+
+/*
+ * sincosf(3)
+ */
+ATF_TC(sincosf_angles);
+ATF_TC_HEAD(sincosf_angles, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test some selected angles");
+}
+
+ATF_TC_BODY(sincosf_angles, tc)
+{
+ const float eps = FLT_EPSILON;
+ size_t i;
+
+ for (i = 0; i < __arraycount(sin_angles); i++) {
+ ATF_CHECK_MSG(sin_angles[i].angle == cos_angles[i].angle,
+ "%zu %d %d", i, sin_angles[i].angle, cos_angles[i].angle);
+ int deg = sin_angles[i].angle;
+ ATF_CHECK_MSG(sin_angles[i].x == cos_angles[i].x,
+ "%zu %g %g", i, sin_angles[i].x, cos_angles[i].x);
+ float theta = sin_angles[i].x;
+ float sin_theta = sin_angles[i].fy;
+ float cos_theta = cos_angles[i].fy;
+ float s, c;
+
+ sincosf(theta, &s, &c);
+ if (cos_theta == 999)
+ cos_theta = cos_angles[i].y;
+ if (sin_theta == 999)
+ sin_theta = sin_angles[i].y;
+
+ if (fabs((s - sin_theta)/sin_theta) > eps) {
+ atf_tc_fail_nonfatal("sin(%d deg = %.8g) = %.8g"
+ " != %.8g",
+ deg, theta, s, sin_theta);
+ }
+ if (fabs((c - cos_theta)/cos_theta) > eps) {
+ atf_tc_fail_nonfatal("cos(%d deg = %.8g) = %.8g"
+ " != %.8g",
+ deg, theta, c, cos_theta);
+ }
+ }
+}
+
+ATF_TC(sincosf_nan);
+ATF_TC_HEAD(sincosf_nan, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cosf(NaN) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosf_nan, tc)
+{
+ const float x = 0.0L / 0.0L;
+ float s, c;
+
+ sincosf(x, &s, &c);
+ ATF_CHECK(isnan(x) && isnan(s) && isnan(c));
+}
+
+ATF_TC(sincosf_inf_neg);
+ATF_TC_HEAD(sincosf_inf_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cosf(-Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosf_inf_neg, tc)
+{
+ const float x = -1.0L / 0.0L;
+ float s, c;
+
+ sincosf(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+
+}
+
+ATF_TC(sincosf_inf_pos);
+ATF_TC_HEAD(sincosf_inf_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosf(+Inf) == (NaN, NaN)");
+}
+
+ATF_TC_BODY(sincosf_inf_pos, tc)
+{
+ const float x = 1.0L / 0.0L;
+ float s, c;
+
+ sincosf(x, &s, &c);
+ ATF_CHECK(isnan(s) && isnan(c));
+}
+
+
+ATF_TC(sincosf_zero_neg);
+ATF_TC_HEAD(sincosf_zero_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosf(-0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincosf_zero_neg, tc)
+{
+ const float x = -0.0L;
+ float s, c;
+
+ sincosf(x, &s, &c);
+
+ ATF_CHECK(s == 0.0 && c == 1.0);
+}
+
+ATF_TC(sincosf_zero_pos);
+ATF_TC_HEAD(sincosf_zero_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sincosf(+0.0) == (0.0, 1.0)");
+}
+
+ATF_TC_BODY(sincosf_zero_pos, tc)
+{
+ const float x = 0.0L;
+
+ float s, c;
+
+ sincosf(x, &s, &c);
+
+ ATF_CHECK(s == 0 && c == 1.0);
+}
+
+ATF_TP_ADD_TCS(tp)
+{
+#ifdef __HAVE_LONG_DOUBLE
+ ATF_TP_ADD_TC(tp, sincosl_angles);
+ ATF_TP_ADD_TC(tp, sincosl_nan);
+ ATF_TP_ADD_TC(tp, sincosl_inf_neg);
+ ATF_TP_ADD_TC(tp, sincosl_inf_pos);
+ ATF_TP_ADD_TC(tp, sincosl_zero_neg);
+ ATF_TP_ADD_TC(tp, sincosl_zero_pos);
+#endif
+
+ ATF_TP_ADD_TC(tp, sincos_angles);
+ ATF_TP_ADD_TC(tp, sincos_nan);
+ ATF_TP_ADD_TC(tp, sincos_inf_neg);
+ ATF_TP_ADD_TC(tp, sincos_inf_pos);
+ ATF_TP_ADD_TC(tp, sincos_zero_neg);
+ ATF_TP_ADD_TC(tp, sincos_zero_pos);
+
+ ATF_TP_ADD_TC(tp, sincosf_angles);
+ ATF_TP_ADD_TC(tp, sincosf_nan);
+ ATF_TP_ADD_TC(tp, sincosf_inf_neg);
+ ATF_TP_ADD_TC(tp, sincosf_inf_pos);
+ ATF_TP_ADD_TC(tp, sincosf_zero_neg);
+ ATF_TP_ADD_TC(tp, sincosf_zero_pos);
+
+ return atf_no_error();
+}
