# Source file src/math/sin.go

```     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  //
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
```

View as plain text