Index: head/lib/msun/Makefile =================================================================== --- head/lib/msun/Makefile +++ head/lib/msun/Makefile @@ -73,7 +73,8 @@ s_nexttowardf.c s_remquo.c s_remquof.c \ s_rint.c s_rintf.c s_round.c s_roundf.c \ s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \ - s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \ + s_signgam.c s_significand.c s_significandf.c s_sin.c \ + s_sincos.c s_sincosf.c s_sinf.c \ s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \ w_cabs.c w_cabsf.c w_drem.c w_dremf.c @@ -104,7 +105,8 @@ s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \ s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \ s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \ - s_scalbnl.c s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c + s_scalbnl.c s_sinl.c s_sincosl.c \ + s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c .endif # C99 complex functions @@ -137,7 +139,8 @@ fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \ lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \ nextafter.3 remainder.3 rint.3 \ - round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \ + round.3 scalbn.3 signbit.3 sin.3 sincos.3 \ + sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \ complex.3 MLINKS+=acos.3 acosf.3 acos.3 acosl.3 @@ -215,6 +218,7 @@ MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3 MLINKS+=sin.3 sinf.3 sin.3 sinl.3 +MLINKS+=sincos.3 sincosf.3 sin.3 sincosl.3 MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \ sqrt.3 sqrtl.3 Index: head/lib/msun/Symbol.map =================================================================== --- head/lib/msun/Symbol.map +++ head/lib/msun/Symbol.map @@ -294,4 +294,7 @@ casinl; catanl; catanhl; + sincos; + sincosf; + sincosl; }; Index: head/lib/msun/man/sincos.3 =================================================================== --- head/lib/msun/man/sincos.3 +++ head/lib/msun/man/sincos.3 @@ -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 . + Index: head/lib/msun/src/k_sincos.h =================================================================== --- head/lib/msun/src/k_sincos.h +++ head/lib/msun/src/k_sincos.h @@ -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 +__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)); +} Index: head/lib/msun/src/k_sincosf.h =================================================================== --- head/lib/msun/src/k_sincosf.h +++ head/lib/msun/src/k_sincosf.h @@ -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 +__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; +} Index: head/lib/msun/src/k_sincosl.h =================================================================== --- head/lib/msun/src/k_sincosl.h +++ head/lib/msun/src/k_sincosl.h @@ -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 +__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 Index: head/lib/msun/src/math.h =================================================================== --- head/lib/msun/src/math.h +++ head/lib/msun/src/math.h @@ -500,6 +500,9 @@ #if __BSD_VISIBLE long double lgammal_r(long double, int *); +void sincos(double, double *, double *); +void sincosf(float, float *, float *); +void sincosl(long double, long double *, long double *); #endif __END_DECLS Index: head/lib/msun/src/math_private.h =================================================================== --- head/lib/msun/src/math_private.h +++ head/lib/msun/src/math_private.h @@ -306,9 +306,21 @@ 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. */ Index: head/lib/msun/src/s_sincos.c =================================================================== --- head/lib/msun/src/s_sincos.c +++ head/lib/msun/src/s_sincos.c @@ -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 +__FBSDID("$FreeBSD$"); + +#include + +#include "math.h" +#define INLINE_REM_PIO2 +#include "math_private.h" +#include "e_rem_pio2.c" +#include "k_sincos.h" + +void +sincos(double x, double *sn, double *cs) +{ + double y[2]; + int32_t n, ix; + + /* High word of x. */ + GET_HIGH_WORD(ix, x); + + /* |x| ~< pi/4 */ + ix &= 0x7fffffff; + if (ix <= 0x3fe921fb) { + if (ix < 0x3e400000) { /* |x| < 2**-27 */ + if ((int)x == 0) { /* Generate inexact. */ + *sn = x; + *cs = 1; + return; + } + } + __kernel_sincos(x, 0, 0, sn, cs); + return; + } + + /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */ + if (ix >= 0x7ff00000) { + *sn = x - x; + *cs = x - x; + return; + } + + /* Argument reduction. */ + n = __ieee754_rem_pio2(x, y); + + switch(n & 3) { + case 0: + __kernel_sincos(y[0], y[1], 1, sn, cs); + break; + case 1: + __kernel_sincos(y[0], y[1], 1, cs, sn); + *cs = -*cs; + break; + case 2: + __kernel_sincos(y[0], y[1], 1, sn, cs); + *sn = -*sn; + *cs = -*cs; + break; + default: + __kernel_sincos(y[0], y[1], 1, cs, sn); + *sn = -*sn; + } +} + +#if (LDBL_MANT_DIG == 53) +__weak_reference(sincos, sincosl); +#endif Index: head/lib/msun/src/s_sincosf.c =================================================================== --- head/lib/msun/src/s_sincosf.c +++ head/lib/msun/src/s_sincosf.c @@ -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 +__FBSDID("$FreeBSD$"); + +#include + +#include "math.h" +#define INLINE_REM_PIO2F +#include "math_private.h" +#include "e_rem_pio2f.c" +#include "k_sincosf.h" + +/* Small multiples of pi/2 rounded to double precision. */ +static const double +p1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */ +p2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */ +p3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */ +p4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */ + +void +sincosf(float x, float *sn, float *cs) +{ + float c, s; + double y; + int32_t n, hx, ix; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7fffffff; + + if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + if ((int)x == 0) { + *sn = x; /* x with inexact if x != 0 */ + *cs = 1; + return; + } + } + __kernel_sincosdf(x, sn, cs); + return; + } + + if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */ + if (ix <= 0x4016cbe3) { /* |x| ~<= 3pi/4 */ + if (hx > 0) { + __kernel_sincosdf(x - p1pio2, cs, sn); + *cs = -*cs; + } else { + __kernel_sincosdf(x + p1pio2, cs, sn); + *sn = -*sn; + } + } else { + if (hx > 0) + __kernel_sincosdf(x - p2pio2, sn, cs); + else + __kernel_sincosdf(x + p2pio2, sn, cs); + *sn = -*sn; + *cs = -*cs; + } + return; + } + + if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */ + if (ix <= 0x40afeddf) { /* |x| ~<= 7*pi/4 */ + if (hx > 0) { + __kernel_sincosdf(x - p3pio2, cs, sn); + *sn = -*sn; + } else { + __kernel_sincosdf(x + p3pio2, cs, sn); + *cs = -*cs; + } + } else { + if (hx > 0) + __kernel_sincosdf(x - p4pio2, sn, cs); + else + __kernel_sincosdf(x + p4pio2, sn, cs); + } + return; + } + + /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */ + if (ix >= 0x7f800000) { + *sn = x - x; + *cs = x - x; + return; + } + + /* Argument reduction. */ + n = __ieee754_rem_pio2f(x, &y); + __kernel_sincosdf(y, &s, &c); + + switch(n & 3) { + case 0: + *sn = s; + *cs = c; + break; + case 1: + *sn = c; + *cs = -s; + break; + case 2: + *sn = -s; + *cs = -c; + break; + default: + *sn = -c; + *cs = s; + } +} + + Index: head/lib/msun/src/s_sincosl.c =================================================================== --- head/lib/msun/src/s_sincosl.c +++ head/lib/msun/src/s_sincosl.c @@ -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 +__FBSDID("$FreeBSD$"); + +#include +#ifdef __i386__ +#include +#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(); +}