Commit 757ec185f0eb91a15c4bdbe0201f0d998f30258c | zig

author	Zhenming-Lin <3575188313@qq.com> 2026-01-28 23:33:35+01:00
committer	Andrew Kelley <andrew@ziglang.org> 2026-01-28 23:33:35+01:00
log	757ec185f0eb91a15c4bdbe0201f0d998f30258c
tree	0334713c4c82f8f0d0879e66e02aa92d5796f56e
parent	3b10383114f7aee96f18f7b3aac7c0e08d840b42

Add `f16`, `f80` and `f128` support for `acos` and `asin` (#30997)

The software impl of `acos` and `asin` depends on the `sqrt` op. Since support for `sqrt` in `f16`, `f80`, and `f128` has been added, the impl of `acos` and `asin` for `f16`, `f80`, and `f128` is now being supplemented.

Reviewed-on: https://codeberg.org/ziglang/zig/pulls/30997
Reviewed-by: Andrew Kelley <andrew@ziglang.org>
Co-authored-by: lzm-build <3575188313@qq.com>
Co-committed-by: lzm-build <3575188313@qq.com>

2 files changed, 675 insertions(+), 168 deletions(-)

lib
- std
  - math
    acos.zig
    asin.zig

lib/std/math/acos.zig+331-78

...	...	@@ -3,10 +3,13 @@
3	3	//
4	4	// https://git.musl-libc.org/cgit/musl/tree/src/math/acosf.c
5	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
6	7
7	8	const std = @import("../std.zig");
8	9	const math = std.math;
9		const expect = std.testing.expect;
	10	const testing = std.testing;
	11	const builtin = @import("builtin");
	12	const native_endian = builtin.cpu.arch.endian();
10	13
11	14	/// Returns the arc-cosine of x.
12	15	///
...	...	@@ -15,71 +18,120 @@ const expect = std.testing.expect;
15	18	pub fn acos(x: anytype) @TypeOf(x) {
16	19	const T = @TypeOf(x);
17	20	return switch (T) {
18		f32 => acos32(x),
19		f64 => acos64(x),
	21	f16 => acosBinary16(x),
	22	f32 => acosBinary32(x),
	23	f64 => acosBinary64(x),
	24	f80 => acosExtended80(x),
	25	f128 => acosBinary128(x),
20	26	else => @compileError("acos not implemented for " ++ @typeName(T)),
21	27	};
22	28	}
23	29
24		fn r32(z: f32) f32 {
25		const pS0 = 1.6666586697e-01;
26		const pS1 = -4.2743422091e-02;
27		const pS2 = -8.6563630030e-03;
28		const qS1 = -7.0662963390e-01;
	30	fn approxBinary16(z: f32) f32 {
	31	const S0: f32 = 1.0000001e0;
	32	const S1: f32 = 1.6664918e-1;
	33	const S2: f32 = 7.55022e-2;
	34	const S3: f32 = 3.9513987e-2;
	35	const S4: f32 = 5.0883885e-2;
	36	return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4)));
	37	}
	38
	39	fn acosBinary16(x: f16) f16 {
	40	const pio2: f32 = math.pi / 2.0;
	41
	42	const hx: u16 = @bitCast(x);
	43	const ix: u16 = hx & 0x7fff;
	44
	45	// \|x\| >= 1 or nan
	46	if (ix >= 0x3c00) {
	47	if (ix == 0x3c00) {
	48	if (hx >> 15 != 0) {
	49	return @floatCast(2.0 * pio2 + 0x1p-120);
	50	}
	51	return 0.0;
	52	}
	53	return 0.0 / (x - x);
	54	}
	55
	56	const xf: f32 = @floatCast(x);
	57
	58	// \|x\| < 0.5
	59	if (ix < 0x3800) {
	60	return @floatCast(pio2 - xf * approxBinary16(xf * xf));
	61	}
	62
	63	// x < -0.5
	64	if (hx >> 15 != 0) {
	65	const z = (1.0 + xf) * 0.5;
	66	const s = @sqrt(z);
	67	const w = approxBinary16(z) * s;
	68	return @floatCast(2.0 * (pio2 - w));
	69	}
	70
	71	// x > 0.5
	72	const z = (1.0 - xf) * 0.5;
	73	const s = @sqrt(z);
	74	const w = approxBinary16(z) * s;
	75	return @floatCast(2.0 * w);
	76	}
	77
	78	fn rationalApproxBinary32(z: f32) f32 {
	79	const pS0: f32 = 1.6666586697e-01;
	80	const pS1: f32 = -4.2743422091e-02;
	81	const pS2: f32 = -8.6563630030e-03;
	82	const qS1: f32 = -7.0662963390e-01;
29	83
30	84	const p = z * (pS0 + z * (pS1 + z * pS2));
31	85	const q = 1.0 + z * qS1;
32	86	return p / q;
33	87	}
34	88
35		fn acos32(x: f32) f32 {
36		const pio2_hi = 1.5707962513e+00;
37		const pio2_lo = 7.5497894159e-08;
	89	fn acosBinary32(x: f32) f32 {
	90	const pio2_hi: f32 = 1.5707962513e+00;
	91	const pio2_lo: f32 = 7.5497894159e-08;
38	92
39		const hx: u32 = @as(u32, @bitCast(x));
40		const ix: u32 = hx & 0x7FFFFFFF;
	93	const hx: u32 = @bitCast(x);
	94	const ix: u32 = hx & 0x7fff_ffff;
41	95
42	96	// \|x\| >= 1 or nan
43		if (ix >= 0x3F800000) {
44		if (ix == 0x3F800000) {
	97	if (ix >= 0x3f800000) {
	98	if (ix == 0x3f800000) {
45	99	if (hx >> 31 != 0) {
46	100	return 2.0 * pio2_hi + 0x1.0p-120;
47		} else {
48		return 0.0;
49	101	}
50		} else {
51		return (x - x) / 0;
	102	return 0.0;
52	103	}
	104	return 0.0 / (x - x);
53	105	}
54	106
55	107	// \|x\| < 0.5
56		if (ix < 0x3F000000) {
57		if (ix <= 0x32800000) { // \|x\| < 2^(-26)
	108	if (ix < 0x3f00_0000) {
	109	// \|x\| < 2^(-26)
	110	if (ix <= 0x3280_0000) {
58	111	return pio2_hi + 0x1.0p-120;
59		} else {
60		return pio2_hi - (x - (pio2_lo - x * r32(x * x)));
61	112	}
	113	return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary32(x * x)));
62	114	}
63	115
64	116	// x < -0.5
65	117	if (hx >> 31 != 0) {
66	118	const z = (1 + x) * 0.5;
67	119	const s = @sqrt(z);
68		const w = r32(z) * s - pio2_lo;
69		return 2 * (pio2_hi - (s + w));
	120	const w = rationalApproxBinary32(z) * s - pio2_lo;
	121	return 2.0 * (pio2_hi - (s + w));
70	122	}
71	123
72	124	// x > 0.5
73	125	const z = (1.0 - x) * 0.5;
74	126	const s = @sqrt(z);
75		const jx = @as(u32, @bitCast(s));
76		const df = @as(f32, @bitCast(jx & 0xFFFFF000));
	127	const hs: u32 = @bitCast(s);
	128	const df: f32 = @bitCast(hs & 0xffff_f000);
77	129	const c = (z - df * df) / (s + df);
78		const w = r32(z) * s + c;
79		return 2 * (df + w);
	130	const w = rationalApproxBinary32(z) * s + c;
	131	return 2.0 * (df + w);
80	132	}
81	133
82		fn r64(z: f64) f64 {
	134	fn rationalApproxBinary64(z: f64) f64 {
83	135	const pS0: f64 = 1.66666666666666657415e-01;
84	136	const pS1: f64 = -3.25565818622400915405e-01;
85	137	const pS2: f64 = 2.01212532134862925881e-01;
...	...	@@ -96,91 +148,292 @@ fn r64(z: f64) f64 {
96	148	return p / q;
97	149	}
98	150
99		fn acos64(x: f64) f64 {
	151	fn acosBinary64(x: f64) f64 {
100	152	const pio2_hi: f64 = 1.57079632679489655800e+00;
101	153	const pio2_lo: f64 = 6.12323399573676603587e-17;
102	154
103		const ux = @as(u64, @bitCast(x));
104		const hx = @as(u32, @intCast(ux >> 32));
105		const ix = hx & 0x7FFFFFFF;
	155	const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32);
	156	const ix: u32 = hx & 0x7fff_ffff;
106	157
107	158	// \|x\| >= 1 or nan
108		if (ix >= 0x3FF00000) {
109		const lx = @as(u32, @intCast(ux & 0xFFFFFFFF));
110
111		// acos(1) = 0, acos(-1) = pi
112		if ((ix - 0x3FF00000) \| lx == 0) {
	159	if (ix >= 0x3ff0_0000) {
	160	const lx: u32 = @truncate(@as(u64, @bitCast(x)));
	161	if ((ix - 0x3ff0_0000 \| lx) == 0) {
113	162	if (hx >> 31 != 0) {
114		return 2 * pio2_hi + 0x1.0p-120;
115		} else {
116		return 0;
	163	return 2.0 * pio2_hi + 0x1.0p-120;
117	164	}
	165	return 0.0;
118	166	}
119
120		return (x - x) / 0;
	167	return 0.0 / (x - x);
121	168	}
122	169
123	170	// \|x\| < 0.5
124		if (ix < 0x3FE00000) {
	171	if (ix < 0x3fe0_0000) {
125	172	// \|x\| < 2^(-57)
126		if (ix <= 0x3C600000) {
	173	if (ix <= 0x3c60_0000) {
127	174	return pio2_hi + 0x1.0p-120;
128		} else {
129		return pio2_hi - (x - (pio2_lo - x * r64(x * x)));
130	175	}
	176	return pio2_hi - (x - (pio2_lo - x * rationalApproxBinary64(x * x)));
131	177	}
132	178
133	179	// x < -0.5
134	180	if (hx >> 31 != 0) {
135	181	const z = (1.0 + x) * 0.5;
136	182	const s = @sqrt(z);
137		const w = r64(z) * s - pio2_lo;
	183	const w = rationalApproxBinary64(z) * s - pio2_lo;
138	184	return 2 * (pio2_hi - (s + w));
139	185	}
140	186
141	187	// x > 0.5
142	188	const z = (1.0 - x) * 0.5;
143	189	const s = @sqrt(z);
144		const jx = @as(u64, @bitCast(s));
145		const df = @as(f64, @bitCast(jx & 0xFFFFFFFF00000000));
	190	const df: f64 = @bitCast(@as(u64, @bitCast(s)) & 0xffff_ffff_0000_0000);
146	191	const c = (z - df * df) / (s + df);
147		const w = r64(z) * s + c;
148		return 2 * (df + w);
	192	const w = rationalApproxBinary64(z) * s + c;
	193	return 2.0 * (df + w);
149	194	}
150	195
151		test acos {
152		try expect(acos(@as(f32, 0.0)) == acos32(0.0));
153		try expect(acos(@as(f64, 0.0)) == acos64(0.0));
	196	fn rationalApproxExtended80(z: f80) f80 {
	197	const pS0: f80 = 1.66666666666666666631e-01;
	198	const pS1: f80 = -4.16313987993683104320e-01;
	199	const pS2: f80 = 3.69068046323246813704e-01;
	200	const pS3: f80 = -1.36213932016738603108e-01;
	201	const pS4: f80 = 1.78324189708471965733e-02;
	202	const pS5: f80 = -2.19216428382605211588e-04;
	203	const pS6: f80 = -7.10526623669075243183e-06;
	204	const qS1: f80 = -2.94788392796209867269e+00;
	205	const qS2: f80 = 3.27309890266528636716e+00;
	206	const qS3: f80 = -1.68285799854822427013e+00;
	207	const qS4: f80 = 3.90699412641738801874e-01;
	208	const qS5: f80 = -3.14365703596053263322e-02;
	209
	210	const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6))))));
	211	const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5))));
	212	return p / q;
154	213	}
155	214
156		test acos32 {
157		const epsilon = 0.000001;
	215	fn acosExtended80(x: f80) f80 {
	216	const pio2_hi: f80 = 1.57079632679489661926;
	217	const pio2_lo: f80 = -2.50827880633416601173e-20;
	218
	219	const hx: u80 = @bitCast(x);
	220	const se: u16 = @truncate(hx >> 64);
	221	const e = se & 0x7fff;
158	222
159		try expect(math.approxEqAbs(f32, acos32(0.0), 1.570796, epsilon));
160		try expect(math.approxEqAbs(f32, acos32(0.2), 1.369438, epsilon));
161		try expect(math.approxEqAbs(f32, acos32(0.3434), 1.220262, epsilon));
162		try expect(math.approxEqAbs(f32, acos32(0.5), 1.047198, epsilon));
163		try expect(math.approxEqAbs(f32, acos32(0.8923), 0.468382, epsilon));
164		try expect(math.approxEqAbs(f32, acos32(-0.2), 1.772154, epsilon));
	223	// \|x\| >= 1 or nan
	224	if (e >= 0x3fff) {
	225	if (x == 1.0) {
	226	return 0.0;
	227	}
	228	if (x == -1.0) {
	229	return 2.0 * pio2_hi + 0x1p-120;
	230	}
	231	return 0.0 / (x - x);
	232	}
	233	// \|x\| < 0.5
	234	if (e < 0x3fff - 1) {
	235	if (e < 0x3fff - math.floatFractionalBits(f80)) {
	236	return pio2_hi + 0x1p-120;
	237	}
	238	return pio2_hi - (rationalApproxExtended80(x * x) * x - pio2_lo + x);
	239	}
	240	// x < -0.5
	241	if (se >> 15 != 0) {
	242	const z = (1 + x) * 0.5;
	243	const s = @sqrt(z);
	244	return 2.0 * (pio2_hi - (rationalApproxExtended80(z) * s - pio2_lo + s));
	245	}
	246	// x > 0.5
	247	const z = (1.0 - x) * 0.5;
	248	const s = @sqrt(z);
	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	return 2.0 * (rationalApproxExtended80(z) * s + c + f);
	253	}
	254
	255	fn rationalApproxBinary128(z: f128) f128 {
	256	const pS0: f128 = 1.66666666666666666666666666666700314e-01;
	257	const pS1: f128 = -7.32816946414566252574527475428622708e-01;
	258	const pS2: f128 = 1.34215708714992334609030036562143589e+00;
	259	const pS3: f128 = -1.32483151677116409805070261790752040e+00;
	260	const pS4: f128 = 7.61206183613632558824485341162121989e-01;
	261	const pS5: f128 = -2.56165783329023486777386833928147375e-01;
	262	const pS6: f128 = 4.80718586374448793411019434585413855e-02;
	263	const pS7: f128 = -4.42523267167024279410230886239774718e-03;
	264	const pS8: f128 = 1.44551535183911458253205638280410064e-04;
	265	const pS9: f128 = -2.10558957916600254061591040482706179e-07;
	266	const qS1: f128 = -4.84690167848739751544716485245697428e+00;
	267	const qS2: f128 = 9.96619113536172610135016921140206980e+00;
	268	const qS3: f128 = -1.13177895428973036660836798461641458e+01;
	269	const qS4: f128 = 7.74004374389488266169304117714658761e+00;
	270	const qS5: f128 = -3.25871986053534084709023539900339905e+00;
	271	const qS6: f128 = 8.27830318881232209752469022352928864e-01;
	272	const qS7: f128 = -1.18768052702942805423330715206348004e-01;
	273	const qS8: f128 = 8.32600764660522313269101537926539470e-03;
	274	const qS9: f128 = -1.99407384882605586705979504567947007e-04;
	275
	276	const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9)))))))));
	277	const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9))))))));
	278	return p / q;
	279	}
	280
	281	fn acosBinary128(x: f128) f128 {
	282	const pio2_hi: f128 = 1.57079632679489661923132169163975140;
	283	const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35;
	284
	285	const hx: u128 = @bitCast(x);
	286	const se: u16 = @truncate(hx >> 112);
	287	const e = se & 0x7fff;
	288
	289	// \|x\| >= 1 or nan
	290	if (e >= 0x3fff) {
	291	if (x == 1.0) {
	292	return 0.0;
	293	}
	294	if (x == -1.0) {
	295	return 2 * pio2_hi + 0x1p-120;
	296	}
	297	return 0.0 / (x - x);
	298	}
	299	// \|x\| < 0.5
	300	if (e < 0x3fff - 1) {
	301	if (e < 0x3fff - math.floatFractionalBits(f128)) {
	302	return pio2_hi + 0x1p-120;
	303	}
	304	return pio2_hi - (rationalApproxBinary128(x * x) * x - pio2_lo + x);
	305	}
	306	// x < -0.5
	307	if (se >> 15 != 0) {
	308	const z = (1 + x) * 0.5;
	309	const s = @sqrt(z);
	310	return 2 * (pio2_hi - (rationalApproxBinary128(z) * s - pio2_lo + s));
	311	}
	312	// x > 0.5
	313	const z = (1.0 - x) * 0.5;
	314	const s = @sqrt(z);
	315	const hs: u128 = @bitCast(s);
	316	const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000);
	317	const c = (z - f * f) / (s + f);
	318	return 2.0 * (rationalApproxBinary128(z) * s + c + f);
	319	}
	320
	321	test "acosBinary16.special" {
	322	try testing.expectApproxEqAbs(acosBinary16(0x0p+0), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0));
	323	try testing.expectApproxEqAbs(acosBinary16(-0x1p+0), 0x1.92p1, math.floatEpsAt(f16, 0x1.92p1));
	324	try testing.expectEqual(acosBinary16(0x1p+0), 0x0p+0);
	325	try testing.expect(math.isNan(acosBinary16(0x1.004p0)));
	326	try testing.expect(math.isNan(acosBinary16(-0x1.004p0)));
	327	try testing.expect(math.isNan(acosBinary16(math.inf(f16))));
	328	try testing.expect(math.isNan(acosBinary16(-math.inf(f16))));
	329	try testing.expect(math.isNan(acosBinary16(math.nan(f16))));
	330	}
	331
	332	test "acosBinary16" {
	333	try testing.expectApproxEqAbs(acosBinary16(0x1.db4p-5), 0x1.834p0, math.floatEpsAt(f16, 0x1.834p0));
	334	try testing.expectApproxEqAbs(acosBinary16(-0x1.068p-2), 0x1.d48p0, math.floatEpsAt(f16, 0x1.d48p0));
	335	try testing.expectApproxEqAbs(acosBinary16(-0x1.2c4p-3), 0x1.b7cp0, math.floatEpsAt(f16, 0x1.b7cp0));
	336	try testing.expectApproxEqAbs(acosBinary16(0x1.65p-3), 0x1.654p0, math.floatEpsAt(f16, 0x1.654p0));
	337	try testing.expectApproxEqAbs(acosBinary16(0x1.dfcp-1), 0x1.6d8p-2, math.floatEpsAt(f16, 0x1.6d8p-2));
	338	try testing.expectApproxEqAbs(acosBinary16(-0x1.764p-1), 0x1.32p1, math.floatEpsAt(f16, 0x1.32p1));
	339	try testing.expectApproxEqAbs(acosBinary16(0x1.b18p-3), 0x1.5b8p0, math.floatEpsAt(f16, 0x1.5b8p0));
	340	try testing.expectApproxEqAbs(acosBinary16(0x1.5acp-3), 0x1.668p0, math.floatEpsAt(f16, 0x1.668p0));
	341	try testing.expectApproxEqAbs(acosBinary16(-0x1.18cp-1), 0x1.134p1, math.floatEpsAt(f16, 0x1.134p1));
	342	try testing.expectApproxEqAbs(acosBinary16(-0x1.03p-1), 0x1.0dp1, math.floatEpsAt(f16, 0x1.0dp1));
165	343	}
166	344
167		test acos64 {
168		const epsilon = 0.000001;
	345	test "acosBinary32.special" {
	346	try testing.expectApproxEqAbs(acosBinary32(0x0p+0), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0));
	347	try testing.expectApproxEqAbs(acosBinary32(-0x1p+0), 0x1.921fb6p+1, math.floatEpsAt(f32, 0x1.921fb6p+1));
	348	try testing.expectEqual(acosBinary32(0x1p+0), 0x0p+0);
	349	try testing.expect(math.isNan(acosBinary32(0x1.000002p+0)));
	350	try testing.expect(math.isNan(acosBinary32(-0x1.000002p+0)));
	351	try testing.expect(math.isNan(acosBinary32(math.inf(f32))));
	352	try testing.expect(math.isNan(acosBinary32(-math.inf(f32))));
	353	try testing.expect(math.isNan(acosBinary32(math.nan(f32))));
	354	}
	355
	356	test "acosBinary32" {
	357	try testing.expectApproxEqAbs(acosBinary32(-0x1.13284cp-2), 0x1.d7c4e6p+0, math.floatEpsAt(f32, 0x1.d7c4e6p+0));
	358	try testing.expectApproxEqAbs(acosBinary32(0x1.6ca8ep-1), 0x1.8e6756p-1, math.floatEpsAt(f32, 0x1.8e6756p-1));
	359	try testing.expectApproxEqAbs(acosBinary32(0x1.c2ca6p-1), 0x1.f9d74cp-2, math.floatEpsAt(f32, 0x1.f9d74cp-2));
	360	try testing.expectApproxEqAbs(acosBinary32(-0x1.55f12p-1), 0x1.26abdcp+1, math.floatEpsAt(f32, 0x1.26abdcp+1));
	361	try testing.expectApproxEqAbs(acosBinary32(-0x1.15679ep-2), 0x1.d85a44p+0, math.floatEpsAt(f32, 0x1.d85a44p+0));
	362	try testing.expectApproxEqAbs(acosBinary32(-0x1.41e132p-5), 0x1.9c2f68p+0, math.floatEpsAt(f32, 0x1.9c2f68p+0));
	363	try testing.expectApproxEqAbs(acosBinary32(0x1.281b0ep-1), 0x1.e881bp-1, math.floatEpsAt(f32, 0x1.e881bp-1));
	364	try testing.expectApproxEqAbs(acosBinary32(0x1.b5ce34p-1), 0x1.1713f6p-1, math.floatEpsAt(f32, 0x1.1713f6p-1));
	365	try testing.expectApproxEqAbs(acosBinary32(-0x1.583482p-3), 0x1.bd5accp+0, math.floatEpsAt(f32, 0x1.bd5accp+0));
	366	try testing.expectApproxEqAbs(acosBinary32(-0x1.ea8224p-1), 0x1.6ce7d8p+1, math.floatEpsAt(f32, 0x1.6ce7d8p+1));
	367	}
	368
	369	test "acosBinary64.special" {
	370	try testing.expectApproxEqAbs(acosBinary64(0x0p+0), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
	371	try testing.expectApproxEqAbs(acosBinary64(-0x1p+0), 0x1.921fb54442d18p+1, math.floatEpsAt(f64, 0x1.921fb54442d18p+1));
	372	try testing.expectEqual(acosBinary64(0x1p+0), 0x0p+0);
	373	try testing.expect(math.isNan(acosBinary64(0x1.0000000000001p+0)));
	374	try testing.expect(math.isNan(acosBinary64(-0x1.0000000000001p+0)));
	375	try testing.expect(math.isNan(acosBinary64(math.inf(f64))));
	376	try testing.expect(math.isNan(acosBinary64(-math.inf(f64))));
	377	try testing.expect(math.isNan(acosBinary64(math.nan(f64))));
	378	}
	379
	380	test "acosBinary64" {
	381	try testing.expectApproxEqAbs(acosBinary64(-0x1.13284b2b5006dp-2), 0x1.d7c4e61020905p+0, math.floatEpsAt(f64, 0x1.d7c4e61020905p+0));
	382	try testing.expectApproxEqAbs(acosBinary64(0x1.6ca8dfb825911p-1), 0x1.8e6756e27c366p-1, math.floatEpsAt(f64, 0x1.8e6756e27c366p-1));
	383	try testing.expectApproxEqAbs(acosBinary64(0x1.c2ca609de7505p-1), 0x1.f9d748eaf956p-2, math.floatEpsAt(f64, 0x1.f9d748eaf956p-2));
	384	try testing.expectApproxEqAbs(acosBinary64(-0x1.55f11fba96889p-1), 0x1.26abdc68d07aap+1, math.floatEpsAt(f64, 0x1.26abdc68d07aap+1));
	385	try testing.expectApproxEqAbs(acosBinary64(-0x1.15679e27084ddp-2), 0x1.d85a44ea44fe4p+0, math.floatEpsAt(f64, 0x1.d85a44ea44fe4p+0));
	386	try testing.expectApproxEqAbs(acosBinary64(-0x1.41e131b093c41p-5), 0x1.9c2f688eee8abp+0, math.floatEpsAt(f64, 0x1.9c2f688eee8abp+0));
	387	try testing.expectApproxEqAbs(acosBinary64(0x1.281b0d18455f5p-1), 0x1.e881b1d4eb2a1p-1, math.floatEpsAt(f64, 0x1.e881b1d4eb2a1p-1));
	388	try testing.expectApproxEqAbs(acosBinary64(0x1.b5ce34a51b239p-1), 0x1.1713f567a87efp-1, math.floatEpsAt(f64, 0x1.1713f567a87efp-1));
	389	try testing.expectApproxEqAbs(acosBinary64(-0x1.583481079de4dp-3), 0x1.bd5acbe8fcc59p+0, math.floatEpsAt(f64, 0x1.bd5acbe8fcc59p+0));
	390	try testing.expectApproxEqAbs(acosBinary64(-0x1.ea8223103b871p-1), 0x1.6ce7d66f628e5p+1, math.floatEpsAt(f64, 0x1.6ce7d66f628e5p+1));
	391	}
	392
	393	test "acosExtended80.special" {
	394	try testing.expectApproxEqAbs(acosExtended80(0x0p+0), 0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0));
	395	try testing.expectApproxEqAbs(acosExtended80(-0x1p+0), 0x1.921fb54442d1846ap+1, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+1));
	396	try testing.expectEqual(acosExtended80(0x1p+0), 0x0p+0);
	397	try testing.expect(math.isNan(acosExtended80(0x1.0000000000000002p+0)));
	398	try testing.expect(math.isNan(acosExtended80(-0x1.0000000000000002p+0)));
	399	try testing.expect(math.isNan(acosExtended80(math.inf(f80))));
	400	try testing.expect(math.isNan(acosExtended80(-math.inf(f80))));
	401	try testing.expect(math.isNan(acosExtended80(math.nan(f80))));
	402	}
169	403
170		try expect(math.approxEqAbs(f64, acos64(0.0), 1.570796, epsilon));
171		try expect(math.approxEqAbs(f64, acos64(0.2), 1.369438, epsilon));
172		try expect(math.approxEqAbs(f64, acos64(0.3434), 1.220262, epsilon));
173		try expect(math.approxEqAbs(f64, acos64(0.5), 1.047198, epsilon));
174		try expect(math.approxEqAbs(f64, acos64(0.8923), 0.468382, epsilon));
175		try expect(math.approxEqAbs(f64, acos64(-0.2), 1.772154, epsilon));
	404	test "acosExtended80" {
	405	try testing.expectApproxEqAbs(acosExtended80(0x1.72068a321edc8804p-1), 0x1.86b349040d28f794p-1, math.floatEpsAt(f80, 0x1.86b349040d28f794p-1));
	406	try testing.expectApproxEqAbs(acosExtended80(-0x1.06d0a467d22977ecp-2), 0x1.d4923ade73ec379cp0, math.floatEpsAt(f80, 0x1.d4923ade73ec379cp0));
	407	try testing.expectApproxEqAbs(acosExtended80(0x1.77d21385faa9798ap-3), 0x1.62e0e8898c6d04f2p0, math.floatEpsAt(f80, 0x1.62e0e8898c6d04f2p0));
	408	try testing.expectApproxEqAbs(acosExtended80(-0x1.73ee3e8bc2a44dbep-1), 0x1.3123cbcd5dc4bd58p1, math.floatEpsAt(f80, 0x1.3123cbcd5dc4bd58p1));
	409	try testing.expectApproxEqAbs(acosExtended80(0x1.0a2dd1f6ffcf668ap-1), 0x1.062a6d562df2d316p0, math.floatEpsAt(f80, 0x1.062a6d562df2d316p0));
	410	try testing.expectApproxEqAbs(acosExtended80(0x1.8e835c490a3aff9ep-3), 0x1.5ffd68b520aa55fap0, math.floatEpsAt(f80, 0x1.5ffd68b520aa55fap0));
	411	try testing.expectApproxEqAbs(acosExtended80(0x1.add20cdc1565064cp-3), 0x1.5bfe6cabda700684p0, math.floatEpsAt(f80, 0x1.5bfe6cabda700684p0));
	412	try testing.expectApproxEqAbs(acosExtended80(0x1.21986d43727fca72p-8), 0x1.90fe1c993b571924p0, math.floatEpsAt(f80, 0x1.90fe1c993b571924p0));
	413	try testing.expectApproxEqAbs(acosExtended80(0x1.d61e0b3fae6a0564p-2), 0x1.18044ccc626e7f9ep0, math.floatEpsAt(f80, 0x1.18044ccc626e7f9ep0));
	414	try testing.expectApproxEqAbs(acosExtended80(-0x1.171e7c4a41883ccap-4), 0x1.a39513b6c16532b4p0, math.floatEpsAt(f80, 0x1.a39513b6c16532b4p0));
176	415	}
177	416
178		test "acos32.special" {
179		try expect(math.isNan(acos32(-2)));
180		try expect(math.isNan(acos32(1.5)));
	417	test "acosBinary128.special" {
	418	try testing.expectApproxEqAbs(acosBinary128(0x0p+0), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0));
	419	try testing.expectApproxEqAbs(acosBinary128(-0x1p+0), 0x1.921fb54442d18469898cc51701b8p1, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p1));
	420	try testing.expectEqual(acosBinary128(0x1p+0), 0x0p+0);
	421	try testing.expect(math.isNan(acosBinary128(0x1.0000000000000000000000000001p0)));
	422	try testing.expect(math.isNan(acosBinary128(-0x1.0000000000000000000000000001p0)));
	423	try testing.expect(math.isNan(acosBinary128(math.inf(f128))));
	424	try testing.expect(math.isNan(acosBinary128(-math.inf(f128))));
	425	try testing.expect(math.isNan(acosBinary128(math.nan(f128))));
181	426	}
182	427
183		test "acos64.special" {
184		try expect(math.isNan(acos64(-2)));
185		try expect(math.isNan(acos64(1.5)));
	428	test "acosBinary128" {
	429	try testing.expectApproxEqAbs(acosBinary128(-0x1.511bdb99a3c4373bedf834ef4f68p-1), 0x1.250e9a58f049eeafa99db4360c88p1, math.floatEpsAt(f128, 0x1.250e9a58f049eeafa99db4360c88p1));
	430	try testing.expectApproxEqAbs(acosBinary128(-0x1.5879cc3ad6dfd2a52e9891c69808p-1), 0x1.2786664b1c676c99437b68590004p1, math.floatEpsAt(f128, 0x1.2786664b1c676c99437b68590004p1));
	431	try testing.expectApproxEqAbs(acosBinary128(0x1.3f988ba64a7eb97a751c5f0b3077p-1), 0x1.cb190cd361c7c03a09c470b4caebp-1, math.floatEpsAt(f128, 0x1.cb190cd361c7c03a09c470b4caebp-1));
	432	try testing.expectApproxEqAbs(acosBinary128(-0x1.3f2d96c7768e4c4fa02315727959p-1), 0x1.1f373be697880111758f582b1a96p1, math.floatEpsAt(f128, 0x1.1f373be697880111758f582b1a96p1));
	433	try testing.expectApproxEqAbs(acosBinary128(0x1.fad303c2e28c1f4d8f9fd0e5686fp-2), 0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0, math.floatEpsAt(f128, 0x1.0d92fd2a0a6ca3e4853c1de9ea6ap0));
	434	try testing.expectApproxEqAbs(acosBinary128(0x1.ddde322bd1a2ee50c5ba30c9c617p-2), 0x1.15d4b306e16fbf9ea4f29e82b154p0, math.floatEpsAt(f128, 0x1.15d4b306e16fbf9ea4f29e82b154p0));
	435	try testing.expectApproxEqAbs(acosBinary128(-0x1.b02f6adefcbeb1d48666b827ff17p-1), 0x1.49b0a0355a5539052388e8a6dc11p1, math.floatEpsAt(f128, 0x1.49b0a0355a5539052388e8a6dc11p1));
	436	try testing.expectApproxEqAbs(acosBinary128(0x1.c8581cce7cd3f6efab0fc60d9b7dp-2), 0x1.1be0b757f4cef022f5d2422b9c78p0, math.floatEpsAt(f128, 0x1.1be0b757f4cef022f5d2422b9c78p0));
	437	try testing.expectApproxEqAbs(acosBinary128(-0x1.bf887b8c4e33cbef59993056f3dep-1), 0x1.513270e671db2d840f20b0186c2cp1, math.floatEpsAt(f128, 0x1.513270e671db2d840f20b0186c2cp1));
	438	try testing.expectApproxEqAbs(acosBinary128(0x1.0c0f600ab6f9c84c6102942044cep-3), 0x1.70851a509f0e8bfbe780aa8f29f9p0, math.floatEpsAt(f128, 0x1.70851a509f0e8bfbe780aa8f29f9p0));
186	439	}

lib/std/math/asin.zig+344-90

...	...	@@ -3,10 +3,14 @@
3	3	//
4	4	// https://git.musl-libc.org/cgit/musl/tree/src/math/asinf.c
5	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
6	7
7	8	const std = @import("../std.zig");
8	9	const math = std.math;
9		const expect = std.testing.expect;
	10	const mem = std.mem;
	11	const testing = std.testing;
	12	const builtin = @import("builtin");
	13	const native_endian = builtin.cpu.arch.endian();
10	14
11	15	/// Returns the arc-sin of x.
12	16	///
...	...	@@ -16,62 +20,101 @@ const expect = std.testing.expect;
16	20	pub fn asin(x: anytype) @TypeOf(x) {
17	21	const T = @TypeOf(x);
18	22	return switch (T) {
19		f32 => asin32(x),
20		f64 => asin64(x),
	23	f16 => asinBinary16(x),
	24	f32 => asinBinary32(x),
	25	f64 => asinBinary64(x),
	26	f80 => asinExtended80(x),
	27	f128 => asinBinary128(x),
21	28	else => @compileError("asin not implemented for " ++ @typeName(T)),
22	29	};
23	30	}
24	31
25		fn r32(z: f32) f32 {
26		const pS0 = 1.6666586697e-01;
27		const pS1 = -4.2743422091e-02;
28		const pS2 = -8.6563630030e-03;
29		const qS1 = -7.0662963390e-01;
	32	fn approxBinary16(z: f32) f32 {
	33	const S0: f32 = 1.0000001e0;
	34	const S1: f32 = 1.6664918e-1;
	35	const S2: f32 = 7.55022e-2;
	36	const S3: f32 = 3.9513987e-2;
	37	const S4: f32 = 5.0883885e-2;
	38	return S0 + z * (S1 + z * (S2 + z * (S3 + z * S4)));
	39	}
	40
	41	fn asinBinary16(x: f16) f16 {
	42	const pio2: f32 = math.pi / 2.0;
	43
	44	const hx: u16 = @bitCast(x);
	45	const ix = hx & 0x7fff;
	46
	47	// \|x\| >= 1
	48	if (ix >= 0x3c00) {
	49	// \|x\| == 1
	50	if (ix == 0x3c00) {
	51	// asin(+-1) = +-pi/2 with inexact
	52	return @floatCast(x * pio2 + 0x1.0p-120);
	53	}
	54	// asin(\|x\| > 1) is nan
	55	return 0.0 / (x - x);
	56	}
	57
	58	// \|x\| < 0.5
	59	if (ix < 0x3800) {
	60	return @floatCast(x * approxBinary16(x * x));
	61	}
	62
	63	// 1 > \|x\| >= 0.5
	64	const z = (1.0 - @abs(x)) * 0.5;
	65	const s = @sqrt(z);
	66	const x_local = pio2 - 2.0 * s * approxBinary16(z);
	67	if (hx >> 15 != 0) {
	68	return @floatCast(-x_local);
	69	}
	70	return @floatCast(x_local);
	71	}
	72
	73	fn rationalApproxBinary32(z: f32) f32 {
	74	const pS0: f32 = 1.6666586697e-01;
	75	const pS1: f32 = -4.2743422091e-02;
	76	const pS2: f32 = -8.6563630030e-03;
	77	const qS1: f32 = -7.0662963390e-01;
30	78
31	79	const p = z * (pS0 + z * (pS1 + z * pS2));
32	80	const q = 1.0 + z * qS1;
33	81	return p / q;
34	82	}
35	83
36		fn asin32(x: f32) f32 {
37		const pio2 = 1.570796326794896558e+00;
	84	fn asinBinary32(x: f32) f32 {
	85	const pio2: f64 = 1.570796326794896558e+00;
38	86
39		const hx: u32 = @as(u32, @bitCast(x));
40		const ix: u32 = hx & 0x7FFFFFFF;
	87	const hx: u32 = @bitCast(x);
	88	const ix = hx & 0x7fff_ffff;
41	89
42	90	// \|x\| >= 1
43		if (ix >= 0x3F800000) {
44		// \|x\| >= 1
45		if (ix == 0x3F800000) {
46		return x * pio2 + 0x1.0p-120; // asin(+-1) = +-pi/2 with inexact
47		} else {
48		return math.nan(f32); // asin(\|x\| > 1) is nan
	91	if (ix >= 0x3f80_0000) {
	92	// \|x\| == 1
	93	if (ix == 0x3f80_0000) {
	94	// asin(+-1) = +-pi/2 with inexact
	95	return @floatCast(@as(f64, @floatCast(x)) * pio2 + 0x1.0p-120);
49	96	}
	97	// asin(\|x\| > 1) is nan
	98	return 0.0 / (x - x);
50	99	}
51	100
52	101	// \|x\| < 0.5
53		if (ix < 0x3F000000) {
	102	if (ix < 0x3f00_0000) {
54	103	// 0x1p-126 <= \|x\| < 0x1p-12
55		if (ix < 0x39800000 and ix >= 0x00800000) {
	104	if (ix < 0x3980_0000 and ix >= 0x0080_0000) {
56	105	return x;
57		} else {
58		return x + x * r32(x * x);
59	106	}
	107	return x + x * rationalApproxBinary32(x * x);
60	108	}
61	109
62	110	// 1 > \|x\| >= 0.5
63		const z = (1 - @abs(x)) * 0.5;
64		const s = @sqrt(z);
65		const fx = pio2 - 2 * (s + s * r32(z));
66
67		if (hx >> 31 != 0) {
68		return -fx;
69		} else {
70		return fx;
71		}
	111	const z = (1.0 - @abs(x)) * 0.5;
	112	const s: f64 = @floatCast(@sqrt(z));
	113	const x_local: f32 = @floatCast(pio2 - 2.0 * (s + s * @as(f64, @floatCast(rationalApproxBinary32(z)))));
	114	return if (hx >> 31 != 0) -x_local else x_local;
72	115	}
73	116
74		fn r64(z: f64) f64 {
	117	fn rationalApproxBinary64(z: f64) f64 {
75	118	const pS0: f64 = 1.66666666666666657415e-01;
76	119	const pS1: f64 = -3.25565818622400915405e-01;
77	120	const pS2: f64 = 2.01212532134862925881e-01;
...	...	@@ -88,96 +131,307 @@ fn r64(z: f64) f64 {
88	131	return p / q;
89	132	}
90	133
91		fn asin64(x: f64) f64 {
	134	fn asinBinary64(x: f64) f64 {
92	135	const pio2_hi: f64 = 1.57079632679489655800e+00;
93	136	const pio2_lo: f64 = 6.12323399573676603587e-17;
94	137
95		const ux = @as(u64, @bitCast(x));
96		const hx = @as(u32, @intCast(ux >> 32));
97		const ix = hx & 0x7FFFFFFF;
	138	const hx: u32 = @intCast(@as(u64, @bitCast(x)) >> 32);
	139	const ix = hx & 0x7fffffff;
98	140
99	141	// \|x\| >= 1 or nan
100		if (ix >= 0x3FF00000) {
101		const lx = @as(u32, @intCast(ux & 0xFFFFFFFF));
102
	142	if (ix >= 0x3ff0_0000) {
	143	const lx: u32 = @truncate(@as(u64, @bitCast(x)));
103	144	// asin(1) = +-pi/2 with inexact
104		if ((ix - 0x3FF00000) \| lx == 0) {
	145	if ((ix - 0x3ff0_0000 \| lx) == 0) {
105	146	return x * pio2_hi + 0x1.0p-120;
106		} else {
107		return math.nan(f64);
108	147	}
	148	return 0.0 / (x - x);
109	149	}
110	150
111	151	// \|x\| < 0.5
112		if (ix < 0x3FE00000) {
	152	if (ix < 0x3fe0_0000) {
113	153	// if 0x1p-1022 <= \|x\| < 0x1p-26 avoid raising overflow
114		if (ix < 0x3E500000 and ix >= 0x00100000) {
	154	if (ix < 0x3e50_0000 and ix >= 0x0010_0000) {
115	155	return x;
116		} else {
117		return x + x * r64(x * x);
118	156	}
	157	return x + x * rationalApproxBinary64(x * x);
119	158	}
120	159
121	160	// 1 > \|x\| >= 0.5
122		const z = (1 - @abs(x)) * 0.5;
	161	const z = (1.0 - @abs(x)) * 0.5;
123	162	const s = @sqrt(z);
124		const r = r64(z);
125		var fx: f64 = undefined;
126
	163	const r = rationalApproxBinary64(z);
127	164	// \|x\| > 0.975
128		if (ix >= 0x3FEF3333) {
129		fx = pio2_hi - 2 * (s + s * r);
130		} else {
131		const jx = @as(u64, @bitCast(s));
132		const df = @as(f64, @bitCast(jx & 0xFFFFFFFF00000000));
133		const c = (z - df * df) / (s + df);
134		fx = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * df));
	165	if (ix >= 0x3fef_3333) {
	166	const x_local = pio2_hi - (2 * (s + s * r) - pio2_lo);
	167	return if (hx >> 31 != 0) -x_local else x_local;
	168	}
	169	// f+c = sqrt(z)
	170	const hs: u64 = @bitCast(s);
	171	const f: f64 = @bitCast(hs & 0xffff_ffff_0000_0000);
	172	const c: f64 = (z - f * f) / (s + f);
	173	const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
	174	return if (hx >> 31 != 0) -x_local else x_local;
	175	}
	176
	177	fn rationalApproxExtended80(z: f80) f80 {
	178	const pS0: f80 = 1.66666666666666666631e-01;
	179	const pS1: f80 = -4.16313987993683104320e-01;
	180	const pS2: f80 = 3.69068046323246813704e-01;
	181	const pS3: f80 = -1.36213932016738603108e-01;
	182	const pS4: f80 = 1.78324189708471965733e-02;
	183	const pS5: f80 = -2.19216428382605211588e-04;
	184	const pS6: f80 = -7.10526623669075243183e-06;
	185	const qS1: f80 = -2.94788392796209867269e+00;
	186	const qS2: f80 = 3.27309890266528636716e+00;
	187	const qS3: f80 = -1.68285799854822427013e+00;
	188	const qS4: f80 = 3.90699412641738801874e-01;
	189	const qS5: f80 = -3.14365703596053263322e-02;
	190
	191	const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * pS6))))));
	192	const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * qS5))));
	193	return p / q;
	194	}
	195
	196	fn asinExtended80(x: f80) f80 {
	197	const pio2_hi: f80 = 1.57079632679489661926;
	198	const pio2_lo: f80 = -2.50827880633416601173e-20;
	199
	200	const hx: u80 = @bitCast(x);
	201	const se: u16 = @truncate(hx >> 64);
	202	const e = se & 0x7fff;
	203	const sign = se >> 15 != 0;
	204
	205	// \|x\| >= 1 or nan
	206	if (e >= 0x3fff) {
	207	// asin(+-1)=+-pi/2 with inexact
	208	if (x == 1.0 or x == -1.0) {
	209	return x * pio2_hi + 0x1p-120;
	210	}
	211	return 0.0 / (x - x);
135	212	}
136	213
137		if (hx >> 31 != 0) {
138		return -fx;
139		} else {
140		return fx;
	214	// \|x\| < 0.5
	215	if (e < 0x3fff - 1) {
	216	if (e < 0x3fff - (math.floatMantissaBits(f80) + 1) / 2) {
	217	// return x with inexact if x!=0
	218	mem.doNotOptimizeAway(x + 0x1p120);
	219	return x;
	220	}
	221	return x + x * rationalApproxExtended80(x * x);
141	222	}
	223
	224	// 1 > \|x\| >= 0.5
	225	const z = (1.0 - @abs(x)) * 0.5;
	226	const s = @sqrt(z);
	227	const r = rationalApproxExtended80(z);
	228
	229	const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff);
	230	if ((m >> 56) >= 0xf7) {
	231	const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
	232	return if (sign) -x_local else x_local;
	233	}
	234
	235	const hs: u80 = @bitCast(s);
	236	const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_0000_0000);
	237	const c = (z - f * f) / (s + f);
	238	const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
	239	return if (sign) -x_local else x_local;
142	240	}
143	241
144		test asin {
145		try expect(asin(@as(f32, 0.0)) == asin32(0.0));
146		try expect(asin(@as(f64, 0.0)) == asin64(0.0));
	242	fn rationalApproxBinary128(z: f128) f128 {
	243	const pS0: f128 = 1.66666666666666666666666666666700314e-01;
	244	const pS1: f128 = -7.32816946414566252574527475428622708e-01;
	245	const pS2: f128 = 1.34215708714992334609030036562143589e+00;
	246	const pS3: f128 = -1.32483151677116409805070261790752040e+00;
	247	const pS4: f128 = 7.61206183613632558824485341162121989e-01;
	248	const pS5: f128 = -2.56165783329023486777386833928147375e-01;
	249	const pS6: f128 = 4.80718586374448793411019434585413855e-02;
	250	const pS7: f128 = -4.42523267167024279410230886239774718e-03;
	251	const pS8: f128 = 1.44551535183911458253205638280410064e-04;
	252	const pS9: f128 = -2.10558957916600254061591040482706179e-07;
	253	const qS1: f128 = -4.84690167848739751544716485245697428e+00;
	254	const qS2: f128 = 9.96619113536172610135016921140206980e+00;
	255	const qS3: f128 = -1.13177895428973036660836798461641458e+01;
	256	const qS4: f128 = 7.74004374389488266169304117714658761e+00;
	257	const qS5: f128 = -3.25871986053534084709023539900339905e+00;
	258	const qS6: f128 = 8.27830318881232209752469022352928864e-01;
	259	const qS7: f128 = -1.18768052702942805423330715206348004e-01;
	260	const qS8: f128 = 8.32600764660522313269101537926539470e-03;
	261	const qS9: f128 = -1.99407384882605586705979504567947007e-04;
	262
	263	const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * (pS5 + z * (pS6 + z * (pS7 + z * (pS8 + z * pS9)))))))));
	264	const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * (qS4 + z * (qS5 + z * (qS6 + z * (qS7 + z * (qS8 + z * qS9))))))));
	265	return p / q;
147	266	}
148	267
149		test asin32 {
150		const epsilon = 0.000001;
	268	fn asinBinary128(x: f128) f128 {
	269	const pio2_hi: f128 = 1.57079632679489661923132169163975140;
	270	const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35;
	271
	272	const hx: u128 = @bitCast(x);
	273	const se: u16 = @truncate(hx >> 112);
	274	const e = se & 0x7fff;
	275	const sign = se >> 15 != 0;
	276
	277	// \|x\| >= 1 or nan
	278	if (e >= 0x3fff) {
	279	// asin(+-1)=+-pi/2 with inexact
	280	if (x == 1.0 or x == -1.0) {
	281	return x * pio2_hi + 0x1p-120;
	282	}
	283	return 0.0 / (x - x);
	284	}
151	285
152		try expect(math.approxEqAbs(f32, asin32(0.0), 0.0, epsilon));
153		try expect(math.approxEqAbs(f32, asin32(0.2), 0.201358, epsilon));
154		try expect(math.approxEqAbs(f32, asin32(-0.2), -0.201358, epsilon));
155		try expect(math.approxEqAbs(f32, asin32(0.3434), 0.350535, epsilon));
156		try expect(math.approxEqAbs(f32, asin32(0.5), 0.523599, epsilon));
157		try expect(math.approxEqAbs(f32, asin32(0.8923), 1.102415, epsilon));
	286	// \|x\| < 0.5
	287	if (e < 0x3fff - 1) {
	288	if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) {
	289	// return x with inexact if x!=0
	290	mem.doNotOptimizeAway(x + 0x1p120);
	291	return x;
	292	}
	293	return x + x * rationalApproxBinary128(x * x);
	294	}
	295
	296	// 1 > \|x\| >= 0.5
	297	const z = (1.0 - @abs(x)) * 0.5;
	298	const s = @sqrt(z);
	299	const r = rationalApproxBinary128(z);
	300
	301	const top: u16 = @truncate((hx >> 96) & 0x0000_ffff);
	302	if (top >= 0xee00) {
	303	const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
	304	return if (sign) -x_local else x_local;
	305	}
	306
	307	const hs: u128 = @bitCast(s);
	308	const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000);
	309	const c = (z - f * f) / (s + f);
	310	const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
	311	return if (sign) -x_local else x_local;
158	312	}
159	313
160		test asin64 {
161		const epsilon = 0.000001;
	314	test "asinBinary16.special" {
	315	try testing.expectApproxEqAbs(asinBinary16(0x1p+0), 0x1.92p0, math.floatEpsAt(f16, 0x1.92p0));
	316	try testing.expectApproxEqAbs(asinBinary16(-0x1p+0), -0x1.92p0, math.floatEpsAt(f16, -0x1.92p0));
	317	try testing.expectEqual(asinBinary16(0x0p+0), 0x0p+0);
	318	try testing.expectEqual(asinBinary16(-0x0p+0), 0x0p+0);
	319	try testing.expect(math.isNan(asinBinary16(0x1.004p0)));
	320	try testing.expect(math.isNan(asinBinary16(-0x1.004p0)));
	321	try testing.expect(math.isNan(asinBinary16(math.inf(f16))));
	322	try testing.expect(math.isNan(asinBinary16(-math.inf(f16))));
	323	try testing.expect(math.isNan(asinBinary16(math.nan(f16))));
	324	}
	325
	326	test "asinBinary16" {
	327	try testing.expectApproxEqAbs(asinBinary16(-0x1.e4cp-6), -0x1.e4cp-6, math.floatEpsAt(f16, -0x1.e4cp-6));
	328	try testing.expectApproxEqAbs(asinBinary16(0x1.d68p-1), 0x1.2a8p0, math.floatEpsAt(f16, 0x1.2a8p0));
	329	try testing.expectApproxEqAbs(asinBinary16(-0x1.a4cp-1), -0x1.eep-1, math.floatEpsAt(f16, -0x1.eep-1));
	330	try testing.expectApproxEqAbs(asinBinary16(-0x1.0a4p-2), -0x1.0d4p-2, math.floatEpsAt(f16, -0x1.0d4p-2));
	331	try testing.expectApproxEqAbs(asinBinary16(0x1.28cp-1), 0x1.3c8p-1, math.floatEpsAt(f16, 0x1.3c8p-1));
	332	try testing.expectApproxEqAbs(asinBinary16(0x1.284p-3), 0x1.298p-3, math.floatEpsAt(f16, 0x1.298p-3));
	333	try testing.expectApproxEqAbs(asinBinary16(-0x1.574p-1), -0x1.784p-1, math.floatEpsAt(f16, -0x1.784p-1));
	334	try testing.expectApproxEqAbs(asinBinary16(-0x1.4ccp-1), -0x1.6a4p-1, math.floatEpsAt(f16, -0x1.6a4p-1));
	335	try testing.expectApproxEqAbs(asinBinary16(0x1.a18p-1), 0x1.e84p-1, math.floatEpsAt(f16, 0x1.e84p-1));
	336	try testing.expectApproxEqAbs(asinBinary16(0x1.7a8p-2), 0x1.83cp-2, math.floatEpsAt(f16, 0x1.83cp-2));
	337	}
	338
	339	test "asinBinary32.special" {
	340	try testing.expectApproxEqAbs(asinBinary32(0x1p+0), 0x1.921fb6p+0, math.floatEpsAt(f32, 0x1.921fb6p+0));
	341	try testing.expectApproxEqAbs(asinBinary32(-0x1p+0), -0x1.921fb6p+0, math.floatEpsAt(f32, -0x1.921fb6p+0));
	342	try testing.expectEqual(asinBinary32(0x0p+0), 0x0p+0);
	343	try testing.expectEqual(asinBinary32(-0x0p+0), 0x0p+0);
	344	try testing.expect(math.isNan(asinBinary32(0x1.000002p+0)));
	345	try testing.expect(math.isNan(asinBinary32(-0x1.000002p+0)));
	346	try testing.expect(math.isNan(asinBinary32(math.inf(f32))));
	347	try testing.expect(math.isNan(asinBinary32(-math.inf(f32))));
	348	try testing.expect(math.isNan(asinBinary32(math.nan(f32))));
	349	}
	350
	351	test "asinBinary32" {
	352	try testing.expectApproxEqAbs(asinBinary32(-0x1.4c2906p-4), -0x1.4c868p-4, math.floatEpsAt(f32, -0x1.4c868p-4));
	353	try testing.expectApproxEqAbs(asinBinary32(0x1.05fcfap-1), 0x1.130648p-1, math.floatEpsAt(f32, 0x1.130648p-1));
	354	try testing.expectApproxEqAbs(asinBinary32(0x1.fab976p-2), 0x1.090abcp-1, math.floatEpsAt(f32, 0x1.090abcp-1));
	355	try testing.expectApproxEqAbs(asinBinary32(0x1.8b4b8cp-1), 0x1.c39fa2p-1, math.floatEpsAt(f32, 0x1.c39fa2p-1));
	356	try testing.expectApproxEqAbs(asinBinary32(0x1.7117c2p-1), 0x1.9c332p-1, math.floatEpsAt(f32, 0x1.9c332p-1));
	357	try testing.expectApproxEqAbs(asinBinary32(0x1.e5e112p-5), 0x1.e62a1cp-5, math.floatEpsAt(f32, 0x1.e62a1cp-5));
	358	try testing.expectApproxEqAbs(asinBinary32(-0x1.07673p-2), -0x1.0a65dep-2, math.floatEpsAt(f32, -0x1.0a65dep-2));
	359	try testing.expectApproxEqAbs(asinBinary32(-0x1.2108dep-2), -0x1.25046p-2, math.floatEpsAt(f32, -0x1.25046p-2));
	360	try testing.expectApproxEqAbs(asinBinary32(-0x1.4e6e6cp-1), -0x1.6c6f0cp-1, math.floatEpsAt(f32, -0x1.6c6f0cp-1));
	361	try testing.expectApproxEqAbs(asinBinary32(0x1.22a16ap-1), 0x1.350f7ap-1, math.floatEpsAt(f32, 0x1.350f7ap-1));
	362	}
	363
	364	test "asinBinary64.special" {
	365	try testing.expectApproxEqAbs(asinBinary64(0x1p+0), 0x1.921fb54442d18p+0, math.floatEpsAt(f64, 0x1.921fb54442d18p+0));
	366	try testing.expectApproxEqAbs(asinBinary64(-0x1p+0), -0x1.921fb54442d18p+0, math.floatEpsAt(f64, -0x1.921fb54442d18p+0));
	367	try testing.expectEqual(asinBinary64(0x0p+0), 0x0p+0);
	368	try testing.expectEqual(asinBinary64(-0x0p+0), 0x0p+0);
	369	try testing.expect(math.isNan(asinBinary64(0x1.000002p+0)));
	370	try testing.expect(math.isNan(asinBinary64(-0x1.000002p+0)));
	371	try testing.expect(math.isNan(asinBinary64(math.inf(f64))));
	372	try testing.expect(math.isNan(asinBinary64(-math.inf(f64))));
	373	try testing.expect(math.isNan(asinBinary64(math.nan(f64))));
	374	}
	375
	376	test "asinBinary64" {
	377	try testing.expectApproxEqAbs(asinBinary64(0x1.e674fba3e40d5p-2), 0x1.fae86c5941692p-2, math.floatEpsAt(f64, 0x1.fae86c5941692p-2));
	378	try testing.expectApproxEqAbs(asinBinary64(-0x1.30fd0566fd979p-1), -0x1.46b6ad730c93ap-1, math.floatEpsAt(f64, -0x1.46b6ad730c93ap-1));
	379	try testing.expectApproxEqAbs(asinBinary64(0x1.6444a25abfeaap-2), 0x1.6be0be8074eep-2, math.floatEpsAt(f64, 0x1.6be0be8074eep-2));
	380	try testing.expectApproxEqAbs(asinBinary64(0x1.40a53228d1a13p-1), 0x1.5a7e98f53f717p-1, math.floatEpsAt(f64, 0x1.5a7e98f53f717p-1));
	381	try testing.expectApproxEqAbs(asinBinary64(-0x1.ccc6d64845cfdp-1), -0x1.1ea2602d14e8p0, math.floatEpsAt(f64, -0x1.1ea2602d14e8p0));
	382	try testing.expectApproxEqAbs(asinBinary64(-0x1.94bd91b7fc74bp-1), -0x1.d2c2634193158p-1, math.floatEpsAt(f64, -0x1.d2c2634193158p-1));
	383	try testing.expectApproxEqAbs(asinBinary64(-0x1.8d741b5797fccp-2), -0x1.982d5f1895d2p-2, math.floatEpsAt(f64, -0x1.982d5f1895d2p-2));
	384	try testing.expectApproxEqAbs(asinBinary64(-0x1.3e8e7e15881c5p-3), -0x1.3fdaf7dfdc864p-3, math.floatEpsAt(f64, -0x1.3fdaf7dfdc864p-3));
	385	try testing.expectApproxEqAbs(asinBinary64(-0x1.88222d8ab8ca9p-2), -0x1.9269540735b7bp-2, math.floatEpsAt(f64, -0x1.9269540735b7bp-2));
	386	try testing.expectApproxEqAbs(asinBinary64(-0x1.41c0e9babcbd2p-2), -0x1.474c4c6625527p-2, math.floatEpsAt(f64, -0x1.474c4c6625527p-2));
	387	}
	388
	389	test "asinExtended80.special" {
	390	try testing.expectApproxEqAbs(asinExtended80(0x1p+0), 0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, 0x1.921fb54442d1846ap+0));
	391	try testing.expectApproxEqAbs(asinExtended80(-0x1p+0), -0x1.921fb54442d1846ap+0, math.floatEpsAt(f80, -0x1.921fb54442d1846ap+0));
	392	try testing.expectEqual(asinExtended80(0x0p+0), 0x0p+0);
	393	try testing.expectEqual(asinExtended80(-0x0p+0), 0x0p+0);
	394	try testing.expect(math.isNan(asinExtended80(0x1.0000000000000002p+0)));
	395	try testing.expect(math.isNan(asinExtended80(-0x1.0000000000000002p+0)));
	396	try testing.expect(math.isNan(asinExtended80(math.inf(f80))));
	397	try testing.expect(math.isNan(asinExtended80(-math.inf(f80))));
	398	try testing.expect(math.isNan(asinExtended80(math.nan(f80))));
	399	}
162	400
163		try expect(math.approxEqAbs(f64, asin64(0.0), 0.0, epsilon));
164		try expect(math.approxEqAbs(f64, asin64(0.2), 0.201358, epsilon));
165		try expect(math.approxEqAbs(f64, asin64(-0.2), -0.201358, epsilon));
166		try expect(math.approxEqAbs(f64, asin64(0.3434), 0.350535, epsilon));
167		try expect(math.approxEqAbs(f64, asin64(0.5), 0.523599, epsilon));
168		try expect(math.approxEqAbs(f64, asin64(0.8923), 1.102415, epsilon));
	401	test "asinExtended80" {
	402	try testing.expectApproxEqAbs(asinExtended80(0x1.63cf98bc52ce0da8p-9), 0x1.63cfb560149daa9p-9, math.floatEpsAt(f80, 0x1.63cfb560149daa9p-9));
	403	try testing.expectApproxEqAbs(asinExtended80(-0x1.0473756f7ae930dp-1), -0x1.113cbacd8cd1b96cp-1, math.floatEpsAt(f80, -0x1.113cbacd8cd1b96cp-1));
	404	try testing.expectApproxEqAbs(asinExtended80(-0x1.2310057e005cc288p-2), -0x1.2721b231d197b064p-2, math.floatEpsAt(f80, -0x1.2721b231d197b064p-2));
	405	try testing.expectApproxEqAbs(asinExtended80(0x1.f13b03bd685d96eap-1), 0x1.547c408c5d2b05aap0, math.floatEpsAt(f80, 0x1.547c408c5d2b05aap0));
	406	try testing.expectApproxEqAbs(asinExtended80(-0x1.d5c507e3ef84041cp-1), -0x1.296b76bfadbb5cecp0, math.floatEpsAt(f80, -0x1.296b76bfadbb5cecp0));
	407	try testing.expectApproxEqAbs(asinExtended80(0x1.8222cbc9147153d8p-1), 0x1.b572da8729a84f2ap-1, math.floatEpsAt(f80, 0x1.b572da8729a84f2ap-1));
	408	try testing.expectApproxEqAbs(asinExtended80(-0x1.42c9e6b4a088a246p-11), -0x1.42c9e80ac0524dap-11, math.floatEpsAt(f80, -0x1.42c9e80ac0524dap-11));
	409	try testing.expectApproxEqAbs(asinExtended80(-0x1.8f78d49deadb521cp-3), -0x1.920ca86aef6c3028p-3, math.floatEpsAt(f80, -0x1.920ca86aef6c3028p-3));
	410	try testing.expectApproxEqAbs(asinExtended80(-0x1.ab98792783515774p-2), -0x1.b91cb4f7204d92fp-2, math.floatEpsAt(f80, -0x1.b91cb4f7204d92fp-2));
	411	try testing.expectApproxEqAbs(asinExtended80(-0x1.104fe30cef6800aap-1), -0x1.1f20815fdc4c5304p-1, math.floatEpsAt(f80, -0x1.1f20815fdc4c5304p-1));
169	412	}
170	413
171		test "asin32.special" {
172		try expect(math.isPositiveZero(asin32(0.0)));
173		try expect(math.isNegativeZero(asin32(-0.0)));
174		try expect(math.isNan(asin32(-2)));
175		try expect(math.isNan(asin32(1.5)));
	414	test "asinBinary128.special" {
	415	try testing.expectApproxEqAbs(asinBinary128(0x1p+0), 0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, 0x1.921fb54442d18469898cc51701b8p0));
	416	try testing.expectApproxEqAbs(asinBinary128(-0x1p+0), -0x1.921fb54442d18469898cc51701b8p0, math.floatEpsAt(f128, -0x1.921fb54442d18469898cc51701b8p0));
	417	try testing.expectEqual(asinBinary128(0x0p+0), 0x0p+0);
	418	try testing.expectEqual(asinBinary128(-0x0p+0), 0x0p+0);
	419	try testing.expect(math.isNan(asinBinary128(0x1.0000000000000000000000000001p0)));
	420	try testing.expect(math.isNan(asinBinary128(-0x1.0000000000000000000000000001p0)));
	421	try testing.expect(math.isNan(asinBinary128(math.inf(f128))));
	422	try testing.expect(math.isNan(asinBinary128(-math.inf(f128))));
	423	try testing.expect(math.isNan(asinBinary128(math.nan(f128))));
176	424	}
177	425
178		test "asin64.special" {
179		try expect(math.isPositiveZero(asin64(0.0)));
180		try expect(math.isNegativeZero(asin64(-0.0)));
181		try expect(math.isNan(asin64(-2)));
182		try expect(math.isNan(asin64(1.5)));
	426	test "asinBinary128" {
	427	try testing.expectApproxEqAbs(asinBinary128(0x1.85868ce287ca0196b01c25fec5ffp-3), 0x1.87e9c740d7837f8e8fa667988fbep-3, math.floatEpsAt(f128, 0x1.87e9c740d7837f8e8fa667988fbep-3));
	428	try testing.expectApproxEqAbs(asinBinary128(0x1.8718d6d30b4daed08d04ef59f478p-1), 0x1.bd11a474e864213b48e0f005f1f4p-1, math.floatEpsAt(f128, 0x1.bd11a474e864213b48e0f005f1f4p-1));
	429	try testing.expectApproxEqAbs(asinBinary128(0x1.11a67640cd7f0ba5d5e362f3abfap-1), 0x1.20b56f8b42649fe72d1f8d68a378p-1, math.floatEpsAt(f128, 0x1.20b56f8b42649fe72d1f8d68a378p-1));
	430	try testing.expectApproxEqAbs(asinBinary128(-0x1.bd13bf14a9dce22188e52650daa7p-1), -0x1.0dc3a7ddb9736e5ad699bf338566p0, math.floatEpsAt(f128, -0x1.0dc3a7ddb9736e5ad699bf338566p0));
	431	try testing.expectApproxEqAbs(asinBinary128(-0x1.dee0bc217fc462af57c484eefa71p-2), -0x1.f250716038f70fa50a5826c03802p-2, math.floatEpsAt(f128, -0x1.f250716038f70fa50a5826c03802p-2));
	432	try testing.expectApproxEqAbs(asinBinary128(-0x1.ea7df9139371c10b9d6fd2bbccd3p-1), -0x1.47a8b4cdd327f90056722feddbabp0, math.floatEpsAt(f128, -0x1.47a8b4cdd327f90056722feddbabp0));
	433	try testing.expectApproxEqAbs(asinBinary128(0x1.04aaea6de3b5a616460702f26dfcp-2), 0x1.079178d52be662dec67e2cd7f6e9p-2, math.floatEpsAt(f128, 0x1.079178d52be662dec67e2cd7f6e9p-2));
	434	try testing.expectApproxEqAbs(asinBinary128(-0x1.c7ea85e6b61be666435a7d99444cp-1), -0x1.192df5a8d71702cf1e27014887b2p0, math.floatEpsAt(f128, -0x1.192df5a8d71702cf1e27014887b2p0));
	435	try testing.expectApproxEqAbs(asinBinary128(-0x1.6e210214e40edf6c8479998189d1p-1), -0x1.97f1092fd94ac0fdfddae2e1222bp-1, math.floatEpsAt(f128, -0x1.97f1092fd94ac0fdfddae2e1222bp-1));
	436	try testing.expectApproxEqAbs(asinBinary128(-0x1.95061bf93ed6986a45d20f0e1064p-3), -0x1.97b62bc5ae6512093828828325e1p-3, math.floatEpsAt(f128, -0x1.97b62bc5ae6512093828828325e1p-3));
183	437	}