Black Lives Matter. Support the Equal Justice Initiative.

# Source file src/math/sin.go

## Documentation: math

```     1  // Copyright 2011 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
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  //	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  	)
125  	// special cases
126  	switch {
127  	case IsNaN(x) || IsInf(x, 0):
128  		return NaN()
129  	}
130
131  	// make argument positive
132  	sign := false
133  	x = Abs(x)
134
135  	var j uint64
136  	var y, z float64
137  	if x >= reduceThreshold {
138  		j, z = trigReduce(x)
139  	} else {
140  		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
141  		y = float64(j)           // integer part of x/(Pi/4), as float
142
143  		// map zeros to origin
144  		if j&1 == 1 {
145  			j++
146  			y++
147  		}
148  		j &= 7                               // octant modulo 2Pi radians (360 degrees)
149  		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
150  	}
151
152  	if j > 3 {
153  		j -= 4
154  		sign = !sign
155  	}
156  	if j > 1 {
157  		sign = !sign
158  	}
159
160  	zz := z * z
161  	if j == 1 || j == 2 {
162  		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
163  	} else {
164  		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])
165  	}
166  	if sign {
167  		y = -y
168  	}
169  	return y
170  }
171
172  // Sin returns the sine of the radian argument x.
173  //
174  // Special cases are:
175  //	Sin(±0) = ±0
176  //	Sin(±Inf) = NaN
177  //	Sin(NaN) = NaN
178  func Sin(x float64) float64
179
180  func sin(x float64) float64 {
181  	const (
182  		PI4A = 7.85398125648498535156e-1  // 0x3fe921fb40000000, Pi/4 split into three parts
183  		PI4B = 3.77489470793079817668e-8  // 0x3e64442d00000000,
184  		PI4C = 2.69515142907905952645e-15 // 0x3ce8469898cc5170,
185  	)
186  	// special cases
187  	switch {
188  	case x == 0 || IsNaN(x):
189  		return x // return ±0 || NaN()
190  	case IsInf(x, 0):
191  		return NaN()
192  	}
193
194  	// make argument positive but save the sign
195  	sign := false
196  	if x < 0 {
197  		x = -x
198  		sign = true
199  	}
200
201  	var j uint64
202  	var y, z float64
203  	if x >= reduceThreshold {
204  		j, z = trigReduce(x)
205  	} else {
206  		j = uint64(x * (4 / Pi)) // integer part of x/(Pi/4), as integer for tests on the phase angle
207  		y = float64(j)           // integer part of x/(Pi/4), as float
208
209  		// map zeros to origin
210  		if j&1 == 1 {
211  			j++
212  			y++
213  		}
214  		j &= 7                               // octant modulo 2Pi radians (360 degrees)
215  		z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
216  	}
217  	// reflect in x axis
218  	if j > 3 {
219  		sign = !sign
220  		j -= 4
221  	}
222  	zz := z * z
223  	if j == 1 || j == 2 {
224  		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])
225  	} else {
226  		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
227  	}
228  	if sign {
229  		y = -y
230  	}
231  	return y
232  }
233
```

View as plain text