]>
Commit | Line | Data |
---|---|---|
8faf50e0 XL |
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 | ||
dc9dc135 XL |
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. | |
48663c56 | 67 | #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] |
8faf50e0 XL |
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 | } | |
48663c56 XL |
108 | } else if ix < 0x40038000 { |
109 | /* |x| < 2.4375 */ | |
110 | x = (x - 1.5) / (1. + 1.5 * x); | |
111 | 2 | |
8faf50e0 | 112 | } else { |
48663c56 XL |
113 | /* 2.4375 <= |x| < 2^66 */ |
114 | x = -1. / x; | |
115 | 3 | |
8faf50e0 XL |
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), | |
48663c56 XL |
152 | ] |
153 | .iter() | |
8faf50e0 XL |
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 | } |