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