1 // Copyright 2011 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
4 5 package math
6 7 /*
8 Floating-point sine and cosine.
9 */
10 11 // The original C code, the long comment, and the constants
12 // below were from http://netlib.sandia.gov/cephes/cmath/sin.c,
13 // available from http://www.netlib.org/cephes/cmath.tgz.
14 // The go code is a simplified version of the original C.
15 //
16 // sin.c
17 //
18 // Circular sine
19 //
20 // SYNOPSIS:
21 //
22 // double x, y, sin();
23 // y = sin( x );
24 //
25 // DESCRIPTION:
26 //
27 // Range reduction is into intervals of pi/4. The reduction error is nearly
28 // eliminated by contriving an extended precision modular arithmetic.
29 //
30 // Two polynomial approximating functions are employed.
31 // Between 0 and pi/4 the sine is approximated by
32 // x + x**3 P(x**2).
33 // Between pi/4 and pi/2 the cosine is represented as
34 // 1 - x**2 Q(x**2).
35 //
36 // ACCURACY:
37 //
38 // Relative error:
39 // arithmetic domain # trials peak rms
40 // DEC 0, 10 150000 3.0e-17 7.8e-18
41 // IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
42 //
43 // Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9. The loss
44 // is not gradual, but jumps suddenly to about 1 part in 10e7. Results may
45 // be meaningless for x > 2**49 = 5.6e14.
46 //
47 // cos.c
48 //
49 // Circular cosine
50 //
51 // SYNOPSIS:
52 //
53 // double x, y, cos();
54 // y = cos( x );
55 //
56 // DESCRIPTION:
57 //
58 // Range reduction is into intervals of pi/4. The reduction error is nearly
59 // eliminated by contriving an extended precision modular arithmetic.
60 //
61 // Two polynomial approximating functions are employed.
62 // Between 0 and pi/4 the cosine is approximated by
63 // 1 - x**2 Q(x**2).
64 // Between pi/4 and pi/2 the sine is represented as
65 // x + x**3 P(x**2).
66 //
67 // ACCURACY:
68 //
69 // Relative error:
70 // arithmetic domain # trials peak rms
71 // IEEE -1.07e9,+1.07e9 130000 2.1e-16 5.4e-17
72 // DEC 0,+1.07e9 17000 3.0e-17 7.2e-18
73 //
74 // Cephes Math Library Release 2.8: June, 2000
75 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
76 //
77 // The readme file at http://netlib.sandia.gov/cephes/ says:
78 // Some software in this archive may be from the book _Methods and
79 // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
80 // International, 1989) or from the Cephes Mathematical Library, a
81 // commercial product. In either event, it is copyrighted by the author.
82 // What you see here may be used freely but it comes with no support or
83 // guarantee.
84 //
85 // The two known misprints in the book are repaired here in the
86 // source listings for the gamma function and the incomplete beta
87 // integral.
88 //
89 // Stephen L. Moshier
90 // moshier@na-net.ornl.gov
91 92 // sin coefficients
93 var _sin = [...]float64{
94 1.58962301576546568060e-10, // 0x3de5d8fd1fd19ccd
95 -2.50507477628578072866e-8, // 0xbe5ae5e5a9291f5d
96 2.75573136213857245213e-6, // 0x3ec71de3567d48a1
97 -1.98412698295895385996e-4, // 0xbf2a01a019bfdf03
98 8.33333333332211858878e-3, // 0x3f8111111110f7d0
99 -1.66666666666666307295e-1, // 0xbfc5555555555548
100 }
101 102 // cos coefficients
103 var _cos = [...]float64{
104 -1.13585365213876817300e-11, // 0xbda8fa49a0861a9b
105 2.08757008419747316778e-9, // 0x3e21ee9d7b4e3f05
106 -2.75573141792967388112e-7, // 0xbe927e4f7eac4bc6
107 2.48015872888517045348e-5, // 0x3efa01a019c844f5
108 -1.38888888888730564116e-3, // 0xbf56c16c16c14f91
109 4.16666666666665929218e-2, // 0x3fa555555555554b
110 }
111 112 // Cos returns the cosine of the radian argument x.
113 //
114 // Special cases are:
115 //
116 // Cos(±Inf) = NaN
117 // Cos(NaN) = NaN
118 func Cos(x float64) float64 {
119 if haveArchCos {
120 return archCos(x)
121 }
122 return cos(x)
123 }
124 125 func cos(x float64) float64 {
126 const (
127 PI4A = 7.85398125648498535156e-1 // 0x3fe921fb40000000, Pi/4 split into three parts
128 PI4B = 3.77489470793079817668e-8 // 0x3e64442d00000000,
129 PI4C = 2.69515142907905952645e-15 // 0x3ce8469898cc5170,
130 )
131 // special cases
132 switch {
133 case IsNaN(x) || IsInf(x, 0):
134 return NaN()
135 }
136 137 // make argument positive
138 sign := false
139 x = Abs(x)
140 141 var j uint64
142 var y, z float64
143 if x >= reduceThreshold {
144 j, z = trigReduce(x)
145 } else {
146 j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
147 y = float64(j) // integer part of x/(Pi/4), as float
148 149 // map zeros to origin
150 if j&1 == 1 {
151 j++
152 y++
153 }
154 j &= 7 // octant modulo 2Pi radians (360 degrees)
155 z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
156 }
157 158 if j > 3 {
159 j -= 4
160 sign = !sign
161 }
162 if j > 1 {
163 sign = !sign
164 }
165 166 zz := z * z
167 if j == 1 || j == 2 {
168 y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
169 } else {
170 y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
171 }
172 if sign {
173 y = -y
174 }
175 return y
176 }
177 178 // Sin returns the sine of the radian argument x.
179 //
180 // Special cases are:
181 //
182 // Sin(±0) = ±0
183 // Sin(±Inf) = NaN
184 // Sin(NaN) = NaN
185 func Sin(x float64) float64 {
186 if haveArchSin {
187 return archSin(x)
188 }
189 return sin(x)
190 }
191 192 func sin(x float64) float64 {
193 const (
194 PI4A = 7.85398125648498535156e-1 // 0x3fe921fb40000000, Pi/4 split into three parts
195 PI4B = 3.77489470793079817668e-8 // 0x3e64442d00000000,
196 PI4C = 2.69515142907905952645e-15 // 0x3ce8469898cc5170,
197 )
198 // special cases
199 switch {
200 case x == 0 || IsNaN(x):
201 return x // return ±0 || NaN()
202 case IsInf(x, 0):
203 return NaN()
204 }
205 206 // make argument positive but save the sign
207 sign := false
208 if x < 0 {
209 x = -x
210 sign = true
211 }
212 213 var j uint64
214 var y, z float64
215 if x >= reduceThreshold {
216 j, z = trigReduce(x)
217 } else {
218 j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
219 y = float64(j) // integer part of x/(Pi/4), as float
220 221 // map zeros to origin
222 if j&1 == 1 {
223 j++
224 y++
225 }
226 j &= 7 // octant modulo 2Pi radians (360 degrees)
227 z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
228 }
229 // reflect in x axis
230 if j > 3 {
231 sign = !sign
232 j -= 4
233 }
234 zz := z * z
235 if j == 1 || j == 2 {
236 y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
237 } else {
238 y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
239 }
240 if sign {
241 y = -y
242 }
243 return y
244 }
245