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/acosf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/acos.c
6// https://git.musl-libc.org/cgit/musl/tree/src/math/acosl.c
7//
8// Ported from ARM-software, which is licensed under the MIT license:
9// https://github.com/ARM-software/optimized-routines/blob/master/LICENSE
10//
11// https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/acosf.c
12// https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/acos.c
13
14const std = @import("../std.zig");
15const math = std.math;
16const testing = std.testing;
17const builtin = @import("builtin");
18const native_endian = builtin.cpu.arch.endian();
19
20/// Returns the arc-cosine of x.
21///
22/// Special cases:
23/// - acos(x) = nan if x < -1 or x > 1
24pub fn acos(x: anytype) @TypeOf(x) {
25 const T = @TypeOf(x);
26 switch (@typeInfo(T)) {
27 .float => |info| switch (info.bits) {
28 16 => return acosBinary16(x),
29 32 => return acosBinary32(x),
30 64 => return acosBinary64(x),
31 80 => return acosExtended80(x),
32 128 => return acosBinary128(x),
33 else => comptime unreachable,
34 },
35 .vector => |info| switch (info.child) {
36 f32 => return acosBinary32Vec(info.len, x),
37 f64 => return acosBinary64Vec(info.len, x),
38 else => @compileError("unimplemented"),
39 },
40 else => comptime unreachable,
41 }
42}
43
44fn approxBinary16(z: f32) f32 {
45 const S0: f32 = 1.0000001e0;
46 const S1: f32 = 1.6664918e-1;
47 const S2: f32 = 7.55022e-2;
48 const S3: f32 = 3.9513987e-2;
49 const S4: f32 = 5.0883885e-2;
50 return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4)));
51}
52
53fn acosBinary16(x: f16) f16 {
54 const pio2: f32 = math.pi / 2.0;
55
56 const hx: u16 = @bitCast(x);
57 const ix: u16 = hx & 0x7fff;
58
59 // |x| >= 1 or nan
60 if (ix >= 0x3c00) {
61 if (ix == 0x3c00) {
62 if (hx >> 15 != 0) {
63 return @floatCast(2.0 * pio2 + 0x1p-120);
64 }
65 return 0.0;
66 }
67 return 0.0 / (x - x);
68 }
69
70 const xf: f32 = @floatCast(x);
71
72 // |x| < 0.5
73 if (ix < 0x3800) {
74 return @floatCast(pio2 - xf * approxBinary16(xf * xf));
75 }
76
77 // x < -0.5
78 if (hx >> 15 != 0) {
79 const z = (1.0 + xf) * 0.5;
80 const s = @sqrt(z);
81 const w = approxBinary16(z) * s;
82 return @floatCast(2.0 * (pio2 - w));
83 }
84
85 // x > 0.5
86 const z = (1.0 - xf) * 0.5;
87 const s = @sqrt(z);
88 const w = approxBinary16(z) * s;
89 return @floatCast(2.0 * w);
90}
91
92fn rationalApproxBinary32(z: f32) f32 {
93 const pS0: f32 = 1.6666586697e-01;
94 const pS1: f32 = -4.2743422091e-02;
95 const pS2: f32 = -8.6563630030e-03;
96 const qS1: f32 = -7.0662963390e-01;
97
98 // f64 is used instead of f32 to avoid
99 // a vectorization on x86_64. The vectorization
100 // causes extra floating point execeptions
101 // that are prohibited by libc-test.
102 const p: f64 = @as(f64, @floatCast(z)) * (pS0 + z * (pS1 + z * pS2));
103 const q: f64 = 1.0 + z * qS1;
104 return @floatCast(p / q);
105}
106
107fn acosBinary32(x: f32) f32 {
108 const pio2_hi: f32 = 1.5707962513e+00;
109 const pio2_lo: f32 = 7.5497894159e-08;
110
111 const hx: u32 = @bitCast(x);
112 const ix: u32 = hx & 0x7fff_ffff;
113
114 // |x| >= 1 or nan
115 if (ix >= 0x3f800000) {
116 if (ix == 0x3f800000) {
117 if (hx >> 31 != 0) {
118 return 2.0 * pio2_hi + 0x1.0p-120;
119 }
120 return 0.0;
121 }
122 return 0.0 / (x - x);
123 }
124
125 // |x| < 0.5
126 if (ix < 0x3f00_0000) {
127 // |x| < 2^(-26)
128 if (ix <= 0x3280_0000) {
129 return pio2_hi + 0x1.0p-120;
130 }
131 return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary32(x * x)));
132 }
133
134 // x < -0.5
135 if (hx >> 31 != 0) {
136 const z = (1 + x) * 0.5;
137 const s = @sqrt(z);
138 const w = rationalApproxBinary32(z) * s - pio2_lo;
139 return 2.0 * (pio2_hi - (s + w));
140 }
141
142 // x > 0.5
143 const z = (1.0 - x) * 0.5;
144 const s = @sqrt(z);
145 const hs: u32 = @bitCast(s);
146 const df: f32 = @bitCast(hs & 0xffff_f000);
147 const c = (z - df * df) / (s + df);
148 const w = rationalApproxBinary32(z) * s + c;
149 return 2.0 * (df + w);
150}
151
152fn rationalApproxBinary64(z: f64) f64 {
153 const pS0: f64 = 1.66666666666666657415e-01;
154 const pS1: f64 = -3.25565818622400915405e-01;
155 const pS2: f64 = 2.01212532134862925881e-01;
156 const pS3: f64 = -4.00555345006794114027e-02;
157 const pS4: f64 = 7.91534994289814532176e-04;
158 const pS5: f64 = 3.47933107596021167570e-05;
159 const qS1: f64 = -2.40339491173441421878e+00;
160 const qS2: f64 = 2.02094576023350569471e+00;
161 const qS3: f64 = -6.88283971605453293030e-01;
162 const qS4: f64 = 7.70381505559019352791e-02;
163
164 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
165 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
166 return p / q;
167}
168
169fn acosBinary64(x: f64) f64 {
170 const pio2_hi: f64 = 1.57079632679489655800e+00;
171 const pio2_lo: f64 = 6.12323399573676603587e-17;
172
173 const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32);
174 const ix: u32 = hx & 0x7fff_ffff;
175
176 // |x| >= 1 or nan
177 if (ix >= 0x3ff0_0000) {
178 const lx: u32 = @truncate(@as(u64, @bitCast(x)));
179 if ((ix - 0x3ff0_0000 | lx) == 0) {
180 if (hx >> 31 != 0) {
181 return 2.0 * pio2_hi + 0x1.0p-120;
182 }
183 return 0.0;
184 }
185 return 0.0 / (x - x);
186 }
187
188 // |x| < 0.5
189 if (ix < 0x3fe0_0000) {
190 // |x| < 2^(-57)
191 if (ix <= 0x3c60_0000) {
192 return pio2_hi + 0x1.0p-120;
193 }
194 return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary64(x * x)));
195 }
196
197 // x < -0.5
198 if (hx >> 31 != 0) {
199 const z = (1.0 + x) * 0.5;
200 const s = @sqrt(z);
201 const w = rationalApproxBinary64(z) * s - pio2_lo;
202 return 2 * (pio2_hi - (s + w));
203 }
204
205 // x > 0.5
206 const z = (1.0 - x) * 0.5;
207 const s = @sqrt(z);
208 const df: f64 = @bitCast(@as(u64, @bitCast(s)) & 0xffff_ffff_0000_0000);
209 const c = (z - df * df) / (s + df);
210 const w = rationalApproxBinary64(z) * s + c;
211 return 2.0 * (df + w);
212}
213
214fn rationalApproxExtended80(z: f80) f80 {
215 const pS0: f80 = 1.66666666666666666631e-01;
216 const pS1: f80 = -4.16313987993683104320e-01;
217 const pS2: f80 = 3.69068046323246813704e-01;
218 const pS3: f80 = -1.36213932016738603108e-01;
219 const pS4: f80 = 1.78324189708471965733e-02;
220 const pS5: f80 = -2.19216428382605211588e-04;
221 const pS6: f80 = -7.10526623669075243183e-06;
222 const qS1: f80 = -2.94788392796209867269e+00;
223 const qS2: f80 = 3.27309890266528636716e+00;
224 const qS3: f80 = -1.68285799854822427013e+00;
225 const qS4: f80 = 3.90699412641738801874e-01;
226 const qS5: f80 = -3.14365703596053263322e-02;
227
228 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6))))));
229 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5))));
230 return p / q;
231}
232
233fn acosExtended80(x: f80) f80 {
234 const pio2_hi: f80 = 1.57079632679489661926;
235 const pio2_lo: f80 = -2.50827880633416601173e-20;
236
237 const hx: u80 = @bitCast(x);
238 const se: u16 = @truncate(hx >> 64);
239 const e = se & 0x7fff;
240
241 // |x| >= 1 or nan
242 if (e >= 0x3fff) {
243 if (x == 1.0) {
244 return 0.0;
245 }
246 if (x == -1.0) {
247 return 2.0 * pio2_hi + 0x1p-120;
248 }
249 return 0.0 / (x - x);
250 }
251 // |x| < 0.5
252 if (e < 0x3fff - 1) {
253 if (e < 0x3fff - math.floatFractionalBits(f80)) {
254 return pio2_hi + 0x1p-120;
255 }
256 return pio2_hi - (rationalApproxExtended80(x * x) * x - pio2_lo + x);
257 }
258 // x < -0.5
259 if (se >> 15 != 0) {
260 const z = (1 + x) * 0.5;
261 const s = @sqrt(z);
262 return 2.0 * (pio2_hi - (rationalApproxExtended80(z) * s - pio2_lo + s));
263 }
264 // x > 0.5
265 const z = (1.0 - x) * 0.5;
266 const s = @sqrt(z);
267 const hs: u80 = @bitCast(s);
268 const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_0000_0000);
269 const c = (z - f * f) / (s + f);
270 return 2.0 * (rationalApproxExtended80(z) * s + c + f);
271}
272
273fn rationalApproxBinary128(z: f128) f128 {
274 const pS0: f128 = 1.66666666666666666666666666666700314e-01;
275 const pS1: f128 = -7.32816946414566252574527475428622708e-01;
276 const pS2: f128 = 1.34215708714992334609030036562143589e+00;
277 const pS3: f128 = -1.32483151677116409805070261790752040e+00;
278 const pS4: f128 = 7.61206183613632558824485341162121989e-01;
279 const pS5: f128 = -2.56165783329023486777386833928147375e-01;
280 const pS6: f128 = 4.80718586374448793411019434585413855e-02;
281 const pS7: f128 = -4.42523267167024279410230886239774718e-03;
282 const pS8: f128 = 1.44551535183911458253205638280410064e-04;
283 const pS9: f128 = -2.10558957916600254061591040482706179e-07;
284 const qS1: f128 = -4.84690167848739751544716485245697428e+00;
285 const qS2: f128 = 9.96619113536172610135016921140206980e+00;
286 const qS3: f128 = -1.13177895428973036660836798461641458e+01;
287 const qS4: f128 = 7.74004374389488266169304117714658761e+00;
288 const qS5: f128 = -3.25871986053534084709023539900339905e+00;
289 const qS6: f128 = 8.27830318881232209752469022352928864e-01;
290 const qS7: f128 = -1.18768052702942805423330715206348004e-01;
291 const qS8: f128 = 8.32600764660522313269101537926539470e-03;
292 const qS9: f128 = -1.99407384882605586705979504567947007e-04;
293
294 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9)))))))));
295 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9))))))));
296 return p / q;
297}
298
299fn acosBinary128(x: f128) f128 {
300 const pio2_hi: f128 = 1.57079632679489661923132169163975140;
301 const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35;
302
303 const hx: u128 = @bitCast(x);
304 const se: u16 = @truncate(hx >> 112);
305 const e = se & 0x7fff;
306
307 // |x| >= 1 or nan
308 if (e >= 0x3fff) {
309 if (x == 1.0) {
310 return 0.0;
311 }
312 if (x == -1.0) {
313 return 2 * pio2_hi + 0x1p-120;
314 }
315 return 0.0 / (x - x);
316 }
317 // |x| < 0.5
318 if (e < 0x3fff - 1) {
319 if (e < 0x3fff - math.floatFractionalBits(f128)) {
320 return pio2_hi + 0x1p-120;
321 }
322 return pio2_hi - (rationalApproxBinary128(x * x) * x - pio2_lo + x);
323 }
324 // x < -0.5
325 if (se >> 15 != 0) {
326 const z = (1 + x) * 0.5;
327 const s = @sqrt(z);
328 return 2 * (pio2_hi - (rationalApproxBinary128(z) * s - pio2_lo + s));
329 }
330 // x > 0.5
331 const z = (1.0 - x) * 0.5;
332 const s = @sqrt(z);
333 const hs: u128 = @bitCast(s);
334 const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000);
335 const c = (z - f * f) / (s + f);
336 return 2.0 * (rationalApproxBinary128(z) * s + c + f);
337}
338
339test "acosBinary16.special" {
340 try testing.expectApproxEqAbs(0x1.92p0, acosBinary16(0x0p+0), math.floatEpsAt(f16, 0x1.92p0));
341 try testing.expectApproxEqAbs(0x1.92p1, acosBinary16(-0x1p+0), math.floatEpsAt(f16, 0x1.92p1));
342 try testing.expectEqual(0x0p+0, acosBinary16(0x1p+0));
343 try testing.expect(math.isNan(acosBinary16(0x1.004p0)));
344 try testing.expect(math.isNan(acosBinary16(-0x1.004p0)));
345 try testing.expect(math.isNan(acosBinary16(math.inf(f16))));
346 try testing.expect(math.isNan(acosBinary16(-math.inf(f16))));
347 try testing.expect(math.isNan(acosBinary16(math.nan(f16))));
348}
349
350test "acosBinary16" {
351 try testing.expectApproxEqAbs(0x1.834p0, acosBinary16(0x1.db4p-5), math.floatEpsAt(f16, 0x1.834p0));
352 try testing.expectApproxEqAbs(0x1.d48p0, acosBinary16(-0x1.068p-2), math.floatEpsAt(f16, 0x1.d48p0));
353 try testing.expectApproxEqAbs(0x1.b7cp0, acosBinary16(-0x1.2c4p-3), math.floatEpsAt(f16, 0x1.b7cp0));
354 try testing.expectApproxEqAbs(0x1.654p0, acosBinary16(0x1.65p-3), math.floatEpsAt(f16, 0x1.654p0));
355 try testing.expectApproxEqAbs(0x1.6d8p-2, acosBinary16(0x1.dfcp-1), math.floatEpsAt(f16, 0x1.6d8p-2));
356 try testing.expectApproxEqAbs(0x1.32p1, acosBinary16(-0x1.764p-1), math.floatEpsAt(f16, 0x1.32p1));
357 try testing.expectApproxEqAbs(0x1.5b8p0, acosBinary16(0x1.b18p-3), math.floatEpsAt(f16, 0x1.5b8p0));
358 try testing.expectApproxEqAbs(0x1.668p0, acosBinary16(0x1.5acp-3), math.floatEpsAt(f16, 0x1.668p0));
359 try testing.expectApproxEqAbs(0x1.134p1, acosBinary16(-0x1.18cp-1), math.floatEpsAt(f16, 0x1.134p1));
360 try testing.expectApproxEqAbs(0x1.0dp1, acosBinary16(-0x1.03p-1), math.floatEpsAt(f16, 0x1.0dp1));
361}
362
363test "acosBinary32.special" {
364 try testing.expectApproxEqAbs(0x1.921fb6p+0, acosBinary32(0x0p+0), math.floatEpsAt(f32, 0x1.921fb6p+0));
365 try testing.expectApproxEqAbs(0x1.921fb6p+1, acosBinary32(-0x1p+0), math.floatEpsAt(f32, 0x1.921fb6p+1));
366 try testing.expectEqual(0x0p+0, acosBinary32(0x1p+0));
367 try testing.expect(math.isNan(acosBinary32(0x1.000002p+0)));
368 try testing.expect(math.isNan(acosBinary32(-0x1.000002p+0)));
369 try testing.expect(math.isNan(acosBinary32(math.inf(f32))));
370 try testing.expect(math.isNan(acosBinary32(-math.inf(f32))));
371 try testing.expect(math.isNan(acosBinary32(math.nan(f32))));
372}
373
374test "acosBinary32" {
375 try testing.expectApproxEqAbs(0x1.d7c4e6p+0, acosBinary32(-0x1.13284cp-2), math.floatEpsAt(f32, 0x1.d7c4e6p+0));
376 try testing.expectApproxEqAbs(0x1.8e6756p-1, acosBinary32(0x1.6ca8ep-1), math.floatEpsAt(f32, 0x1.8e6756p-1));
377 try testing.expectApproxEqAbs(0x1.f9d74cp-2, acosBinary32(0x1.c2ca6p-1), math.floatEpsAt(f32, 0x1.f9d74cp-2));
378 try testing.expectApproxEqAbs(0x1.26abdcp+1, acosBinary32(-0x1.55f12p-1), math.floatEpsAt(f32, 0x1.26abdcp+1));
379 try testing.expectApproxEqAbs(0x1.d85a44p+0, acosBinary32(-0x1.15679ep-2), math.floatEpsAt(f32, 0x1.d85a44p+0));
380 try testing.expectApproxEqAbs(0x1.9c2f68p+0, acosBinary32(-0x1.41e132p-5), math.floatEpsAt(f32, 0x1.9c2f68p+0));
381 try testing.expectApproxEqAbs(0x1.e881bp-1, acosBinary32(0x1.281b0ep-1), math.floatEpsAt(f32, 0x1.e881bp-1));
382 try testing.expectApproxEqAbs(0x1.1713f6p-1, acosBinary32(0x1.b5ce34p-1), math.floatEpsAt(f32, 0x1.1713f6p-1));
383 try testing.expectApproxEqAbs(0x1.bd5accp+0, acosBinary32(-0x1.583482p-3), math.floatEpsAt(f32, 0x1.bd5accp+0));
384 try testing.expectApproxEqAbs(0x1.6ce7d8p+1, acosBinary32(-0x1.ea8224p-1), math.floatEpsAt(f32, 0x1.6ce7d8p+1));
385}
386
387test "acosBinary64.special" {
388 try testing.expectApproxEqAbs(0x1.921fb54442d18p+0, acosBinary64(0x0p+0), math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
389 try testing.expectApproxEqAbs(0x1.921fb54442d18p+1, acosBinary64(-0x1p+0), math.floatEpsAt(f64, 0x1.921fb54442d18p+1));
390 try testing.expectEqual(0x0p+0, acosBinary64(0x1p+0));
391 try testing.expect(math.isNan(acosBinary64(0x1.0000000000001p+0)));
392 try testing.expect(math.isNan(acosBinary64(-0x1.0000000000001p+0)));
393 try testing.expect(math.isNan(acosBinary64(math.inf(f64))));
394 try testing.expect(math.isNan(acosBinary64(-math.inf(f64))));
395 try testing.expect(math.isNan(acosBinary64(math.nan(f64))));
396}
397
398test "acosBinary64" {
399 try testing.expectApproxEqAbs(0x1.d7c4e61020905p+0, acosBinary64(-0x1.13284b2b5006dp-2), math.floatEpsAt(f64, 0x1.d7c4e61020905p+0));
400 try testing.expectApproxEqAbs(0x1.8e6756e27c366p-1, acosBinary64(0x1.6ca8dfb825911p-1), math.floatEpsAt(f64, 0x1.8e6756e27c366p-1));
401 try testing.expectApproxEqAbs(0x1.f9d748eaf956p-2, acosBinary64(0x1.c2ca609de7505p-1), math.floatEpsAt(f64, 0x1.f9d748eaf956p-2));
402 try testing.expectApproxEqAbs(0x1.26abdc68d07aap+1, acosBinary64(-0x1.55f11fba96889p-1), math.floatEpsAt(f64, 0x1.26abdc68d07aap+1));
403 try testing.expectApproxEqAbs(0x1.d85a44ea44fe4p+0, acosBinary64(-0x1.15679e27084ddp-2), math.floatEpsAt(f64, 0x1.d85a44ea44fe4p+0));
404 try testing.expectApproxEqAbs(0x1.9c2f688eee8abp+0, acosBinary64(-0x1.41e131b093c41p-5), math.floatEpsAt(f64, 0x1.9c2f688eee8abp+0));
405 try testing.expectApproxEqAbs(0x1.e881b1d4eb2a1p-1, acosBinary64(0x1.281b0d18455f5p-1), math.floatEpsAt(f64, 0x1.e881b1d4eb2a1p-1));
406 try testing.expectApproxEqAbs(0x1.1713f567a87efp-1, acosBinary64(0x1.b5ce34a51b239p-1), math.floatEpsAt(f64, 0x1.1713f567a87efp-1));
407 try testing.expectApproxEqAbs(0x1.bd5acbe8fcc59p+0, acosBinary64(-0x1.583481079de4dp-3), math.floatEpsAt(f64, 0x1.bd5acbe8fcc59p+0));
408 try testing.expectApproxEqAbs(0x1.6ce7d66f628e5p+1, acosBinary64(-0x1.ea8223103b871p-1), math.floatEpsAt(f64, 0x1.6ce7d66f628e5p+1));
409}
410
411test "acosExtended80.special" {
412 try testing.expectApproxEqAbs(0x1.921fb54442d1846ap+0, acosExtended80(0x0p+0), math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0));
413 try testing.expectApproxEqAbs(0x1.921fb54442d1846ap+1, acosExtended80(-0x1p+0), math.floatEpsAt(f80, 0x1.921fb54442d1846ap+1));
414 try testing.expectEqual(0x0p+0, acosExtended80(0x1p+0));
415 try testing.expect(math.isNan(acosExtended80(0x1.0000000000000002p+0)));
416 try testing.expect(math.isNan(acosExtended80(-0x1.0000000000000002p+0)));
417 try testing.expect(math.isNan(acosExtended80(math.inf(f80))));
418 try testing.expect(math.isNan(acosExtended80(-math.inf(f80))));
419 try testing.expect(math.isNan(acosExtended80(math.nan(f80))));
420}
421
422test "acosExtended80" {
423 try testing.expectApproxEqAbs(0x1.86b349040d28f794p-1, acosExtended80(0x1.72068a321edc8804p-1), math.floatEpsAt(f80, 0x1.86b349040d28f794p-1));
424 try testing.expectApproxEqAbs(0x1.d4923ade73ec379cp0, acosExtended80(-0x1.06d0a467d22977ecp-2), math.floatEpsAt(f80, 0x1.d4923ade73ec379cp0));
425 try testing.expectApproxEqAbs(0x1.62e0e8898c6d04f2p0, acosExtended80(0x1.77d21385faa9798ap-3), math.floatEpsAt(f80, 0x1.62e0e8898c6d04f2p0));
426 try testing.expectApproxEqAbs(0x1.3123cbcd5dc4bd58p1, acosExtended80(-0x1.73ee3e8bc2a44dbep-1), math.floatEpsAt(f80, 0x1.3123cbcd5dc4bd58p1));
427 try testing.expectApproxEqAbs(0x1.062a6d562df2d316p0, acosExtended80(0x1.0a2dd1f6ffcf668ap-1), math.floatEpsAt(f80, 0x1.062a6d562df2d316p0));
428 try testing.expectApproxEqAbs(0x1.5ffd68b520aa55fap0, acosExtended80(0x1.8e835c490a3aff9ep-3), math.floatEpsAt(f80, 0x1.5ffd68b520aa55fap0));
429 try testing.expectApproxEqAbs(0x1.5bfe6cabda700684p0, acosExtended80(0x1.add20cdc1565064cp-3), math.floatEpsAt(f80, 0x1.5bfe6cabda700684p0));
430 try testing.expectApproxEqAbs(0x1.90fe1c993b571924p0, acosExtended80(0x1.21986d43727fca72p-8), math.floatEpsAt(f80, 0x1.90fe1c993b571924p0));
431 try testing.expectApproxEqAbs(0x1.18044ccc626e7f9ep0, acosExtended80(0x1.d61e0b3fae6a0564p-2), math.floatEpsAt(f80, 0x1.18044ccc626e7f9ep0));
432 try testing.expectApproxEqAbs(0x1.a39513b6c16532b4p0, acosExtended80(-0x1.171e7c4a41883ccap-4), math.floatEpsAt(f80, 0x1.a39513b6c16532b4p0));
433}
434
435test "acosBinary128.special" {
436 try testing.expectApproxEqAbs(0x1.921fb54442d18469898cc51701b8p0, acosBinary128(0x0p+0), math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0));
437 try testing.expectApproxEqAbs(0x1.921fb54442d18469898cc51701b8p1, acosBinary128(-0x1p+0), math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p1));
438 try testing.expectEqual(0x0p+0, acosBinary128(0x1p+0));
439 try testing.expect(math.isNan(acosBinary128(0x1.0000000000000000000000000001p0)));
440 try testing.expect(math.isNan(acosBinary128(-0x1.0000000000000000000000000001p0)));
441 try testing.expect(math.isNan(acosBinary128(math.inf(f128))));
442 try testing.expect(math.isNan(acosBinary128(-math.inf(f128))));
443 try testing.expect(math.isNan(acosBinary128(math.nan(f128))));
444}
445
446test "acosBinary128" {
447 if (builtin.cpu.arch.isSPARC()) return error.SkipZigTest;
448
449 try testing.expectApproxEqAbs(0x1.250e9a58f049eeafa99db4360c88p1, acosBinary128(-0x1.511bdb99a3c4373bedf834ef4f68p-1), math.floatEpsAt(f128, 0x1.250e9a58f049eeafa99db4360c88p1));
450 try testing.expectApproxEqAbs(0x1.2786664b1c676c99437b68590004p1, acosBinary128(-0x1.5879cc3ad6dfd2a52e9891c69808p-1), math.floatEpsAt(f128, 0x1.2786664b1c676c99437b68590004p1));
451 try testing.expectApproxEqAbs(0x1.cb190cd361c7c03a09c470b4caebp-1, acosBinary128(0x1.3f988ba64a7eb97a751c5f0b3077p-1), math.floatEpsAt(f128, 0x1.cb190cd361c7c03a09c470b4caebp-1));
452 try testing.expectApproxEqAbs(0x1.1f373be697880111758f582b1a96p1, acosBinary128(-0x1.3f2d96c7768e4c4fa02315727959p-1), math.floatEpsAt(f128, 0x1.1f373be697880111758f582b1a96p1));
453 try testing.expectApproxEqAbs(0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0, acosBinary128(0x1.fad303c2e28c1f4d8f9fd0e5686fp-2), math.floatEpsAt(f128, 0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0));
454 try testing.expectApproxEqAbs(0x1.15d4b306e16fbf9ea4f29e82b154p0, acosBinary128(0x1.ddde322bd1a2ee50c5ba30c9c617p-2), math.floatEpsAt(f128, 0x1.15d4b306e16fbf9ea4f29e82b154p0));
455 try testing.expectApproxEqAbs(0x1.49b0a0355a5539052388e8a6dc11p1, acosBinary128(-0x1.b02f6adefcbeb1d48666b827ff17p-1), math.floatEpsAt(f128, 0x1.49b0a0355a5539052388e8a6dc11p1));
456 try testing.expectApproxEqAbs(0x1.1be0b757f4cef022f5d2422b9c78p0, acosBinary128(0x1.c8581cce7cd3f6efab0fc60d9b7dp-2), math.floatEpsAt(f128, 0x1.1be0b757f4cef022f5d2422b9c78p0));
457 try testing.expectApproxEqAbs(0x1.513270e671db2d840f20b0186c2cp1, acosBinary128(-0x1.bf887b8c4e33cbef59993056f3dep-1), math.floatEpsAt(f128, 0x1.513270e671db2d840f20b0186c2cp1));
458 try testing.expectApproxEqAbs(0x1.70851a509f0e8bfbe780aa8f29f9p0, acosBinary128(0x1.0c0f600ab6f9c84c6102942044cep-3), math.floatEpsAt(f128, 0x1.70851a509f0e8bfbe780aa8f29f9p0));
459}
460
461fn acosBinary32Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f32)) @TypeOf(x) {
462 const pi: @Vector(vec_len, f32) = @splat(math.pi);
463 const pi_over_2: @Vector(vec_len, f32) = @splat(math.pi / 2.0);
464 const zero: @Vector(vec_len, f32) = @splat(0.0);
465 const half: @Vector(vec_len, f32) = @splat(0.5);
466 const neg_one: @Vector(vec_len, f32) = @splat(-1.0);
467 const two: @Vector(vec_len, f32) = @splat(2.0);
468 const c0: @Vector(vec_len, f32) = @splat(0x1.55555ep-3);
469 const c1: @Vector(vec_len, f32) = @splat(0x1.33261ap-4);
470 const c2: @Vector(vec_len, f32) = @splat(0x1.70d7dcp-5);
471 const c3: @Vector(vec_len, f32) = @splat(0x1.b059dp-6);
472 const c4: @Vector(vec_len, f32) = @splat(0x1.3af7d8p-5);
473
474 const ax = @abs(x);
475 const ax_lt_half = ax < half;
476 const is_neg = x < zero;
477 const z2 = @select(f32, ax_lt_half, x * x, @mulAdd(@Vector(vec_len, f32), -half, ax, half));
478 const z = @select(f32, ax_lt_half, ax, @sqrt(z2));
479 const z3 = z2 * z;
480 const p3_4 = @mulAdd(@Vector(vec_len, f32), z2, c4, c3);
481 const p2_4 = @mulAdd(@Vector(vec_len, f32), z2, p3_4, c2);
482 const p1_4 = @mulAdd(@Vector(vec_len, f32), z2, p2_4, c1);
483 const p0_4 = @mulAdd(@Vector(vec_len, f32), z2, p1_4, c0);
484 const p = @mulAdd(@Vector(vec_len, f32), z3, p0_4, z);
485 const mul = @select(f32, ax_lt_half, neg_one, two);
486 const add = @select(f32, ax_lt_half, pi_over_2, @select(f32, is_neg, pi, zero));
487 return @mulAdd(@Vector(vec_len, f32), mul, @select(f32, is_neg, -p, p), add);
488}
489
490fn acosBinary64Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f64)) @TypeOf(x) {
491 const pi: @Vector(vec_len, f64) = @splat(math.pi);
492 const pi_over_2: @Vector(vec_len, f64) = @splat(math.pi / 2.0);
493 const zero: @Vector(vec_len, f64) = @splat(0.0);
494 const half: @Vector(vec_len, f64) = @splat(0.5);
495 const neg_one: @Vector(vec_len, f64) = @splat(-1.0);
496 const two: @Vector(vec_len, f64) = @splat(2.0);
497 const c0: @Vector(vec_len, f64) = @splat(0x1.555555555554ep-3);
498 const c1: @Vector(vec_len, f64) = @splat(0x1.3333333337233p-4);
499 const c2: @Vector(vec_len, f64) = @splat(0x1.6db6db67f6d9fp-5);
500 const c3: @Vector(vec_len, f64) = @splat(0x1.f1c71fbd29fbbp-6);
501 const c4: @Vector(vec_len, f64) = @splat(0x1.6e8b264d467d6p-6);
502 const c5: @Vector(vec_len, f64) = @splat(0x1.1c5997c357e9dp-6);
503 const c6: @Vector(vec_len, f64) = @splat(0x1.c86a22cd9389dp-7);
504 const c7: @Vector(vec_len, f64) = @splat(0x1.856073c22ebbep-7);
505 const c8: @Vector(vec_len, f64) = @splat(0x1.fd1151acb6bedp-8);
506 const c9: @Vector(vec_len, f64) = @splat(0x1.087182f799c1dp-6);
507 const c10: @Vector(vec_len, f64) = @splat(-0x1.6602748120927p-7);
508 const c11: @Vector(vec_len, f64) = @splat(0x1.cfa0dd1f9478p-6);
509
510 const ax = @abs(x);
511 const ax_lt_half = ax < half;
512 const is_neg = x < zero;
513 const z2 = @select(f64, ax_lt_half, x * x, @mulAdd(@Vector(vec_len, f64), -half, ax, half));
514 const z = @select(f64, ax_lt_half, ax, @sqrt(z2));
515 const z3 = z2 * z;
516 const z4 = z2 * z2;
517 const z8 = z4 * z4;
518 const p0_1 = @mulAdd(@Vector(vec_len, f64), z2, c1, c0);
519 const p2_3 = @mulAdd(@Vector(vec_len, f64), z2, c3, c2);
520 const p0_3 = @mulAdd(@Vector(vec_len, f64), z4, p2_3, p0_1);
521 const p4_5 = @mulAdd(@Vector(vec_len, f64), z2, c5, c4);
522 const p6_7 = @mulAdd(@Vector(vec_len, f64), z2, c7, c6);
523 const p4_7 = @mulAdd(@Vector(vec_len, f64), z4, p6_7, p4_5);
524 const p8_9 = @mulAdd(@Vector(vec_len, f64), z2, c9, c8);
525 const p10_11 = @mulAdd(@Vector(vec_len, f64), z2, c11, c10);
526 const p8_11 = @mulAdd(@Vector(vec_len, f64), z4, p10_11, p8_9);
527 const p4_11 = @mulAdd(@Vector(vec_len, f64), z8, p8_11, p4_7);
528 const p0_11 = @mulAdd(@Vector(vec_len, f64), z8, p4_11, p0_3);
529 const p = @mulAdd(@Vector(vec_len, f64), z3, p0_11, z);
530 const mul = @select(f64, ax_lt_half, neg_one, two);
531 const add = @select(f64, ax_lt_half, pi_over_2, @select(f64, is_neg, pi, zero));
532 return @mulAdd(@Vector(vec_len, f64), mul, @select(f64, is_neg, -p, p), add);
533}
534
535test "acosBinary32Vec.special" {
536 const input: @Vector(8, f32) = .{
537 0x0p+0,
538 -0x1p+0,
539 0x1p+0,
540 0x1.000002p+0,
541 -0x1.000002p+0,
542 math.inf(f32),
543 -math.inf(f32),
544 math.nan(f32),
545 };
546 const output = acosBinary32Vec(8, input);
547 try testing.expectApproxEqAbs(0x1.921fb6p+0, output[0], math.floatEpsAt(f32, 0x1.921fb6p+0));
548 try testing.expectApproxEqAbs(0x1.921fb6p+1, output[1], math.floatEpsAt(f32, 0x1.921fb6p+1));
549 try testing.expectEqual(0x0p+0, output[2]);
550 try testing.expect(math.isNan(output[3]));
551 try testing.expect(math.isNan(output[4]));
552 try testing.expect(math.isNan(output[5]));
553 try testing.expect(math.isNan(output[6]));
554 try testing.expect(math.isNan(output[7]));
555}
556
557test "acosBinary32Vec" {
558 const input: @Vector(10, f32) = .{
559 -0x1.13284cp-2,
560 0x1.6ca8ep-1,
561 0x1.c2ca6p-1,
562 -0x1.55f12p-1,
563 -0x1.15679ep-2,
564 -0x1.41e132p-5,
565 0x1.281b0ep-1,
566 0x1.b5ce34p-1,
567 -0x1.583482p-3,
568 -0x1.ea8224p-1,
569 };
570 const output = acosBinary32Vec(10, input);
571 try testing.expectApproxEqAbs(0x1.d7c4e6p+0, output[0], math.floatEpsAt(f32, 0x1.d7c4e6p+0));
572 try testing.expectApproxEqAbs(0x1.8e6756p-1, output[1], math.floatEpsAt(f32, 0x1.8e6756p-1));
573 try testing.expectApproxEqAbs(0x1.f9d74cp-2, output[2], math.floatEpsAt(f32, 0x1.f9d74cp-2));
574 try testing.expectApproxEqAbs(0x1.26abdcp+1, output[3], math.floatEpsAt(f32, 0x1.26abdcp+1));
575 try testing.expectApproxEqAbs(0x1.d85a44p+0, output[4], math.floatEpsAt(f32, 0x1.d85a44p+0));
576 try testing.expectApproxEqAbs(0x1.9c2f68p+0, output[5], math.floatEpsAt(f32, 0x1.9c2f68p+0));
577 try testing.expectApproxEqAbs(0x1.e881bp-1, output[6], math.floatEpsAt(f32, 0x1.e881bp-1));
578 try testing.expectApproxEqAbs(0x1.1713f6p-1, output[7], math.floatEpsAt(f32, 0x1.1713f6p-1));
579 try testing.expectApproxEqAbs(0x1.bd5accp+0, output[8], math.floatEpsAt(f32, 0x1.bd5accp+0));
580 try testing.expectApproxEqAbs(0x1.6ce7d8p+1, output[9], math.floatEpsAt(f32, 0x1.6ce7d8p+1));
581}
582
583test "acosBinary64Vec.special" {
584 const input: @Vector(8, f64) = .{
585 0x0p+0,
586 -0x1p+0,
587 0x1p+0,
588 0x1.0000000000001p+0,
589 -0x1.0000000000001p+0,
590 math.inf(f64),
591 -math.inf(f64),
592 math.nan(f64),
593 };
594 const output = acosBinary64Vec(8, input);
595 try testing.expectApproxEqAbs(0x1.921fb54442d18p+0, output[0], math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
596 try testing.expectApproxEqAbs(0x1.921fb54442d18p+1, output[1], math.floatEpsAt(f64, 0x1.921fb54442d18p+1));
597 try testing.expectEqual(0x0p+0, output[2]);
598 try testing.expect(math.isNan(output[3]));
599 try testing.expect(math.isNan(output[4]));
600 try testing.expect(math.isNan(output[5]));
601 try testing.expect(math.isNan(output[6]));
602 try testing.expect(math.isNan(output[7]));
603}
604
605test "acosBinary64Vec" {
606 const input: @Vector(10, f64) = .{
607 -0x1.13284b2b5006dp-2,
608 0x1.6ca8dfb825911p-1,
609 0x1.c2ca609de7505p-1,
610 -0x1.55f11fba96889p-1,
611 -0x1.15679e27084ddp-2,
612 -0x1.41e131b093c41p-5,
613 0x1.281b0d18455f5p-1,
614 0x1.b5ce34a51b239p-1,
615 -0x1.583481079de4dp-3,
616 -0x1.ea8223103b871p-1,
617 };
618 const output = acosBinary64Vec(10, input);
619 try testing.expectApproxEqAbs(0x1.d7c4e61020905p+0, output[0], math.floatEpsAt(f64, 0x1.d7c4e61020905p+0));
620 try testing.expectApproxEqAbs(0x1.8e6756e27c366p-1, output[1], math.floatEpsAt(f64, 0x1.8e6756e27c366p-1));
621 try testing.expectApproxEqAbs(0x1.f9d748eaf956p-2, output[2], math.floatEpsAt(f64, 0x1.f9d748eaf956p-2));
622 try testing.expectApproxEqAbs(0x1.26abdc68d07aap+1, output[3], math.floatEpsAt(f64, 0x1.26abdc68d07aap+1));
623 try testing.expectApproxEqAbs(0x1.d85a44ea44fe4p+0, output[4], math.floatEpsAt(f64, 0x1.d85a44ea44fe4p+0));
624 try testing.expectApproxEqAbs(0x1.9c2f688eee8abp+0, output[5], math.floatEpsAt(f64, 0x1.9c2f688eee8abp+0));
625 try testing.expectApproxEqAbs(0x1.e881b1d4eb2a1p-1, output[6], math.floatEpsAt(f64, 0x1.e881b1d4eb2a1p-1));
626 try testing.expectApproxEqAbs(0x1.1713f567a87efp-1, output[7], math.floatEpsAt(f64, 0x1.1713f567a87efp-1));
627 try testing.expectApproxEqAbs(0x1.bd5acbe8fcc59p+0, output[8], math.floatEpsAt(f64, 0x1.bd5acbe8fcc59p+0));
628 try testing.expectApproxEqAbs(0x1.6ce7d66f628e5p+1, output[9], math.floatEpsAt(f64, 0x1.6ce7d66f628e5p+1));
629}