1+ /* @(#)s_sincos.c 5.1 13/07/15 */
2+ /*
3+ * ====================================================
4+ * Copyright (C) 2013 Elliot Saba. All rights reserved.
5+ *
6+ * Developed at the University of Washington.
7+ * Permission to use, copy, modify, and distribute this
8+ * software is freely granted, provided that this notice
9+ * is preserved.
10+ * ====================================================
11+ */
12+ #include "cdefs-compat.h"
13+
14+ /* sincos(x, s, c)
15+ * Several applications need sine and cosine of the same
16+ * angle x. This function computes both at the same time,
17+ * and stores the results in *sin and *cos.
18+ *
19+ * kernel function:
20+ * __kernel_sin ... sine function on [-pi/4,pi/4]
21+ * __kernel_cos ... cose function on [-pi/4,pi/4]
22+ * __ieee754_rem_pio2 ... argument reduction routine
23+ *
24+ * Method.
25+ * Borrow liberally from s_sin.c and s_cos.c, merging
26+ * efforts where applicable and returning their values in
27+ * appropriate variables, thereby slightly reducing the
28+ * amount of work relative to just calling sin/cos(x)
29+ * separately
30+ *
31+ * Special cases:
32+ * Let trig be any of sin, cos, or tan.
33+ * sincos(+-INF, s, c) is NaN, with signals;
34+ * sincos(NaN, s, c) is that NaN;
35+ */
36+
37+ #include <float.h>
38+
39+ #include "openlibm.h"
40+ //#define INLINE_REM_PIO2
41+ #include "math_private.h"
42+ //#include "e_rem_pio2.c"
43+
44+ /* Constants used in polynomial approximation of sin/cos */
45+ static const double
46+ one = 1.00000000000000000000e+00 , /* 0x3FF00000, 0x00000000 */
47+ half = 5.00000000000000000000e-01 , /* 0x3FE00000, 0x00000000 */
48+ S1 = -1.66666666666666324348e-01 , /* 0xBFC55555, 0x55555549 */
49+ S2 = 8.33333333332248946124e-03 , /* 0x3F811111, 0x1110F8A6 */
50+ S3 = -1.98412698298579493134e-04 , /* 0xBF2A01A0, 0x19C161D5 */
51+ S4 = 2.75573137070700676789e-06 , /* 0x3EC71DE3, 0x57B1FE7D */
52+ S5 = -2.50507602534068634195e-08 , /* 0xBE5AE5E6, 0x8A2B9CEB */
53+ S6 = 1.58969099521155010221e-10 , /* 0x3DE5D93A, 0x5ACFD57C */
54+ C1 = 4.16666666666666019037e-02 , /* 0x3FA55555, 0x5555554C */
55+ C2 = -1.38888888888741095749e-03 , /* 0xBF56C16C, 0x16C15177 */
56+ C3 = 2.48015872894767294178e-05 , /* 0x3EFA01A0, 0x19CB1590 */
57+ C4 = -2.75573143513906633035e-07 , /* 0xBE927E4F, 0x809C52AD */
58+ C5 = 2.08757232129817482790e-09 , /* 0x3E21EE9E, 0xBDB4B1C4 */
59+ C6 = -1.13596475577881948265e-11 ; /* 0xBDA8FAE9, 0xBE8838D4 */
60+
61+ void
62+ __kernel_sincos ( double x , double y , int iy , double * k_s , double * k_c )
63+ {
64+ /* Inline calculation of sin/cos, as we can save
65+ some work, and we will always need to calculate
66+ both values, no matter the result of switch */
67+ double z , w , r , v , hz ;
68+ z = x * x ;
69+ w = z * z ;
70+
71+ /* cos-specific computation; equivalent to calling
72+ __kernel_cos(x,y) and storing in k_c*/
73+ r = z * (C1 + z * (C2 + z * C3 )) + w * w * (C4 + z * (C5 + z * C6 ));
74+ hz = 0.5 * z ;
75+ v = one - hz ;
76+
77+ * k_c = v + (((one - v )- hz ) + (z * r - x * y ));
78+
79+ /* sin-specific computation; equivalent to calling
80+ __kernel_sin(x,y,1) and storing in k_s*/
81+ r = S2 + z * (S3 + z * S4 ) + z * w * (S5 + z * S6 );
82+ v = z * x ;
83+ if (iy == 0 )
84+ * k_s = x + v * (S1 + z * r );
85+ else
86+ * k_s = x - ((z * (half * y - v * r )- y )- v * S1 );
87+ }
88+
89+ void
90+ sincos (double x , double * s , double * c )
91+ {
92+ double y [2 ];
93+ int32_t ix ;
94+
95+ /* Store high word of x in ix */
96+ GET_HIGH_WORD (ix ,x );
97+
98+ /* |x| ~< pi/4 */
99+ ix &= 0x7fffffff ;
100+ if (ix <= 0x3fe921fb ) {
101+ /* Check for small x for sin and cos */
102+ if (ix < 0x3e46a09e ) {
103+ /* Check for exact zero */
104+ if ( (int )x == 0 ) {
105+ * s = x ;
106+ * c = 1.0 ;
107+ return ;
108+ }
109+ }
110+ /* Call kernel function with 0 extra */
111+ __kernel_sincos (x ,0.0 ,0 , s , c );
112+ } else if ( ix >= 0x7ff00000 ) {
113+ /* sincos(Inf or NaN) is NaN */
114+ * s = x - x ;
115+ * c = x - x ;
116+ }
117+
118+ /*argument reduction needed*/
119+ else {
120+ double k_c , k_s ;
121+ printf ( "bleh?\n" );
122+
123+ /* Calculate remainer, then sub out to kernel */
124+ int32_t n = __ieee754_rem_pio2 (x ,y );
125+ __kernel_sincos ( y [0 ], y [1 ], 1 , & k_s , & k_c );
126+
127+ /* Figure out permutation of sin/cos outputs to true outputs */
128+ switch (n & 3 ) {
129+ case 0 :
130+ * c = k_c ;
131+ * s = k_s ;
132+ break ;
133+ case 1 :
134+ * c = - k_s ;
135+ * s = k_c ;
136+ break ;
137+ case 2 :
138+ * c = - k_c ;
139+ * s = - k_s ;
140+ break ;
141+ default :
142+ * c = k_s ;
143+ * s = - k_c ;
144+ break ;
145+ }
146+ }
147+ }
148+
149+ #if (LDBL_MANT_DIG == 53 )
150+ __weak_reference (sincos , sincosl );
151+ #endif
0 commit comments