d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Implement sincos, sincosf, and sincosl.

Browse Source 
The primary benefit of these functions is that argument
reduction is done once instead of twice in independent
calls to sin() and cos().

* lib/msun/Makefile:
  . Add s_sincos[fl].c to the build.
  . Add sincos.3 documentation.
  . Add appropriate MLINKS.

* lib/msun/Symbol.map:
  . Expose sincos[fl] symbols in dynamic libm.so.

* lib/msun/man/sincos.3:
  . Documentation for sincos[fl].

* lib/msun/src/k_sincos.h:
  . Kernel for sincos() function.  This merges the individual kernels
    for sin() and cos().  The merger offered an opportunity to re-arrange
    the individual kernels for better performance.

* lib/msun/src/k_sincosf.h:
   . Kernel for sincosf() function.  This merges the individual kernels
     for sinf() and cosf(). The merger offered an opportunity to re-arrange
     the individual kernels for better performance.

* lib/msun/src/k_sincosl.h:
   . Kernel for sincosl() function.  This merges the individual kernels
     for sinl() and cosl(). The merger offered an opportunity to re-arrange
     the individual kernels for better performance.

* lib/msun/src/math.h:
  . Add prototytpes for sincos[fl]().

* lib/msun/src/math_private.h:
  . Add RETURNV macros.  This is needed to reset fpsetprec on I386
    hardware for a function with type void.

* lib/msun/src/s_sincos.c:
  . Implementation of sincos() where sin() and cos() were merged into
    one routine and possibly re-arranged for better performance.

* lib/msun/src/s_sincosf.c:
  . Implementation of sincosf() where sinf() and cosf() were merged into
    one routine and possibly re-arranged for better performance.

* lib/msun/src/s_sincosl.c:
  . Implementation of sincosl() where sinl() and cosl() were merged into
    one routine and possibly re-arranged for better performance.

PR:		215977, 218300
Submitted by:	Steven G. Kargl <sgk@troutmask.apl.washington.edu>
MFC after:	1 month
Differential Revision:	https://reviews.freebsd.org/D10765
This commit is contained in:
mmel 2017-05-28 06:13:38 +00:00
parent 299229a5a6
commit d02641556c
11 changed files with 648 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
lib/msun/Makefile
View File

@ -73,7 +73,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
s_nexttowardf.c s_remquo.c s_remquof.c \
s_rint.c s_rintf.c s_round.c s_roundf.c \
s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
s_signgam.c s_significand.c s_significandf.c s_sin.c \
s_sincos.c s_sincosf.c s_sinf.c \
s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
w_cabs.c w_cabsf.c w_drem.c w_dremf.c
@ -104,7 +105,8 @@ COMMON_SRCS+= catrigl.c \
s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
s_scalbnl.c s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
s_scalbnl.c s_sinl.c s_sincosl.c \
s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
.endif
# C99 complex functions
@ -137,7 +139,8 @@ MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
nextafter.3 remainder.3 rint.3 \
round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
round.3 scalbn.3 signbit.3 sin.3 sincos.3 \
sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
complex.3
MLINKS+=acos.3 acosf.3 acos.3 acosl.3
@ -215,6 +218,7 @@ MLINKS+=round.3 roundf.3 round.3 roundl.3
MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
MLINKS+=sin.3 sinf.3 sin.3 sinl.3
MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3
MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
sqrt.3 sqrtl.3

3
lib/msun/Symbol.map
View File

@ -294,4 +294,7 @@ FBSD_1.5 {
casinl;
catanl;
catanhl;
sincos;
sincosf;
sincosl;
};

82
lib/msun/man/sincos.3 Normal file
View File

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

52
lib/msun/src/k_sincos.h Normal file
View File

@ -0,0 +1,52 @@
/*-
* ====================================================
* 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>
__FBSDID("$FreeBSD$");
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));
}

43
lib/msun/src/k_sincosf.h Normal file
View File

@ -0,0 +1,43 @@
/*-
* ====================================================
* 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>
__FBSDID("$FreeBSD$");
/* |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;
}

134
lib/msun/src/k_sincosl.h Normal file
View File

@ -0,0 +1,134 @@
/*-
* ====================================================
* 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>
__FBSDID("$FreeBSD$");
#if LDBL_MANT_DIG == 64 /* ld80 version of k_sincosl.c. */
#if defined(__amd64__) || defined(__i386__)
/* Long double constants are slow on these arches, and broken on i386. */
static const volatile double
C1hi = 0.041666666666666664, /* 0x15555555555555.0p-57 */
C1lo = 2.2598839032744733e-18, /* 0x14d80000000000.0p-111 */
S1hi = -0.16666666666666666, /* -0x15555555555555.0p-55 */
S1lo = -9.2563760475949941e-18; /* -0x15580000000000.0p-109 */
#define S1 ((long double)S1hi + S1lo)
#define C1 ((long double)C1hi + C1lo)
#else
static const long double
C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
S1 = -0.166666666666666666671L, /* -0xaaaaaaaaaaaaaaab.0p-66 */
#endif
static const double
C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */
C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */
C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */
C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */
C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */
C7 = 4.7383039476436467e-14, /* 0x1aac9d9af5c43e.0p-97 */
S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
static inline void
__kernel_sincosl(long double x, long double y, int iy, long double *sn,
long double *cs)
{
long double hz, r, v, w, z;
z = x * x;
v = z * x;
/*
* XXX Replace Horner scheme with an algorithm suitable for CPUs
* with more complex pipelines.
*/
r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * S8)))));
if (iy == 0)
*sn = x + v * (S1 + z * r);
else
*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
hz = z / 2;
w = 1 - hz;
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
z * C7))))));
*cs = w + (((1 - w) - hz) + (z * r - x * y));
}
#elif LDBL_MANT_DIG == 113 /* ld128 version of k_sincosl.c. */
static const long double
C1 = 0.04166666666666666666666666666666658424671L,
C2 = -0.001388888888888888888888888888863490893732L,
C3 = 0.00002480158730158730158730158600795304914210L,
C4 = -0.2755731922398589065255474947078934284324e-6L,
C5 = 0.2087675698786809897659225313136400793948e-8L,
C6 = -0.1147074559772972315817149986812031204775e-10L,
C7 = 0.4779477332386808976875457937252120293400e-13L,
S1 = -0.16666666666666666666666666666666666606732416116558L,
S2 = 0.0083333333333333333333333333333331135404851288270047L,
S3 = -0.00019841269841269841269841269839935785325638310428717L,
S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
static const double
C8 = -0.1561920696721507929516718307820958119868e-15,
C9 = 0.4110317413744594971475941557607804508039e-18,
C10 = -0.8896592467191938803288521958313920156409e-21,
C11 = 0.1601061435794535138244346256065192782581e-23,
S9 = -0.82206352458348947812512122163446202498005154296863e-17,
S10 = 0.19572940011906109418080609928334380560135358385256e-19,
S11 = -0.38680813379701966970673724299207480965452616911420e-22,
S12 = 0.64038150078671872796678569586315881020659912139412e-25;
static inline void
__kernel_sincosl(long double x, long double y, int iy, long double *sn,
long double *cs)
{
long double hz, r, v, w, z;
z = x * x;
v = z * x;
/*
* XXX Replace Horner scheme with an algorithm suitable for CPUs
* with more complex pipelines.
*/
r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
if (iy == 0)
*sn = x + v * (S1 + z * r);
else
*cs = x - ((z * (y / 2 - v * r) - y) - v * S1);
hz = z / 2;
w = 1 - hz;
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
*cs = w + (((1 - w) - hz) + (z * r - x * y));
}
#else
#error "Unsupported long double format"
#endif

3
lib/msun/src/math.h
View File

@ -500,6 +500,9 @@ long double truncl(long double);
#if __BSD_VISIBLE
long double lgammal_r(long double, int *);
void sincos(double, double *, double *);
void sincosf(float, float *, float *);
void sincosl(long double, long double *, long double *);
#endif
__END_DECLS

14
lib/msun/src/math_private.h
View File

@ -306,9 +306,21 @@ do { \
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(x)
#define ENTERI()
#define RETURNI(x) RETURNF(x)
#define ENTERV()
#define RETURNV() return
#endif
/* Default return statement if hack*_t() is not used. */

80
lib/msun/src/s_sincos.c Normal file
View File

@ -0,0 +1,80 @@
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* s_sin.c and s_cos.c merged by Steven G. Kargl. Descriptions of the
* algorithms are contained in the original files.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <float.h>
#include "math.h"
#define INLINE_REM_PIO2
#include "math_private.h"
#include "e_rem_pio2.c"
#include "k_sincos.h"
void
sincos(double x, double *sn, double *cs)
{
double y[2];
int32_t n, ix;
/* High word of x. */
GET_HIGH_WORD(ix, x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if (ix <= 0x3fe921fb) {
if (ix < 0x3e400000) { /* |x| < 2**-27 */
if ((int)x == 0) { /* Generate inexact. */
*sn = x;
*cs = 1;
return;
}
}
__kernel_sincos(x, 0, 0, sn, cs);
return;
}
/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
if (ix >= 0x7ff00000) {
*sn = x - x;
*cs = x - x;
return;
}
/* Argument reduction. */
n = __ieee754_rem_pio2(x, y);
switch(n & 3) {
case 0:
__kernel_sincos(y[0], y[1], 1, sn, cs);
break;
case 1:
__kernel_sincos(y[0], y[1], 1, cs, sn);
*cs = -*cs;
break;
case 2:
__kernel_sincos(y[0], y[1], 1, sn, cs);
*sn = -*sn;
*cs = -*cs;
break;
default:
__kernel_sincos(y[0], y[1], 1, cs, sn);
*sn = -*sn;
}
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(sincos, sincosl);
#endif

126
lib/msun/src/s_sincosf.c Normal file
View File

@ -0,0 +1,126 @@
/*-
* ====================================================
* 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>
__FBSDID("$FreeBSD$");
#include <float.h>
#include "math.h"
#define INLINE_REM_PIO2F
#include "math_private.h"
#include "e_rem_pio2f.c"
#include "k_sincosf.h"
/* Small multiples of pi/2 rounded to double precision. */
static const double
p1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
p2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
p3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
p4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
void
sincosf(float x, float *sn, float *cs)
{
float c, s;
double y;
int32_t n, hx, ix;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7fffffff;
if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */
if (ix < 0x39800000) { /* |x| < 2**-12 */
if ((int)x == 0) {
*sn = x; /* x with inexact if x != 0 */
*cs = 1;
return;
}
}
__kernel_sincosdf(x, sn, cs);
return;
}
if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */
if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */
if (hx > 0) {
__kernel_sincosdf(x - p1pio2, cs, sn);
*cs = -*cs;
} else {
__kernel_sincosdf(x + p1pio2, cs, sn);
*sn = -*sn;
}
} else {
if (hx > 0)
__kernel_sincosdf(x - p2pio2, sn, cs);
else
__kernel_sincosdf(x + p2pio2, sn, cs);
*sn = -*sn;
*cs = -*cs;
}
return;
}
if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */
if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */
if (hx > 0) {
__kernel_sincosdf(x - p3pio2, cs, sn);
*sn = -*sn;
} else {
__kernel_sincosdf(x + p3pio2, cs, sn);
*cs = -*cs;
}
} else {
if (hx > 0)
__kernel_sincosdf(x - p4pio2, sn, cs);
else
__kernel_sincosdf(x + p4pio2, sn, cs);
}
return;
}
/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
if (ix >= 0x7f800000) {
*sn = x - x;
*cs = x - x;
return;
}
/* Argument reduction. */
n = __ieee754_rem_pio2f(x, &y);
__kernel_sincosdf(y, &s, &c);
switch(n & 3) {
case 0:
*sn = s;
*cs = c;
break;
case 1:
*sn = c;
*cs = -s;
break;
case 2:
*sn = -s;
*cs = -c;
break;
default:
*sn = -c;
*cs = s;
}
}

105
lib/msun/src/s_sincosl.c Normal file
View File

@ -0,0 +1,105 @@
/*-
* Copyright (c) 2007, 2010-2013 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
* s_sinl.c and s_cosl.c merged by Steven G. Kargl.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <float.h>
#ifdef __i386__
#include <ieeefp.h>
#endif
#include "math.h"
#include "math_private.h"
#include "k_sincosl.h"
#if LDBL_MANT_DIG == 64
#include "../ld80/e_rem_pio2l.h"
#elif LDBL_MANT_DIG == 113
#include "../ld128/e_rem_pio2l.h"
#else
#error "Unsupported long double format"
#endif
void
sincosl(long double x, long double *sn, long double *cs)
{
union IEEEl2bits z;
int e0, sgn;
long double y[2];
z.e = x;
sgn = z.bits.sign;
z.bits.sign = 0;
ENTERV();
/* Optimize the case where x is already within range. */
if (z.e < M_PI_4) {
/*
* If x = +-0 or x is a subnormal number, then sin(x) = x and
* cos(x) = 1.
*/
if (z.bits.exp == 0) {
*sn = x;
*cs = 1;
} else
__kernel_sincosl(x, 0, 0, sn, cs);
RETURNV();
}
/* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
if (z.bits.exp == 32767) {
*sn = x - x;
*cs = x - x;
RETURNV();
}
/* Range reduction. */
e0 = __ieee754_rem_pio2l(x, y);
switch (e0 & 3) {
case 0:
__kernel_sincosl(y[0], y[1], 1, sn, cs);
break;
case 1:
__kernel_sincosl(y[0], y[1], 1, cs, sn);
*cs = -*cs;
break;
case 2:
__kernel_sincosl(y[0], y[1], 1, sn, cs);
*sn = -*sn;
*cs = -*cs;
break;
default:
__kernel_sincosl(y[0], y[1], 1, cs, sn);
*sn = -*sn;
}
RETURNV();
}