1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */ |
2 | /* |
3 | * ==================================================== |
4 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 | * |
6 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
7 | * Permission to use, copy, modify, and distribute this |
8 | * software is freely granted, provided that this notice |
9 | * is preserved. |
10 | * ==================================================== |
11 | */ |
12 | /* atan(x) |
13 | * Method |
14 | * 1. Reduce x to positive by atan(x) = -atan(-x). |
15 | * 2. According to the integer k=4t+0.25 chopped, t=x, the argument |
16 | * is further reduced to one of the following intervals and the |
17 | * arctangent of t is evaluated by the corresponding formula: |
18 | * |
19 | * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) |
20 | * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) ) |
21 | * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) ) |
22 | * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) ) |
23 | * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) |
24 | * |
25 | * Constants: |
26 | * The hexadecimal values are the intended ones for the following |
27 | * constants. The decimal values may be used, provided that the |
28 | * compiler will convert from decimal to binary accurately enough |
29 | * to produce the hexadecimal values shown. |
30 | */ |
31 | |
32 | use super::fabs; |
33 | use core::f64; |
34 | |
35 | const ATANHI: [f64; 4] = [ |
36 | 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ |
37 | 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ |
38 | 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ |
39 | 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ |
40 | ]; |
41 | |
42 | const ATANLO: [f64; 4] = [ |
43 | 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ |
44 | 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ |
45 | 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ |
46 | 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ |
47 | ]; |
48 | |
49 | const AT: [f64; 11] = [ |
50 | 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ |
51 | -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ |
52 | 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ |
53 | -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */ |
54 | 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */ |
55 | -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */ |
56 | 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */ |
57 | -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */ |
58 | 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */ |
59 | -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */ |
60 | 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ |
61 | ]; |
62 | |
63 | /// Arctangent (f64) |
64 | /// |
65 | /// Computes the inverse tangent (arc tangent) of the input value. |
66 | /// Returns a value in radians, in the range of -pi/2 to pi/2. |
67 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
68 | pub fn atan(x: f64) -> f64 { |
69 | let mut x = x; |
70 | let mut ix = (x.to_bits() >> 32) as u32; |
71 | let sign = ix >> 31; |
72 | ix &= 0x7fff_ffff; |
73 | if ix >= 0x4410_0000 { |
74 | if x.is_nan() { |
75 | return x; |
76 | } |
77 | |
78 | let z = ATANHI[3] + f64::from_bits(0x0380_0000); // 0x1p-120f |
79 | return if sign != 0 { -z } else { z }; |
80 | } |
81 | |
82 | let id = if ix < 0x3fdc_0000 { |
83 | /* |x| < 0.4375 */ |
84 | if ix < 0x3e40_0000 { |
85 | /* |x| < 2^-27 */ |
86 | if ix < 0x0010_0000 { |
87 | /* raise underflow for subnormal x */ |
88 | force_eval!(x as f32); |
89 | } |
90 | |
91 | return x; |
92 | } |
93 | |
94 | -1 |
95 | } else { |
96 | x = fabs(x); |
97 | if ix < 0x3ff30000 { |
98 | /* |x| < 1.1875 */ |
99 | if ix < 0x3fe60000 { |
100 | /* 7/16 <= |x| < 11/16 */ |
101 | x = (2. * x - 1.) / (2. + x); |
102 | 0 |
103 | } else { |
104 | /* 11/16 <= |x| < 19/16 */ |
105 | x = (x - 1.) / (x + 1.); |
106 | 1 |
107 | } |
108 | } else if ix < 0x40038000 { |
109 | /* |x| < 2.4375 */ |
110 | x = (x - 1.5) / (1. + 1.5 * x); |
111 | 2 |
112 | } else { |
113 | /* 2.4375 <= |x| < 2^66 */ |
114 | x = -1. / x; |
115 | 3 |
116 | } |
117 | }; |
118 | |
119 | let z = x * x; |
120 | let w = z * z; |
121 | /* break sum from i=0 to 10 AT[i]z**(i+1) into odd and even poly */ |
122 | let s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10]))))); |
123 | let s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9])))); |
124 | |
125 | if id < 0 { |
126 | return x - x * (s1 + s2); |
127 | } |
128 | |
129 | let z = i!(ATANHI, id as usize) - (x * (s1 + s2) - i!(ATANLO, id as usize) - x); |
130 | |
131 | if sign != 0 { |
132 | -z |
133 | } else { |
134 | z |
135 | } |
136 | } |
137 | |
138 | #[cfg (test)] |
139 | mod tests { |
140 | use super::atan; |
141 | use core::f64; |
142 | |
143 | #[test ] |
144 | fn sanity_check() { |
145 | for (input, answer) in [ |
146 | (3.0_f64.sqrt() / 3.0, f64::consts::FRAC_PI_6), |
147 | (1.0, f64::consts::FRAC_PI_4), |
148 | (3.0_f64.sqrt(), f64::consts::FRAC_PI_3), |
149 | (-3.0_f64.sqrt() / 3.0, -f64::consts::FRAC_PI_6), |
150 | (-1.0, -f64::consts::FRAC_PI_4), |
151 | (-3.0_f64.sqrt(), -f64::consts::FRAC_PI_3), |
152 | ] |
153 | .iter() |
154 | { |
155 | assert!( |
156 | (atan(*input) - answer) / answer < 1e-5, |
157 | " \natan( {:.4}/16) = {:.4}, actual: {}" , |
158 | input * 16.0, |
159 | answer, |
160 | atan(*input) |
161 | ); |
162 | } |
163 | } |
164 | |
165 | #[test ] |
166 | fn zero() { |
167 | assert_eq!(atan(0.0), 0.0); |
168 | } |
169 | |
170 | #[test ] |
171 | fn infinity() { |
172 | assert_eq!(atan(f64::INFINITY), f64::consts::FRAC_PI_2); |
173 | } |
174 | |
175 | #[test ] |
176 | fn minus_infinity() { |
177 | assert_eq!(atan(f64::NEG_INFINITY), -f64::consts::FRAC_PI_2); |
178 | } |
179 | |
180 | #[test ] |
181 | fn nan() { |
182 | assert!(atan(f64::NAN).is_nan()); |
183 | } |
184 | } |
185 | |