1 | module math |
2 | |
3 | // The original C code, the long comment, and the constants below were |
4 | // from http://netlib.sandia.gov/cephes/cmath/atan.c, available from |
5 | // http://www.netlib.org/cephes/ctgz. |
6 | // The go code is a version of the original C. |
7 | // |
8 | // atan.c |
9 | // Inverse circular tangent (arctangent) |
10 | // |
11 | // SYNOPSIS: |
12 | // double x, y, atan() |
13 | // y = atan( x ) |
14 | // |
15 | // DESCRIPTION: |
16 | // Returns radian angle between -pi/2.0 and +pi/2.0 whose tangent is x. |
17 | // |
18 | // Range reduction is from three intervals into the interval from zero to 0.66. |
19 | // The approximant uses a rational function of degree 4/5 of the form |
20 | // x + x**3 P(x)/Q(x). |
21 | // |
22 | // ACCURACY: |
23 | // Relative error: |
24 | // arithmetic domain # trials peak rms |
25 | // DEC -10, 10 50000 2.4e-17 8.3e-18 |
26 | // IEEE -10, 10 10^6 1.8e-16 5.0e-17 |
27 | // |
28 | // Cephes Math Library Release 2.8: June, 2000 |
29 | // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier |
30 | // |
31 | // The readme file at http://netlib.sandia.gov/cephes/ says: |
32 | // Some software in this archive may be from the book _Methods and |
33 | // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster |
34 | // International, 1989) or from the Cephes Mathematical Library, a |
35 | // commercial product. In either event, it is copyrighted by the author. |
36 | // What you see here may be used freely but it comes with no support or |
37 | // guarantee. |
38 | // |
39 | // The two known misprints in the book are repaired here in the |
40 | // source listings for the gamma function and the incomplete beta |
41 | // integral. |
42 | // |
43 | // Stephen L. Moshier |
44 | // [email protected] |
45 | // pi/2.0 = PIO2 + morebits |
46 | // tan3pio8 = tan(3*pi/8) |
47 | |
48 | const ( |
49 | morebits = 6.123233995736765886130e-17 |
50 | tan3pio8 = 2.41421356237309504880 |
51 | ) |
52 | |
53 | // xatan evaluates a series valid in the range [0, 0.66]. |
54 | [inline] |
55 | fn xatan(x f64) f64 { |
56 | xatan_p0 := -8.750608600031904122785e-01 |
57 | xatan_p1 := -1.615753718733365076637e+01 |
58 | xatan_p2 := -7.500855792314704667340e+01 |
59 | xatan_p3 := -1.228866684490136173410e+02 |
60 | xatan_p4 := -6.485021904942025371773e+01 |
61 | xatan_q0 := 2.485846490142306297962e+01 |
62 | xatan_q1 := 1.650270098316988542046e+02 |
63 | xatan_q2 := 4.328810604912902668951e+02 |
64 | xatan_q3 := 4.853903996359136964868e+02 |
65 | xatan_q4 := 1.945506571482613964425e+02 |
66 | mut z := x * x |
67 | z = z * ((((xatan_p0 * z + xatan_p1) * z + xatan_p2) * z + xatan_p3) * z + xatan_p4) / (((((z + |
68 | xatan_q0) * z + xatan_q1) * z + xatan_q2) * z + xatan_q3) * z + xatan_q4) |
69 | z = x * z + x |
70 | return z |
71 | } |
72 | |
73 | // satan reduces its argument (known to be positive) |
74 | // to the range [0, 0.66] and calls xatan. |
75 | [inline] |
76 | fn satan(x f64) f64 { |
77 | if x <= 0.66 { |
78 | return xatan(x) |
79 | } |
80 | if x > math.tan3pio8 { |
81 | return pi / 2.0 - xatan(1.0 / x) + f64(math.morebits) |
82 | } |
83 | return pi / 4 + xatan((x - 1.0) / (x + 1.0)) + 0.5 * f64(math.morebits) |
84 | } |
85 | |
86 | // atan returns the arctangent, in radians, of x. |
87 | // |
88 | // special cases are: |
89 | // atan(±0) = ±0 |
90 | // atan(±inf) = ±pi/2.0 |
91 | pub fn atan(x f64) f64 { |
92 | if x == 0 { |
93 | return x |
94 | } |
95 | if x > 0 { |
96 | return satan(x) |
97 | } |
98 | return -satan(-x) |
99 | } |
100 | |
101 | // atan2 returns the arc tangent of y/x, using |
102 | // the signs of the two to determine the quadrant |
103 | // of the return value. |
104 | // |
105 | // special cases are (in order): |
106 | // atan2(y, nan) = nan |
107 | // atan2(nan, x) = nan |
108 | // atan2(+0, x>=0) = +0 |
109 | // atan2(-0, x>=0) = -0 |
110 | // atan2(+0, x<=-0) = +pi |
111 | // atan2(-0, x<=-0) = -pi |
112 | // atan2(y>0, 0) = +pi/2.0 |
113 | // atan2(y<0, 0) = -pi/2.0 |
114 | // atan2(+inf, +inf) = +pi/4 |
115 | // atan2(-inf, +inf) = -pi/4 |
116 | // atan2(+inf, -inf) = 3pi/4 |
117 | // atan2(-inf, -inf) = -3pi/4 |
118 | // atan2(y, +inf) = 0 |
119 | // atan2(y>0, -inf) = +pi |
120 | // atan2(y<0, -inf) = -pi |
121 | // atan2(+inf, x) = +pi/2.0 |
122 | // atan2(-inf, x) = -pi/2.0 |
123 | pub fn atan2(y f64, x f64) f64 { |
124 | // special cases |
125 | if is_nan(y) || is_nan(x) { |
126 | return nan() |
127 | } |
128 | if y == 0.0 { |
129 | if x >= 0 && !signbit(x) { |
130 | return copysign(0, y) |
131 | } |
132 | return copysign(pi, y) |
133 | } |
134 | if x == 0.0 { |
135 | return copysign(pi / 2.0, y) |
136 | } |
137 | if is_inf(x, 0) { |
138 | if is_inf(x, 1) { |
139 | if is_inf(y, 0) { |
140 | return copysign(pi / 4, y) |
141 | } |
142 | return copysign(0, y) |
143 | } |
144 | if is_inf(y, 0) { |
145 | return copysign(3.0 * pi / 4.0, y) |
146 | } |
147 | return copysign(pi, y) |
148 | } |
149 | if is_inf(y, 0) { |
150 | return copysign(pi / 2.0, y) |
151 | } |
152 | // Call atan and determine the quadrant. |
153 | q := atan(y / x) |
154 | if x < 0 { |
155 | if q <= 0 { |
156 | return q + pi |
157 | } |
158 | return q - pi |
159 | } |
160 | return q |
161 | } |
162 | |
163 | /* |
164 | Floating-point arcsine and arccosine. |
165 | |
166 | They are implemented by computing the arctangent |
167 | after appropriate range reduction. |
168 | */ |
169 | |
170 | // asin returns the arcsine, in radians, of x. |
171 | // |
172 | // special cases are: |
173 | // asin(±0) = ±0 |
174 | // asin(x) = nan if x < -1 or x > 1 |
175 | pub fn asin(x_ f64) f64 { |
176 | mut x := x_ |
177 | if x == 0.0 { |
178 | return x // special case |
179 | } |
180 | mut sign := false |
181 | if x < 0.0 { |
182 | x = -x |
183 | sign = true |
184 | } |
185 | if x > 1.0 { |
186 | return nan() // special case |
187 | } |
188 | mut temp := sqrt(1.0 - x * x) |
189 | if x > 0.7 { |
190 | temp = pi / 2.0 - satan(temp / x) |
191 | } else { |
192 | temp = satan(x / temp) |
193 | } |
194 | if sign { |
195 | temp = -temp |
196 | } |
197 | return temp |
198 | } |
199 | |
200 | // acos returns the arccosine, in radians, of x. |
201 | // |
202 | // special case is: |
203 | // acos(x) = nan if x < -1 or x > 1 |
204 | [inline] |
205 | pub fn acos(x f64) f64 { |
206 | if x < -1.0 || x > 1.0 { |
207 | return nan() |
208 | } |
209 | if x > 0.5 { |
210 | return f64(2.0) * asin(sqrt(0.5 - 0.5 * x)) |
211 | } |
212 | mut z := pi / f64(4.0) - asin(x) |
213 | z = z + math.morebits |
214 | z = z + pi / f64(4.0) |
215 | return z |
216 | } |