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/asinf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/asin.c
6// https://git.musl-libc.org/cgit/musl/tree/src/math/asinl.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/asinf.c
12// https://github.com/ARM-software/optimized-routines/blob/master/math/aarch64/advsimd/asin.c
13
14const std = @import("../std.zig");
15const math = std.math;
16const mem = std.mem;
17const testing = std.testing;
18const builtin = @import("builtin");
19const native_endian = builtin.cpu.arch.endian();
20
21/// Returns the arc-sin of x.
22///
23/// Special Cases:
24/// - asin(+-0) = +-0
25/// - asin(x) = nan if x < -1 or x > 1
26pub fn asin(x: anytype) @TypeOf(x) {
27 const T = @TypeOf(x);
28 switch (@typeInfo(T)) {
29 .float => |info| switch (info.bits) {
30 16 => return asinBinary16(x),
31 32 => return asinBinary32(x),
32 64 => return asinBinary64(x),
33 80 => return asinExtended80(x),
34 128 => return asinBinary128(x),
35 else => comptime unreachable,
36 },
37 .vector => |info| switch (info.child) {
38 f32 => return asinBinary32Vec(info.len, x),
39 f64 => return asinBinary64Vec(info.len, x),
40 else => @compileError("unimplemented"),
41 },
42 else => comptime unreachable,
43 }
44}
45
46fn approxBinary16(z: f32) f32 {
47 const S0: f32 = 1.0000001e0;
48 const S1: f32 = 1.6664918e-1;
49 const S2: f32 = 7.55022e-2;
50 const S3: f32 = 3.9513987e-2;
51 const S4: f32 = 5.0883885e-2;
52 return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4)));
53}
54
55fn asinBinary16(x: f16) f16 {
56 const pio2: f32 = math.pi / 2.0;
57
58 const hx: u16 = @bitCast(x);
59 const ix = hx & 0x7fff;
60
61 // |x| >= 1
62 if (ix >= 0x3c00) {
63 // |x| == 1
64 if (ix == 0x3c00) {
65 // asin(+-1) = +-pi/2 with inexact
66 return @floatCast(x * pio2 + 0x1.0p-120);
67 }
68 // asin(|x| > 1) is nan
69 return 0.0 / (x - x);
70 }
71
72 // |x| < 0.5
73 if (ix < 0x3800) {
74 return @floatCast(x * approxBinary16(x * x));
75 }
76
77 // 1 > |x| >= 0.5
78 const z = (1.0 - @abs(x)) * 0.5;
79 const s = @sqrt(z);
80 const x_local = pio2 - 2.0 * s * approxBinary16(z);
81 if (hx >> 15 != 0) {
82 return @floatCast(-x_local);
83 }
84 return @floatCast(x_local);
85}
86
87fn rationalApproxBinary32(z: f32) f32 {
88 const pS0: f32 = 1.6666586697e-01;
89 const pS1: f32 = -4.2743422091e-02;
90 const pS2: f32 = -8.6563630030e-03;
91 const qS1: f32 = -7.0662963390e-01;
92
93 const p = z * (pS0 + z * (pS1 + z * pS2));
94 const q = 1.0 + z * qS1;
95 return p / q;
96}
97
98fn asinBinary32(x: f32) f32 {
99 const pio2: f64 = 1.570796326794896558e+00;
100
101 const hx: u32 = @bitCast(x);
102 const ix = hx & 0x7fff_ffff;
103
104 // |x| >= 1
105 if (ix >= 0x3f80_0000) {
106 // |x| == 1
107 if (ix == 0x3f80_0000) {
108 // asin(+-1) = +-pi/2 with inexact
109 return @floatCast(@as(f64, @floatCast(x)) * pio2 + 0x1.0p-120);
110 }
111 // asin(|x| > 1) is nan
112 return 0.0 / (x - x);
113 }
114
115 // |x| < 0.5
116 if (ix < 0x3f00_0000) {
117 // 0x1p-126 <= |x| < 0x1p-12
118 if (ix < 0x3980_0000 and ix >= 0x0080_0000) {
119 return x;
120 }
121 return x + x * rationalApproxBinary32(x * x);
122 }
123
124 // 1 > |x| >= 0.5
125 const z = (1.0 - @abs(x)) * 0.5;
126 const s: f64 = @floatCast(@sqrt(z));
127 const x_local: f32 = @floatCast(pio2 - 2.0 * (s + s * @as(f64, @floatCast(rationalApproxBinary32(z)))));
128 return if (hx >> 31 != 0) -x_local else x_local;
129}
130
131fn rationalApproxBinary64(z: f64) f64 {
132 const pS0: f64 = 1.66666666666666657415e-01;
133 const pS1: f64 = -3.25565818622400915405e-01;
134 const pS2: f64 = 2.01212532134862925881e-01;
135 const pS3: f64 = -4.00555345006794114027e-02;
136 const pS4: f64 = 7.91534994289814532176e-04;
137 const pS5: f64 = 3.47933107596021167570e-05;
138 const qS1: f64 = -2.40339491173441421878e+00;
139 const qS2: f64 = 2.02094576023350569471e+00;
140 const qS3: f64 = -6.88283971605453293030e-01;
141 const qS4: f64 = 7.70381505559019352791e-02;
142
143 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
144 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
145 return p / q;
146}
147
148fn asinBinary64(x: f64) f64 {
149 const pio2_hi: f64 = 1.57079632679489655800e+00;
150 const pio2_lo: f64 = 6.12323399573676603587e-17;
151
152 const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32);
153 const ix = hx & 0x7fffffff;
154
155 // |x| >= 1 or nan
156 if (ix >= 0x3ff0_0000) {
157 const lx: u32 = @truncate(@as(u64, @bitCast(x)));
158 // asin(1) = +-pi/2 with inexact
159 if ((ix - 0x3ff0_0000 | lx) == 0) {
160 return x * pio2_hi + 0x1.0p-120;
161 }
162 return 0.0 / (x - x);
163 }
164
165 // |x| < 0.5
166 if (ix < 0x3fe0_0000) {
167 // if 0x1p-1022 <= |x| < 0x1p-26 avoid raising overflow
168 if (ix < 0x3e50_0000 and ix >= 0x0010_0000) {
169 return x;
170 }
171 return x + x * rationalApproxBinary64(x * x);
172 }
173
174 // 1 > |x| >= 0.5
175 const z = (1.0 - @abs(x)) * 0.5;
176 const s = @sqrt(z);
177 const r = rationalApproxBinary64(z);
178 // |x| > 0.975
179 if (ix >= 0x3fef_3333) {
180 const x_local = pio2_hi - (2 * (s + s * r) - pio2_lo);
181 return if (hx >> 31 != 0) -x_local else x_local;
182 }
183 // f+c = sqrt(z)
184 const hs: u64 = @bitCast(s);
185 const f: f64 = @bitCast(hs & 0xffff_ffff_0000_0000);
186 const c: f64 = (z - f * f) / (s + f);
187 const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
188 return if (hx >> 31 != 0) -x_local else x_local;
189}
190
191fn rationalApproxExtended80(z: f80) f80 {
192 const pS0: f80 = 1.66666666666666666631e-01;
193 const pS1: f80 = -4.16313987993683104320e-01;
194 const pS2: f80 = 3.69068046323246813704e-01;
195 const pS3: f80 = -1.36213932016738603108e-01;
196 const pS4: f80 = 1.78324189708471965733e-02;
197 const pS5: f80 = -2.19216428382605211588e-04;
198 const pS6: f80 = -7.10526623669075243183e-06;
199 const qS1: f80 = -2.94788392796209867269e+00;
200 const qS2: f80 = 3.27309890266528636716e+00;
201 const qS3: f80 = -1.68285799854822427013e+00;
202 const qS4: f80 = 3.90699412641738801874e-01;
203 const qS5: f80 = -3.14365703596053263322e-02;
204
205 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6))))));
206 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5))));
207 return p / q;
208}
209
210fn asinExtended80(x: f80) f80 {
211 const pio2_hi: f80 = 1.57079632679489661926;
212 const pio2_lo: f80 = -2.50827880633416601173e-20;
213
214 const hx: u80 = @bitCast(x);
215 const se: u16 = @truncate(hx >> 64);
216 const e = se & 0x7fff;
217 const sign = se >> 15 != 0;
218
219 // |x| >= 1 or nan
220 if (e >= 0x3fff) {
221 // asin(+-1)=+-pi/2 with inexact
222 if (x == 1.0 or x == -1.0) {
223 return x * pio2_hi + 0x1p-120;
224 }
225 return 0.0 / (x - x);
226 }
227
228 // |x| < 0.5
229 if (e < 0x3fff - 1) {
230 if (e < 0x3fff - (math.floatMantissaBits(f80) + 1) / 2) {
231 // return x with inexact if x!=0
232 mem.doNotOptimizeAway(x + 0x1p120);
233 return x;
234 }
235 return x + x * rationalApproxExtended80(x * x);
236 }
237
238 // 1 > |x| >= 0.5
239 const z = (1.0 - @abs(x)) * 0.5;
240 const s = @sqrt(z);
241 const r = rationalApproxExtended80(z);
242
243 const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff);
244 if ((m >> 56) >= 0xf7) {
245 const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
246 return if (sign) -x_local else x_local;
247 }
248
249 const hs: u80 = @bitCast(s);
250 const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_0000_0000);
251 const c = (z - f * f) / (s + f);
252 const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
253 return if (sign) -x_local else x_local;
254}
255
256fn rationalApproxBinary128(z: f128) f128 {
257 const pS0: f128 = 1.66666666666666666666666666666700314e-01;
258 const pS1: f128 = -7.32816946414566252574527475428622708e-01;
259 const pS2: f128 = 1.34215708714992334609030036562143589e+00;
260 const pS3: f128 = -1.32483151677116409805070261790752040e+00;
261 const pS4: f128 = 7.61206183613632558824485341162121989e-01;
262 const pS5: f128 = -2.56165783329023486777386833928147375e-01;
263 const pS6: f128 = 4.80718586374448793411019434585413855e-02;
264 const pS7: f128 = -4.42523267167024279410230886239774718e-03;
265 const pS8: f128 = 1.44551535183911458253205638280410064e-04;
266 const pS9: f128 = -2.10558957916600254061591040482706179e-07;
267 const qS1: f128 = -4.84690167848739751544716485245697428e+00;
268 const qS2: f128 = 9.96619113536172610135016921140206980e+00;
269 const qS3: f128 = -1.13177895428973036660836798461641458e+01;
270 const qS4: f128 = 7.74004374389488266169304117714658761e+00;
271 const qS5: f128 = -3.25871986053534084709023539900339905e+00;
272 const qS6: f128 = 8.27830318881232209752469022352928864e-01;
273 const qS7: f128 = -1.18768052702942805423330715206348004e-01;
274 const qS8: f128 = 8.32600764660522313269101537926539470e-03;
275 const qS9: f128 = -1.99407384882605586705979504567947007e-04;
276
277 const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9)))))))));
278 const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9))))))));
279 return p / q;
280}
281
282fn asinBinary128(x: f128) f128 {
283 const pio2_hi: f128 = 1.57079632679489661923132169163975140;
284 const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35;
285
286 const hx: u128 = @bitCast(x);
287 const se: u16 = @truncate(hx >> 112);
288 const e = se & 0x7fff;
289 const sign = se >> 15 != 0;
290
291 // |x| >= 1 or nan
292 if (e >= 0x3fff) {
293 // asin(+-1)=+-pi/2 with inexact
294 if (x == 1.0 or x == -1.0) {
295 return x * pio2_hi + 0x1p-120;
296 }
297 return 0.0 / (x - x);
298 }
299
300 // |x| < 0.5
301 if (e < 0x3fff - 1) {
302 if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) {
303 // return x with inexact if x!=0
304 mem.doNotOptimizeAway(x + 0x1p120);
305 return x;
306 }
307 return x + x * rationalApproxBinary128(x * x);
308 }
309
310 // 1 > |x| >= 0.5
311 const z = (1.0 - @abs(x)) * 0.5;
312 const s = @sqrt(z);
313 const r = rationalApproxBinary128(z);
314
315 const top: u16 = @truncate((hx >> 96) & 0x0000_ffff);
316 if (top >= 0xee00) {
317 const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
318 return if (sign) -x_local else x_local;
319 }
320
321 const hs: u128 = @bitCast(s);
322 const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000);
323 const c = (z - f * f) / (s + f);
324 const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
325 return if (sign) -x_local else x_local;
326}
327
328test "asinBinary16.special" {
329 try testing.expectApproxEqAbs(0x1.92p0, asinBinary16(0x1p+0), math.floatEpsAt(f16, 0x1.92p0));
330 try testing.expectApproxEqAbs(-0x1.92p0, asinBinary16(-0x1p+0), math.floatEpsAt(f16, -0x1.92p0));
331 try testing.expectEqual(0x0p+0, asinBinary16(0x0p+0));
332 try testing.expectEqual(0x0p+0, asinBinary16(-0x0p+0));
333 try testing.expect(math.isNan(asinBinary16(0x1.004p0)));
334 try testing.expect(math.isNan(asinBinary16(-0x1.004p0)));
335 try testing.expect(math.isNan(asinBinary16(math.inf(f16))));
336 try testing.expect(math.isNan(asinBinary16(-math.inf(f16))));
337 try testing.expect(math.isNan(asinBinary16(math.nan(f16))));
338}
339
340test "asinBinary16" {
341 try testing.expectApproxEqAbs(-0x1.e4cp-6, asinBinary16(-0x1.e4cp-6), math.floatEpsAt(f16, -0x1.e4cp-6));
342 try testing.expectApproxEqAbs(0x1.2a8p0, asinBinary16(0x1.d68p-1), math.floatEpsAt(f16, 0x1.2a8p0));
343 try testing.expectApproxEqAbs(-0x1.eep-1, asinBinary16(-0x1.a4cp-1), math.floatEpsAt(f16, -0x1.eep-1));
344 try testing.expectApproxEqAbs(-0x1.0d4p-2, asinBinary16(-0x1.0a4p-2), math.floatEpsAt(f16, -0x1.0d4p-2));
345 try testing.expectApproxEqAbs(0x1.3c8p-1, asinBinary16(0x1.28cp-1), math.floatEpsAt(f16, 0x1.3c8p-1));
346 try testing.expectApproxEqAbs(0x1.298p-3, asinBinary16(0x1.284p-3), math.floatEpsAt(f16, 0x1.298p-3));
347 try testing.expectApproxEqAbs(-0x1.784p-1, asinBinary16(-0x1.574p-1), math.floatEpsAt(f16, -0x1.784p-1));
348 try testing.expectApproxEqAbs(-0x1.6a4p-1, asinBinary16(-0x1.4ccp-1), math.floatEpsAt(f16, -0x1.6a4p-1));
349 try testing.expectApproxEqAbs(0x1.e84p-1, asinBinary16(0x1.a18p-1), math.floatEpsAt(f16, 0x1.e84p-1));
350 try testing.expectApproxEqAbs(0x1.83cp-2, asinBinary16(0x1.7a8p-2), math.floatEpsAt(f16, 0x1.83cp-2));
351}
352
353test "asinBinary32.special" {
354 try testing.expectApproxEqAbs(0x1.921fb6p+0, asinBinary32(0x1p+0), math.floatEpsAt(f32, 0x1.921fb6p+0));
355 try testing.expectApproxEqAbs(-0x1.921fb6p+0, asinBinary32(-0x1p+0), math.floatEpsAt(f32, -0x1.921fb6p+0));
356 try testing.expectEqual(0x0p+0, asinBinary32(0x0p+0));
357 try testing.expectEqual(0x0p+0, asinBinary32(-0x0p+0));
358 try testing.expect(math.isNan(asinBinary32(0x1.000002p+0)));
359 try testing.expect(math.isNan(asinBinary32(-0x1.000002p+0)));
360 try testing.expect(math.isNan(asinBinary32(math.inf(f32))));
361 try testing.expect(math.isNan(asinBinary32(-math.inf(f32))));
362 try testing.expect(math.isNan(asinBinary32(math.nan(f32))));
363}
364
365test "asinBinary32" {
366 try testing.expectApproxEqAbs(-0x1.4c868p-4, asinBinary32(-0x1.4c2906p-4), math.floatEpsAt(f32, -0x1.4c868p-4));
367 try testing.expectApproxEqAbs(0x1.130648p-1, asinBinary32(0x1.05fcfap-1), math.floatEpsAt(f32, 0x1.130648p-1));
368 try testing.expectApproxEqAbs(0x1.090abcp-1, asinBinary32(0x1.fab976p-2), math.floatEpsAt(f32, 0x1.090abcp-1));
369 try testing.expectApproxEqAbs(0x1.c39fa2p-1, asinBinary32(0x1.8b4b8cp-1), math.floatEpsAt(f32, 0x1.c39fa2p-1));
370 try testing.expectApproxEqAbs(0x1.9c332p-1, asinBinary32(0x1.7117c2p-1), math.floatEpsAt(f32, 0x1.9c332p-1));
371 try testing.expectApproxEqAbs(0x1.e62a1cp-5, asinBinary32(0x1.e5e112p-5), math.floatEpsAt(f32, 0x1.e62a1cp-5));
372 try testing.expectApproxEqAbs(-0x1.0a65dep-2, asinBinary32(-0x1.07673p-2), math.floatEpsAt(f32, -0x1.0a65dep-2));
373 try testing.expectApproxEqAbs(-0x1.25046p-2, asinBinary32(-0x1.2108dep-2), math.floatEpsAt(f32, -0x1.25046p-2));
374 try testing.expectApproxEqAbs(-0x1.6c6f0cp-1, asinBinary32(-0x1.4e6e6cp-1), math.floatEpsAt(f32, -0x1.6c6f0cp-1));
375 try testing.expectApproxEqAbs(0x1.350f7ap-1, asinBinary32(0x1.22a16ap-1), math.floatEpsAt(f32, 0x1.350f7ap-1));
376}
377
378test "asinBinary64.special" {
379 try testing.expectApproxEqAbs(0x1.921fb54442d18p+0, asinBinary64(0x1p+0), math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
380 try testing.expectApproxEqAbs(-0x1.921fb54442d18p+0, asinBinary64(-0x1p+0), math.floatEpsAt(f64, -0x1.921fb54442d18p+0));
381 try testing.expectEqual(0x0p+0, asinBinary64(0x0p+0));
382 try testing.expectEqual(0x0p+0, asinBinary64(-0x0p+0));
383 try testing.expect(math.isNan(asinBinary64(0x1.000002p+0)));
384 try testing.expect(math.isNan(asinBinary64(-0x1.000002p+0)));
385 try testing.expect(math.isNan(asinBinary64(math.inf(f64))));
386 try testing.expect(math.isNan(asinBinary64(-math.inf(f64))));
387 try testing.expect(math.isNan(asinBinary64(math.nan(f64))));
388}
389
390test "asinBinary64" {
391 try testing.expectApproxEqAbs(0x1.fae86c5941692p-2, asinBinary64(0x1.e674fba3e40d5p-2), math.floatEpsAt(f64, 0x1.fae86c5941692p-2));
392 try testing.expectApproxEqAbs(-0x1.46b6ad730c93ap-1, asinBinary64(-0x1.30fd0566fd979p-1), math.floatEpsAt(f64, -0x1.46b6ad730c93ap-1));
393 try testing.expectApproxEqAbs(0x1.6be0be8074eep-2, asinBinary64(0x1.6444a25abfeaap-2), math.floatEpsAt(f64, 0x1.6be0be8074eep-2));
394 try testing.expectApproxEqAbs(0x1.5a7e98f53f717p-1, asinBinary64(0x1.40a53228d1a13p-1), math.floatEpsAt(f64, 0x1.5a7e98f53f717p-1));
395 try testing.expectApproxEqAbs(-0x1.1ea2602d14e8p0, asinBinary64(-0x1.ccc6d64845cfdp-1), math.floatEpsAt(f64, -0x1.1ea2602d14e8p0));
396 try testing.expectApproxEqAbs(-0x1.d2c2634193158p-1, asinBinary64(-0x1.94bd91b7fc74bp-1), math.floatEpsAt(f64, -0x1.d2c2634193158p-1));
397 try testing.expectApproxEqAbs(-0x1.982d5f1895d2p-2, asinBinary64(-0x1.8d741b5797fccp-2), math.floatEpsAt(f64, -0x1.982d5f1895d2p-2));
398 try testing.expectApproxEqAbs(-0x1.3fdaf7dfdc864p-3, asinBinary64(-0x1.3e8e7e15881c5p-3), math.floatEpsAt(f64, -0x1.3fdaf7dfdc864p-3));
399 try testing.expectApproxEqAbs(-0x1.9269540735b7bp-2, asinBinary64(-0x1.88222d8ab8ca9p-2), math.floatEpsAt(f64, -0x1.9269540735b7bp-2));
400 try testing.expectApproxEqAbs(-0x1.474c4c6625527p-2, asinBinary64(-0x1.41c0e9babcbd2p-2), math.floatEpsAt(f64, -0x1.474c4c6625527p-2));
401}
402
403test "asinExtended80.special" {
404 try testing.expectApproxEqAbs(0x1.921fb54442d1846ap+0, asinExtended80(0x1p+0), math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0));
405 try testing.expectApproxEqAbs(-0x1.921fb54442d1846ap+0, asinExtended80(-0x1p+0), math.floatEpsAt(f80, -0x1.921fb54442d1846ap+0));
406 try testing.expectEqual(0x0p+0, asinExtended80(0x0p+0));
407 try testing.expectEqual(0x0p+0, asinExtended80(-0x0p+0));
408 try testing.expect(math.isNan(asinExtended80(0x1.0000000000000002p+0)));
409 try testing.expect(math.isNan(asinExtended80(-0x1.0000000000000002p+0)));
410 try testing.expect(math.isNan(asinExtended80(math.inf(f80))));
411 try testing.expect(math.isNan(asinExtended80(-math.inf(f80))));
412 try testing.expect(math.isNan(asinExtended80(math.nan(f80))));
413}
414
415test "asinExtended80" {
416 try testing.expectApproxEqAbs(0x1.63cfb560149daa9p-9, asinExtended80(0x1.63cf98bc52ce0da8p-9), math.floatEpsAt(f80, 0x1.63cfb560149daa9p-9));
417 try testing.expectApproxEqAbs(-0x1.113cbacd8cd1b96cp-1, asinExtended80(-0x1.0473756f7ae930dp-1), math.floatEpsAt(f80, -0x1.113cbacd8cd1b96cp-1));
418 try testing.expectApproxEqAbs(-0x1.2721b231d197b064p-2, asinExtended80(-0x1.2310057e005cc288p-2), math.floatEpsAt(f80, -0x1.2721b231d197b064p-2));
419 try testing.expectApproxEqAbs(0x1.547c408c5d2b05aap0, asinExtended80(0x1.f13b03bd685d96eap-1), math.floatEpsAt(f80, 0x1.547c408c5d2b05aap0));
420 try testing.expectApproxEqAbs(-0x1.296b76bfadbb5cecp0, asinExtended80(-0x1.d5c507e3ef84041cp-1), math.floatEpsAt(f80, -0x1.296b76bfadbb5cecp0));
421 try testing.expectApproxEqAbs(0x1.b572da8729a84f2ap-1, asinExtended80(0x1.8222cbc9147153d8p-1), math.floatEpsAt(f80, 0x1.b572da8729a84f2ap-1));
422 try testing.expectApproxEqAbs(-0x1.42c9e80ac0524dap-11, asinExtended80(-0x1.42c9e6b4a088a246p-11), math.floatEpsAt(f80, -0x1.42c9e80ac0524dap-11));
423 try testing.expectApproxEqAbs(-0x1.920ca86aef6c3028p-3, asinExtended80(-0x1.8f78d49deadb521cp-3), math.floatEpsAt(f80, -0x1.920ca86aef6c3028p-3));
424 try testing.expectApproxEqAbs(-0x1.b91cb4f7204d92fp-2, asinExtended80(-0x1.ab98792783515774p-2), math.floatEpsAt(f80, -0x1.b91cb4f7204d92fp-2));
425 try testing.expectApproxEqAbs(-0x1.1f20815fdc4c5304p-1, asinExtended80(-0x1.104fe30cef6800aap-1), math.floatEpsAt(f80, -0x1.1f20815fdc4c5304p-1));
426}
427
428test "asinBinary128.special" {
429 try testing.expectApproxEqAbs(0x1.921fb54442d18469898cc51701b8p0, asinBinary128(0x1p+0), math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0));
430 try testing.expectApproxEqAbs(-0x1.921fb54442d18469898cc51701b8p0, asinBinary128(-0x1p+0), math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p0));
431 try testing.expectEqual(0x0p+0, asinBinary128(0x0p+0));
432 try testing.expectEqual(0x0p+0, asinBinary128(-0x0p+0));
433 try testing.expect(math.isNan(asinBinary128(0x1.0000000000000000000000000001p0)));
434 try testing.expect(math.isNan(asinBinary128(-0x1.0000000000000000000000000001p0)));
435 try testing.expect(math.isNan(asinBinary128(math.inf(f128))));
436 try testing.expect(math.isNan(asinBinary128(-math.inf(f128))));
437 try testing.expect(math.isNan(asinBinary128(math.nan(f128))));
438}
439
440test "asinBinary128" {
441 if (builtin.cpu.arch.isSPARC()) return error.SkipZigTest;
442
443 try testing.expectApproxEqAbs(0x1.87e9c740d7837f8e8fa667988fbep-3, asinBinary128(0x1.85868ce287ca0196b01c25fec5ffp-3), math.floatEpsAt(f128, 0x1.87e9c740d7837f8e8fa667988fbep-3));
444 try testing.expectApproxEqAbs(0x1.bd11a474e864213b48e0f005f1f4p-1, asinBinary128(0x1.8718d6d30b4daed08d04ef59f478p-1), math.floatEpsAt(f128, 0x1.bd11a474e864213b48e0f005f1f4p-1));
445 try testing.expectApproxEqAbs(0x1.20b56f8b42649fe72d1f8d68a378p-1, asinBinary128(0x1.11a67640cd7f0ba5d5e362f3abfap-1), math.floatEpsAt(f128, 0x1.20b56f8b42649fe72d1f8d68a378p-1));
446 try testing.expectApproxEqAbs(-0x1.0dc3a7ddb9736e5ad699bf338566p0, asinBinary128(-0x1.bd13bf14a9dce22188e52650daa7p-1), math.floatEpsAt(f128, -0x1.0dc3a7ddb9736e5ad699bf338566p0));
447 try testing.expectApproxEqAbs(-0x1.f250716038f70fa50a5826c03802p-2, asinBinary128(-0x1.dee0bc217fc462af57c484eefa71p-2), math.floatEpsAt(f128, -0x1.f250716038f70fa50a5826c03802p-2));
448 try testing.expectApproxEqAbs(-0x1.47a8b4cdd327f90056722feddbabp0, asinBinary128(-0x1.ea7df9139371c10b9d6fd2bbccd3p-1), math.floatEpsAt(f128, -0x1.47a8b4cdd327f90056722feddbabp0));
449 try testing.expectApproxEqAbs(0x1.079178d52be662dec67e2cd7f6e9p-2, asinBinary128(0x1.04aaea6de3b5a616460702f26dfcp-2), math.floatEpsAt(f128, 0x1.079178d52be662dec67e2cd7f6e9p-2));
450 try testing.expectApproxEqAbs(-0x1.192df5a8d71702cf1e27014887b2p0, asinBinary128(-0x1.c7ea85e6b61be666435a7d99444cp-1), math.floatEpsAt(f128, -0x1.192df5a8d71702cf1e27014887b2p0));
451 try testing.expectApproxEqAbs(-0x1.97f1092fd94ac0fdfddae2e1222bp-1, asinBinary128(-0x1.6e210214e40edf6c8479998189d1p-1), math.floatEpsAt(f128, -0x1.97f1092fd94ac0fdfddae2e1222bp-1));
452 try testing.expectApproxEqAbs(-0x1.97b62bc5ae6512093828828325e1p-3, asinBinary128(-0x1.95061bf93ed6986a45d20f0e1064p-3), math.floatEpsAt(f128, -0x1.97b62bc5ae6512093828828325e1p-3));
453}
454
455fn asinBinary32Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f32)) @TypeOf(x) {
456 const pi_over_2: @Vector(vec_len, f32) = @splat(math.pi / 2.0);
457 const zero: @Vector(vec_len, f32) = @splat(0.0);
458 const half: @Vector(vec_len, f32) = @splat(0.5);
459 const neg_two: @Vector(vec_len, f32) = @splat(-2.0);
460 const c0: @Vector(vec_len, f32) = @splat(0x1.55555ep-3);
461 const c1: @Vector(vec_len, f32) = @splat(0x1.33261ap-4);
462 const c2: @Vector(vec_len, f32) = @splat(0x1.70d7dcp-5);
463 const c3: @Vector(vec_len, f32) = @splat(0x1.b059dp-6);
464 const c4: @Vector(vec_len, f32) = @splat(0x1.3af7d8p-5);
465
466 const ax = @abs(x);
467 const ax_lt_half = ax < half;
468 const z2 = @select(f32, ax_lt_half, x * x, @mulAdd(@Vector(vec_len, f32), -half, ax, half));
469 const z = @select(f32, ax_lt_half, ax, @sqrt(z2));
470 const z3 = z2 * z;
471 const p3_4 = @mulAdd(@Vector(vec_len, f32), z2, c4, c3);
472 const p2_4 = @mulAdd(@Vector(vec_len, f32), z2, p3_4, c2);
473 const p1_4 = @mulAdd(@Vector(vec_len, f32), z2, p2_4, c1);
474 const p0_4 = @mulAdd(@Vector(vec_len, f32), z2, p1_4, c0);
475 const p = @mulAdd(@Vector(vec_len, f32), z3, p0_4, z);
476 const y = @select(f32, ax_lt_half, p, @mulAdd(@Vector(vec_len, f32), p, neg_two, pi_over_2));
477 return @select(f32, x < zero, -y, y);
478}
479
480fn asinBinary64Vec(comptime vec_len: comptime_int, x: @Vector(vec_len, f64)) @TypeOf(x) {
481 const pi_over_2: @Vector(vec_len, f64) = @splat(math.pi / 2.0);
482 const zero: @Vector(vec_len, f64) = @splat(0.0);
483 const half: @Vector(vec_len, f64) = @splat(0.5);
484 const neg_two: @Vector(vec_len, f64) = @splat(-2.0);
485 const c0: @Vector(vec_len, f64) = @splat(0x1.555555555554ep-3);
486 const c1: @Vector(vec_len, f64) = @splat(0x1.3333333337233p-4);
487 const c2: @Vector(vec_len, f64) = @splat(0x1.6db6db67f6d9fp-5);
488 const c3: @Vector(vec_len, f64) = @splat(0x1.f1c71fbd29fbbp-6);
489 const c4: @Vector(vec_len, f64) = @splat(0x1.6e8b264d467d6p-6);
490 const c5: @Vector(vec_len, f64) = @splat(0x1.1c5997c357e9dp-6);
491 const c6: @Vector(vec_len, f64) = @splat(0x1.c86a22cd9389dp-7);
492 const c7: @Vector(vec_len, f64) = @splat(0x1.856073c22ebbep-7);
493 const c8: @Vector(vec_len, f64) = @splat(0x1.fd1151acb6bedp-8);
494 const c9: @Vector(vec_len, f64) = @splat(0x1.087182f799c1dp-6);
495 const c10: @Vector(vec_len, f64) = @splat(-0x1.6602748120927p-7);
496 const c11: @Vector(vec_len, f64) = @splat(0x1.cfa0dd1f9478p-6);
497
498 const ax = @abs(x);
499 const ax_lt_half = ax < half;
500 const z2 = @select(f64, ax_lt_half, x * x, @mulAdd(@Vector(vec_len, f64), -half, ax, half));
501 const z = @select(f64, ax_lt_half, ax, @sqrt(z2));
502 const z3 = z2 * z;
503 const z4 = z2 * z2;
504 const z8 = z4 * z4;
505 const p0_1 = @mulAdd(@Vector(vec_len, f64), z2, c1, c0);
506 const p2_3 = @mulAdd(@Vector(vec_len, f64), z2, c3, c2);
507 const p0_3 = @mulAdd(@Vector(vec_len, f64), z4, p2_3, p0_1);
508 const p4_5 = @mulAdd(@Vector(vec_len, f64), z2, c5, c4);
509 const p6_7 = @mulAdd(@Vector(vec_len, f64), z2, c7, c6);
510 const p4_7 = @mulAdd(@Vector(vec_len, f64), z4, p6_7, p4_5);
511 const p8_9 = @mulAdd(@Vector(vec_len, f64), z2, c9, c8);
512 const p10_11 = @mulAdd(@Vector(vec_len, f64), z2, c11, c10);
513 const p8_11 = @mulAdd(@Vector(vec_len, f64), z4, p10_11, p8_9);
514 const p4_11 = @mulAdd(@Vector(vec_len, f64), z8, p8_11, p4_7);
515 const p0_11 = @mulAdd(@Vector(vec_len, f64), z8, p4_11, p0_3);
516 const p = @mulAdd(@Vector(vec_len, f64), z3, p0_11, z);
517 const y = @select(f64, ax_lt_half, p, @mulAdd(@Vector(vec_len, f64), p, neg_two, pi_over_2));
518 return @select(f64, x < zero, -y, y);
519}
520
521test "asinBinary32Vec.special" {
522 const input: @Vector(9, f32) = .{
523 0x1p+0,
524 -0x1p+0,
525 0x0p+0,
526 -0x0p+0,
527 0x1.000002p+0,
528 -0x1.000002p+0,
529 math.inf(f32),
530 -math.inf(f32),
531 math.nan(f32),
532 };
533 const output = asinBinary32Vec(9, input);
534 try testing.expectApproxEqAbs(0x1.921fb6p+0, output[0], math.floatEpsAt(f32, 0x1.921fb6p+0));
535 try testing.expectApproxEqAbs(-0x1.921fb6p+0, output[1], math.floatEpsAt(f32, -0x1.921fb6p+0));
536 try testing.expectEqual(0x0p+0, output[2]);
537 try testing.expectEqual(0x0p+0, output[3]);
538 try testing.expect(math.isNan(output[4]));
539 try testing.expect(math.isNan(output[5]));
540 try testing.expect(math.isNan(output[6]));
541 try testing.expect(math.isNan(output[7]));
542 try testing.expect(math.isNan(output[8]));
543}
544
545test "asinBinary32Vec" {
546 const input: @Vector(10, f32) = .{
547 -0x1.4c2906p-4,
548 0x1.05fcfap-1,
549 0x1.fab976p-2,
550 0x1.8b4b8cp-1,
551 0x1.7117c2p-1,
552 0x1.e5e112p-5,
553 -0x1.07673p-2,
554 -0x1.2108dep-2,
555 -0x1.4e6e6cp-1,
556 0x1.22a16ap-1,
557 };
558 const output = asinBinary32Vec(10, input);
559 try testing.expectApproxEqAbs(-0x1.4c868p-4, output[0], math.floatEpsAt(f32, -0x1.4c868p-4));
560 try testing.expectApproxEqAbs(0x1.130648p-1, output[1], math.floatEpsAt(f32, 0x1.130648p-1));
561 try testing.expectApproxEqAbs(0x1.090abcp-1, output[2], math.floatEpsAt(f32, 0x1.090abcp-1));
562 try testing.expectApproxEqAbs(0x1.c39fa2p-1, output[3], math.floatEpsAt(f32, 0x1.c39fa2p-1));
563 try testing.expectApproxEqAbs(0x1.9c332p-1, output[4], math.floatEpsAt(f32, 0x1.9c332p-1));
564 try testing.expectApproxEqAbs(0x1.e62a1cp-5, output[5], math.floatEpsAt(f32, 0x1.e62a1cp-5));
565 try testing.expectApproxEqAbs(-0x1.0a65dep-2, output[6], math.floatEpsAt(f32, -0x1.0a65dep-2));
566 try testing.expectApproxEqAbs(-0x1.25046p-2, output[7], math.floatEpsAt(f32, -0x1.25046p-2));
567 try testing.expectApproxEqAbs(-0x1.6c6f0cp-1, output[8], math.floatEpsAt(f32, -0x1.6c6f0cp-1));
568 try testing.expectApproxEqAbs(0x1.350f7ap-1, output[9], math.floatEpsAt(f32, 0x1.350f7ap-1));
569}
570
571test "asinBinary64Vec.special" {
572 const input: @Vector(9, f64) = .{
573 0x1p+0,
574 -0x1p+0,
575 0x0p+0,
576 -0x0p+0,
577 0x1.000002p+0,
578 -0x1.000002p+0,
579 math.inf(f64),
580 -math.inf(f64),
581 math.nan(f64),
582 };
583 const output = asinBinary64Vec(9, input);
584 try testing.expectApproxEqAbs(0x1.921fb54442d18p+0, output[0], math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
585 try testing.expectApproxEqAbs(-0x1.921fb54442d18p+0, output[1], math.floatEpsAt(f64, -0x1.921fb54442d18p+0));
586 try testing.expectEqual(0x0p+0, output[2]);
587 try testing.expectEqual(0x0p+0, output[3]);
588 try testing.expect(math.isNan(output[4]));
589 try testing.expect(math.isNan(output[5]));
590 try testing.expect(math.isNan(output[6]));
591 try testing.expect(math.isNan(output[7]));
592 try testing.expect(math.isNan(output[8]));
593}
594
595test "asinBinary64Vec" {
596 const input: @Vector(10, f64) = .{
597 0x1.e674fba3e40d5p-2,
598 -0x1.30fd0566fd979p-1,
599 0x1.6444a25abfeaap-2,
600 0x1.40a53228d1a13p-1,
601 -0x1.ccc6d64845cfdp-1,
602 -0x1.94bd91b7fc74bp-1,
603 -0x1.8d741b5797fccp-2,
604 -0x1.3e8e7e15881c5p-3,
605 -0x1.88222d8ab8ca9p-2,
606 -0x1.41c0e9babcbd2p-2,
607 };
608 const output = asinBinary64Vec(10, input);
609 try testing.expectApproxEqAbs(0x1.fae86c5941692p-2, output[0], math.floatEpsAt(f64, 0x1.fae86c5941692p-2));
610 try testing.expectApproxEqAbs(-0x1.46b6ad730c93ap-1, output[1], math.floatEpsAt(f64, -0x1.46b6ad730c93ap-1));
611 try testing.expectApproxEqAbs(0x1.6be0be8074eep-2, output[2], math.floatEpsAt(f64, 0x1.6be0be8074eep-2));
612 try testing.expectApproxEqAbs(0x1.5a7e98f53f717p-1, output[3], math.floatEpsAt(f64, 0x1.5a7e98f53f717p-1));
613 try testing.expectApproxEqAbs(-0x1.1ea2602d14e8p0, output[4], math.floatEpsAt(f64, -0x1.1ea2602d14e8p0));
614 try testing.expectApproxEqAbs(-0x1.d2c2634193158p-1, output[5], math.floatEpsAt(f64, -0x1.d2c2634193158p-1));
615 try testing.expectApproxEqAbs(-0x1.982d5f1895d2p-2, output[6], math.floatEpsAt(f64, -0x1.982d5f1895d2p-2));
616 try testing.expectApproxEqAbs(-0x1.3fdaf7dfdc864p-3, output[7], math.floatEpsAt(f64, -0x1.3fdaf7dfdc864p-3));
617 try testing.expectApproxEqAbs(-0x1.9269540735b7bp-2, output[8], math.floatEpsAt(f64, -0x1.9269540735b7bp-2));
618 try testing.expectApproxEqAbs(-0x1.474c4c6625527p-2, output[9], math.floatEpsAt(f64, -0x1.474c4c6625527p-2));
619}