1 | |
---|
2 | /* @(#)z_sine.c 1.0 98/08/13 */ |
---|
3 | /****************************************************************** |
---|
4 | * The following routines are coded directly from the algorithms |
---|
5 | * and coefficients given in "Software Manual for the Elementary |
---|
6 | * Functions" by William J. Cody, Jr. and William Waite, Prentice |
---|
7 | * Hall, 1980. |
---|
8 | ******************************************************************/ |
---|
9 | |
---|
10 | /* |
---|
11 | FUNCTION |
---|
12 | <<sin>>, <<cos>>, <<sine>>, <<sinf>>, <<cosf>>, <<sinef>>---sine or cosine |
---|
13 | INDEX |
---|
14 | sin |
---|
15 | INDEX |
---|
16 | sinf |
---|
17 | INDEX |
---|
18 | cos |
---|
19 | INDEX |
---|
20 | cosf |
---|
21 | SYNOPSIS |
---|
22 | #include <math.h> |
---|
23 | double sin(double <[x]>); |
---|
24 | float sinf(float <[x]>); |
---|
25 | double cos(double <[x]>); |
---|
26 | float cosf(float <[x]>); |
---|
27 | |
---|
28 | DESCRIPTION |
---|
29 | <<sin>> and <<cos>> compute (respectively) the sine and cosine |
---|
30 | of the argument <[x]>. Angles are specified in radians. |
---|
31 | RETURNS |
---|
32 | The sine or cosine of <[x]> is returned. |
---|
33 | |
---|
34 | PORTABILITY |
---|
35 | <<sin>> and <<cos>> are ANSI C. |
---|
36 | <<sinf>> and <<cosf>> are extensions. |
---|
37 | |
---|
38 | QUICKREF |
---|
39 | sin ansi pure |
---|
40 | sinf - pure |
---|
41 | */ |
---|
42 | |
---|
43 | /****************************************************************** |
---|
44 | * sine |
---|
45 | * |
---|
46 | * Input: |
---|
47 | * x - floating point value |
---|
48 | * cosine - indicates cosine value |
---|
49 | * |
---|
50 | * Output: |
---|
51 | * Sine of x. |
---|
52 | * |
---|
53 | * Description: |
---|
54 | * This routine calculates sines and cosines. |
---|
55 | * |
---|
56 | *****************************************************************/ |
---|
57 | |
---|
58 | #include "fdlibm.h" |
---|
59 | #include "zmath.h" |
---|
60 | |
---|
61 | #ifndef _DOUBLE_IS_32BITS |
---|
62 | |
---|
63 | static const double HALF_PI = 1.57079632679489661923; |
---|
64 | static const double ONE_OVER_PI = 0.31830988618379067154; |
---|
65 | static const double r[] = { -0.16666666666666665052, |
---|
66 | 0.83333333333331650314e-02, |
---|
67 | -0.19841269841201840457e-03, |
---|
68 | 0.27557319210152756119e-05, |
---|
69 | -0.25052106798274584544e-07, |
---|
70 | 0.16058936490371589114e-09, |
---|
71 | -0.76429178068910467734e-12, |
---|
72 | 0.27204790957888846175e-14 }; |
---|
73 | |
---|
74 | double |
---|
75 | sine (double x, |
---|
76 | int cosine) |
---|
77 | { |
---|
78 | int sgn, N; |
---|
79 | double y, XN, g, R, res; |
---|
80 | double YMAX = 210828714.0; |
---|
81 | |
---|
82 | switch (numtest (x)) |
---|
83 | { |
---|
84 | case NAN: |
---|
85 | errno = EDOM; |
---|
86 | return (x); |
---|
87 | case INF: |
---|
88 | errno = EDOM; |
---|
89 | return (z_notanum.d); |
---|
90 | } |
---|
91 | |
---|
92 | /* Use sin and cos properties to ease computations. */ |
---|
93 | if (cosine) |
---|
94 | { |
---|
95 | sgn = 1; |
---|
96 | y = fabs (x) + HALF_PI; |
---|
97 | } |
---|
98 | else |
---|
99 | { |
---|
100 | if (x < 0.0) |
---|
101 | { |
---|
102 | sgn = -1; |
---|
103 | y = -x; |
---|
104 | } |
---|
105 | else |
---|
106 | { |
---|
107 | sgn = 1; |
---|
108 | y = x; |
---|
109 | } |
---|
110 | } |
---|
111 | |
---|
112 | /* Check for values of y that will overflow here. */ |
---|
113 | if (y > YMAX) |
---|
114 | { |
---|
115 | errno = ERANGE; |
---|
116 | return (x); |
---|
117 | } |
---|
118 | |
---|
119 | /* Calculate the exponent. */ |
---|
120 | if (y < 0.0) |
---|
121 | N = (int) (y * ONE_OVER_PI - 0.5); |
---|
122 | else |
---|
123 | N = (int) (y * ONE_OVER_PI + 0.5); |
---|
124 | XN = (double) N; |
---|
125 | |
---|
126 | if (N & 1) |
---|
127 | sgn = -sgn; |
---|
128 | |
---|
129 | if (cosine) |
---|
130 | XN -= 0.5; |
---|
131 | |
---|
132 | y = fabs (x) - XN * __PI; |
---|
133 | |
---|
134 | if (-z_rooteps < y && y < z_rooteps) |
---|
135 | res = y; |
---|
136 | |
---|
137 | else |
---|
138 | { |
---|
139 | g = y * y; |
---|
140 | |
---|
141 | /* Calculate the Taylor series. */ |
---|
142 | R = (((((((r[6] * g + r[5]) * g + r[4]) * g + r[3]) * g + r[2]) * g + r[1]) * g + r[0]) * g); |
---|
143 | |
---|
144 | /* Finally, compute the result. */ |
---|
145 | res = y + y * R; |
---|
146 | } |
---|
147 | |
---|
148 | res *= sgn; |
---|
149 | |
---|
150 | return (res); |
---|
151 | } |
---|
152 | |
---|
153 | #endif /* _DOUBLE_IS_32BITS */ |
---|