Add sincos{,f,l} from FreeBSD
This commit is contained in:
christos
parent
ffd100a962
commit
5dd94c5e95
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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.3 atan2.3 atanh.3 \
|
|||
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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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 .
|
||||
|
||||
|
|
@ -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 (Berkeley) 6/4/93";
|
|||
#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);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,80 @@
|
|||
/* e_rem_pio2f.c -- float version of e_rem_pio2.c
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* 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;
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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));
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
||||
/*
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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, ian@cygnus.com.
|
||||
* 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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
}
|
||||
Loading…
Reference in New Issue