master
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/atanf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/atan.c
6
7const std = @import("../std.zig");
8const math = std.math;
9const mem = std.mem;
10const expect = std.testing.expect;
11
12/// Returns the arc-tangent of x.
13///
14/// Special Cases:
15/// - atan(+-0) = +-0
16/// - atan(+-inf) = +-pi/2
17pub fn atan(x: anytype) @TypeOf(x) {
18 const T = @TypeOf(x);
19 return switch (T) {
20 f32 => atan32(x),
21 f64 => atan64(x),
22 else => @compileError("atan not implemented for " ++ @typeName(T)),
23 };
24}
25
26fn atan32(x_: f32) f32 {
27 const atanhi = [_]f32{
28 4.6364760399e-01, // atan(0.5)hi
29 7.8539812565e-01, // atan(1.0)hi
30 9.8279368877e-01, // atan(1.5)hi
31 1.5707962513e+00, // atan(inf)hi
32 };
33
34 const atanlo = [_]f32{
35 5.0121582440e-09, // atan(0.5)lo
36 3.7748947079e-08, // atan(1.0)lo
37 3.4473217170e-08, // atan(1.5)lo
38 7.5497894159e-08, // atan(inf)lo
39 };
40
41 const aT = [_]f32{
42 3.3333328366e-01,
43 -1.9999158382e-01,
44 1.4253635705e-01,
45 -1.0648017377e-01,
46 6.1687607318e-02,
47 };
48
49 var x = x_;
50 var ix: u32 = @as(u32, @bitCast(x));
51 const sign = ix >> 31;
52 ix &= 0x7FFFFFFF;
53
54 // |x| >= 2^26
55 if (ix >= 0x4C800000) {
56 if (math.isNan(x)) {
57 return x;
58 } else {
59 const z = atanhi[3] + 0x1.0p-120;
60 return if (sign != 0) -z else z;
61 }
62 }
63
64 var id: ?usize = undefined;
65
66 // |x| < 0.4375
67 if (ix < 0x3EE00000) {
68 // |x| < 2^(-12)
69 if (ix < 0x39800000) {
70 if (ix < 0x00800000) {
71 mem.doNotOptimizeAway(x * x);
72 }
73 return x;
74 }
75 id = null;
76 } else {
77 x = @abs(x);
78 // |x| < 1.1875
79 if (ix < 0x3F980000) {
80 // 7/16 <= |x| < 11/16
81 if (ix < 0x3F300000) {
82 id = 0;
83 x = (2.0 * x - 1.0) / (2.0 + x);
84 }
85 // 11/16 <= |x| < 19/16
86 else {
87 id = 1;
88 x = (x - 1.0) / (x + 1.0);
89 }
90 } else {
91 // |x| < 2.4375
92 if (ix < 0x401C0000) {
93 id = 2;
94 x = (x - 1.5) / (1.0 + 1.5 * x);
95 }
96 // 2.4375 <= |x| < 2^26
97 else {
98 id = 3;
99 x = -1.0 / x;
100 }
101 }
102 }
103
104 const z = x * x;
105 const w = z * z;
106 const s1 = z * (aT[0] + w * (aT[2] + w * aT[4]));
107 const s2 = w * (aT[1] + w * aT[3]);
108
109 if (id) |id_value| {
110 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x);
111 return if (sign != 0) -zz else zz;
112 } else {
113 return x - x * (s1 + s2);
114 }
115}
116
117fn atan64(x_: f64) f64 {
118 const atanhi = [_]f64{
119 4.63647609000806093515e-01, // atan(0.5)hi
120 7.85398163397448278999e-01, // atan(1.0)hi
121 9.82793723247329054082e-01, // atan(1.5)hi
122 1.57079632679489655800e+00, // atan(inf)hi
123 };
124
125 const atanlo = [_]f64{
126 2.26987774529616870924e-17, // atan(0.5)lo
127 3.06161699786838301793e-17, // atan(1.0)lo
128 1.39033110312309984516e-17, // atan(1.5)lo
129 6.12323399573676603587e-17, // atan(inf)lo
130 };
131
132 const aT = [_]f64{
133 3.33333333333329318027e-01,
134 -1.99999999998764832476e-01,
135 1.42857142725034663711e-01,
136 -1.11111104054623557880e-01,
137 9.09088713343650656196e-02,
138 -7.69187620504482999495e-02,
139 6.66107313738753120669e-02,
140 -5.83357013379057348645e-02,
141 4.97687799461593236017e-02,
142 -3.65315727442169155270e-02,
143 1.62858201153657823623e-02,
144 };
145
146 var x = x_;
147 const ux: u64 = @bitCast(x);
148 var ix: u32 = @intCast(ux >> 32);
149 const sign = ix >> 31;
150 ix &= 0x7FFFFFFF;
151
152 // |x| >= 2^66
153 if (ix >= 0x44100000) {
154 if (math.isNan(x)) {
155 return x;
156 } else {
157 const z = atanhi[3] + 0x1.0p-120;
158 return if (sign != 0) -z else z;
159 }
160 }
161
162 var id: ?usize = undefined;
163
164 // |x| < 0.4375
165 if (ix < 0x3FDC0000) {
166 // |x| < 2^(-27)
167 if (ix < 0x3E400000) {
168 if (ix < 0x00100000) {
169 mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
170 }
171 return x;
172 }
173 id = null;
174 } else {
175 x = @abs(x);
176 // |x| < 1.1875
177 if (ix < 0x3FF30000) {
178 // 7/16 <= |x| < 11/16
179 if (ix < 0x3FE60000) {
180 id = 0;
181 x = (2.0 * x - 1.0) / (2.0 + x);
182 }
183 // 11/16 <= |x| < 19/16
184 else {
185 id = 1;
186 x = (x - 1.0) / (x + 1.0);
187 }
188 } else {
189 // |x| < 2.4375
190 if (ix < 0x40038000) {
191 id = 2;
192 x = (x - 1.5) / (1.0 + 1.5 * x);
193 }
194 // 2.4375 <= |x| < 2^66
195 else {
196 id = 3;
197 x = -1.0 / x;
198 }
199 }
200 }
201
202 const z = x * x;
203 const w = z * z;
204 const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
205 const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
206
207 if (id) |id_value| {
208 const zz = atanhi[id_value] - ((x * (s1 + s2) - atanlo[id_value]) - x);
209 return if (sign != 0) -zz else zz;
210 } else {
211 return x - x * (s1 + s2);
212 }
213}
214
215test atan {
216 try expect(@as(u32, @bitCast(atan(@as(f32, 0.2)))) == @as(u32, @bitCast(atan32(0.2))));
217 try expect(atan(@as(f64, 0.2)) == atan64(0.2));
218}
219
220test atan32 {
221 const epsilon = 0.000001;
222
223 try expect(math.approxEqAbs(f32, atan32(0.2), 0.197396, epsilon));
224 try expect(math.approxEqAbs(f32, atan32(-0.2), -0.197396, epsilon));
225 try expect(math.approxEqAbs(f32, atan32(0.3434), 0.330783, epsilon));
226 try expect(math.approxEqAbs(f32, atan32(0.8923), 0.728545, epsilon));
227 try expect(math.approxEqAbs(f32, atan32(1.5), 0.982794, epsilon));
228}
229
230test atan64 {
231 const epsilon = 0.000001;
232
233 try expect(math.approxEqAbs(f64, atan64(0.2), 0.197396, epsilon));
234 try expect(math.approxEqAbs(f64, atan64(-0.2), -0.197396, epsilon));
235 try expect(math.approxEqAbs(f64, atan64(0.3434), 0.330783, epsilon));
236 try expect(math.approxEqAbs(f64, atan64(0.8923), 0.728545, epsilon));
237 try expect(math.approxEqAbs(f64, atan64(1.5), 0.982794, epsilon));
238}
239
240test "atan32.special" {
241 const epsilon = 0.000001;
242
243 try expect(math.isPositiveZero(atan32(0.0)));
244 try expect(math.isNegativeZero(atan32(-0.0)));
245 try expect(math.approxEqAbs(f32, atan32(math.inf(f32)), math.pi / 2.0, epsilon));
246 try expect(math.approxEqAbs(f32, atan32(-math.inf(f32)), -math.pi / 2.0, epsilon));
247}
248
249test "atan64.special" {
250 const epsilon = 0.000001;
251
252 try expect(math.isPositiveZero(atan64(0.0)));
253 try expect(math.isNegativeZero(atan64(-0.0)));
254 try expect(math.approxEqAbs(f64, atan64(math.inf(f64)), math.pi / 2.0, epsilon));
255 try expect(math.approxEqAbs(f64, atan64(-math.inf(f64)), -math.pi / 2.0, epsilon));
256}