Source file src/pkg/math/sin.go
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 // [email protected]
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 x.
113 //
114 // Special cases are:
115 // Cos(±Inf) = NaN
116 // Cos(NaN) = NaN
117 func Cos(x float64) float64
118
119 func cos(x float64) float64 {
120 const (
121 PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts
122 PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000,
123 PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
124 M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
125 )
126 // special cases
127 switch {
128 case IsNaN(x) || IsInf(x, 0):
129 return NaN()
130 }
131
132 // make argument positive
133 sign := false
134 if x < 0 {
135 x = -x
136 }
137
138 j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
139 y := float64(j) // integer part of x/(Pi/4), as float
140
141 // map zeros to origin
142 if j&1 == 1 {
143 j += 1
144 y += 1
145 }
146 j &= 7 // octant modulo 2Pi radians (360 degrees)
147 if j > 3 {
148 j -= 4
149 sign = !sign
150 }
151 if j > 1 {
152 sign = !sign
153 }
154
155 z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
156 zz := z * z
157 if j == 1 || j == 2 {
158 y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
159 } else {
160 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])
161 }
162 if sign {
163 y = -y
164 }
165 return y
166 }
167
168 // Sin returns the sine of x.
169 //
170 // Special cases are:
171 // Sin(±0) = ±0
172 // Sin(±Inf) = NaN
173 // Sin(NaN) = NaN
174 func Sin(x float64) float64
175
176 func sin(x float64) float64 {
177 const (
178 PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts
179 PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000,
180 PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170,
181 M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
182 )
183 // special cases
184 switch {
185 case x == 0 || IsNaN(x):
186 return x // return ±0 || NaN()
187 case IsInf(x, 0):
188 return NaN()
189 }
190
191 // make argument positive but save the sign
192 sign := false
193 if x < 0 {
194 x = -x
195 sign = true
196 }
197
198 j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
199 y := float64(j) // integer part of x/(Pi/4), as float
200
201 // map zeros to origin
202 if j&1 == 1 {
203 j += 1
204 y += 1
205 }
206 j &= 7 // octant modulo 2Pi radians (360 degrees)
207 // reflect in x axis
208 if j > 3 {
209 sign = !sign
210 j -= 4
211 }
212
213 z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
214 zz := z * z
215 if j == 1 || j == 2 {
216 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])
217 } else {
218 y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
219 }
220 if sign {
221 y = -y
222 }
223 return y
224 }