Commit d34b868bcf759a81275fe51b9c2faeb15bf7bedd | zig

author	Andrew Kelley <andrew@ziglang.org> 2026-04-03 12:04:32+02:00
committer	Andrew Kelley <andrew@ziglang.org> 2026-04-03 12:04:32+02:00
log	d34b868bcf759a81275fe51b9c2faeb15bf7bedd
tree	989d270f4df185757eea5cd3c31cf6479c55cc0c
parent	ce3f254526609a9b2088d81070903107c969bb81
parent	4ccac1de416bff8846352ef3bf2fb592c5a419a9

Merge pull request '`libzigc`: Implement 12 more math functions' (#31604) from mihael/zig:libzigc/more-math-functions-2 into master

Reviewed-on: https://codeberg.org/ziglang/zig/pulls/31604
Reviewed-by: Andrew Kelley <andrew@ziglang.org>

54 files changed, 1257 insertions(+), 1700 deletions(-)

lib
- c
  - math.zig
- compiler_rt
- libc
  - mingw
    math
    arm-common
    sincosl.c
    arm
    sincos.S
    arm64
    rint.c
    rintf.c
    sincos.S
    frexpf.c
    frexpl.c
    x86
    cos.def.h
    cosl.c
    cosl_internal.S
    sin.def.h
    sinl.c
    sinl_internal.S
    tanl.S
  - musl
    src
    math
    __cosl.c
    __sinl.c
    __tanl.c
    aarch64
    lrint.c
    lrintf.c
    rintf.c
    cosl.c
    finite.c
    finitef.c
    frexp.c
    frexpf.c
    frexpl.c
    i386
    lrint.c
    lrintf.c
    rintf.c
    lrint.c
    lrintf.c
    powerpc64
    lrint.c
    lrintf.c
    rintf.c
    s390x
    rintf.c
    sincosl.c
    sinl.c
    tanl.c
    x32
    lrint.s
    lrintf.s
    x86_64
    lrint.c
    lrintf.c
src
- libs
test
- libc.zig

lib/c/math.zig+111-22

...	...	@@ -36,6 +36,8 @@ comptime {
36	36
37	37	if (builtin.target.isMinGW() or builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
38	38	symbol(&coshf, "coshf");
	39	symbol(&frexpf, "frexpf");
	40	symbol(&frexpl, "frexpl");
39	41	symbol(&hypotf, "hypotf");
40	42	symbol(&hypotl, "hypotl");
41	43	symbol(&modff, "modff");
...	...	@@ -46,6 +48,11 @@ comptime {
46	48	symbol(&tanhf, "tanhf");
47	49	}
48	50
	51	if (builtin.target.isMinGW() or builtin.target.isMuslLibC()) {
	52	symbol(&rint, "rint");
	53	symbol(&rintf, "rintf");
	54	}
	55
49	56	if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
50	57	symbol(&acos, "acos");
51	58	symbol(&acosf, "acosf");
...	...	@@ -60,7 +67,12 @@ comptime {
60	67	symbol(&exp10, "exp10");
61	68	symbol(&exp10f, "exp10f");
62	69	symbol(&fdim, "fdim");
	70	symbol(&finite, "finite");
	71	symbol(&finitef, "finitef");
	72	symbol(&frexp, "frexp");
63	73	symbol(&hypot, "hypot");
	74	symbol(&lrint, "lrint");
	75	symbol(&lrintf, "lrintf");
64	76	symbol(&modf, "modf");
65	77	symbol(&pow, "pow");
66	78	symbol(&pow10, "pow10");
...	...	@@ -71,7 +83,6 @@ comptime {
71	83	if (builtin.target.isMuslLibC()) {
72	84	symbol(&copysign, "copysign");
73	85	symbol(&copysignf, "copysignf");
74		symbol(&rint, "rint");
75	86	}
76	87
77	88	symbol(&copysignl, "copysignl");
...	...	@@ -161,6 +172,45 @@ fn fdim(x: f64, y: f64) callconv(.c) f64 {
161	172	return 0;
162	173	}
163	174
	175	fn finite(x: f64) callconv(.c) c_int {
	176	return if (math.isFinite(x)) 1 else 0;
	177	}
	178
	179	fn finitef(x: f32) callconv(.c) c_int {
	180	return if (math.isFinite(x)) 1 else 0;
	181	}
	182
	183	fn frexpGeneric(comptime T: type, x: T, e: *c_int) T {
	184	// libc expects `*e` to be unspecified in this case; an unspecified C value
	185	// should be a valid value of the relevant type, yet Zig's std
	186	// implementation sets it to `undefined` -- which can even be nonsense
	187	// according to the type (int). Therefore, we're setting it to a valid
	188	// int value in Zig -- a zero.
	189	//
	190	// This mirrors the handling of infinities, where libc also expects
	191	// unspecified for the value of `*e` and Zig std sets it to a zero.
	192	if (math.isNan(x)) {
	193	e.* = 0;
	194	return x;
	195	}
	196
	197	const r = math.frexp(x);
	198	e.* = r.exponent;
	199	return r.significand;
	200	}
	201
	202	fn frexp(x: f64, e: *c_int) callconv(.c) f64 {
	203	return frexpGeneric(f64, x, e);
	204	}
	205
	206	fn frexpf(x: f32, e: *c_int) callconv(.c) f32 {
	207	return frexpGeneric(f32, x, e);
	208	}
	209
	210	fn frexpl(x: c_longdouble, e: *c_int) callconv(.c) c_longdouble {
	211	return frexpGeneric(c_longdouble, x, e);
	212	}
	213
164	214	fn hypot(x: f64, y: f64) callconv(.c) f64 {
165	215	return math.hypot(x, y);
166	216	}
...	...	@@ -185,6 +235,14 @@ fn isnanl(x: c_longdouble) callconv(.c) c_int {
185	235	return if (math.isNan(x)) 1 else 0;
186	236	}
187	237
	238	fn lrint(x: f64) callconv(.c) c_long {
	239	return @intFromFloat(rint(x));
	240	}
	241
	242	fn lrintf(x: f32) callconv(.c) c_long {
	243	return @intFromFloat(rintf(x));
	244	}
	245
188	246	fn modfGeneric(comptime T: type, x: T, iptr: *T) T {
189	247	if (math.isNegativeInf(x)) {
190	248	iptr.* = -math.inf(T);
...	...	@@ -299,7 +357,7 @@ fn pow10f(x: f32) callconv(.c) f32 {
299	357	}
300	358
301	359	fn rint(x: f64) callconv(.c) f64 {
302		const toint: f64 = 1.0 / @as(f64, math.floatEps(f64));
	360	const toint: f64 = 1.0 / math.floatEps(f64);
303	361	const a: u64 = @bitCast(x);
304	362	const e = a >> 52 & 0x7ff;
305	363	const s = a >> 63;
...	...	@@ -319,39 +377,70 @@ fn rint(x: f64) callconv(.c) f64 {
319	377	return y;
320	378	}
321	379
322		test "rint" {
	380	fn rintf(x: f32) callconv(.c) f32 {
	381	const toint: f32 = 1.0 / math.floatEps(f32);
	382	const a: u32 = @bitCast(x);
	383	const e = a >> 23 & 0xff;
	384	const s = a >> 31;
	385	var y: f32 = undefined;
	386
	387	if (e >= 0x7f + 23) {
	388	return x;
	389	}
	390
	391	if (s == 1) {
	392	y = x - toint + toint;
	393	} else {
	394	y = x + toint - toint;
	395	}
	396
	397	if (y == 0) {
	398	return if (s == 1) -0.0 else 0;
	399	}
	400	return y;
	401	}
	402
	403	fn testRint(comptime T: type) !void {
	404	const f = switch (T) {
	405	f32 => rintf,
	406	f64 => rint,
	407	else => @compileError("rint not implemented for" ++ @typeName(T)),
	408	};
	409
323	410	// Positive numbers round correctly
324		try expectEqual(@as(f64, 42.0), rint(42.2));
325		try expectEqual(@as(f64, 42.0), rint(41.8));
	411	try expectEqual(@as(T, 42.0), f(42.2));
	412	try expectEqual(@as(T, 42.0), f(41.8));
326	413
327	414	// Negative numbers round correctly
328		try expectEqual(@as(f64, -6.0), rint(-5.9));
329		try expectEqual(@as(f64, -6.0), rint(-6.1));
	415	try expectEqual(@as(T, -6.0), f(-5.9));
	416	try expectEqual(@as(T, -6.0), f(-6.1));
330	417
331	418	// No rounding needed test
332		try expectEqual(@as(f64, 5.0), rint(5.0));
333		try expectEqual(@as(f64, -10.0), rint(-10.0));
334		try expectEqual(@as(f64, 0.0), rint(0.0));
	419	try expectEqual(@as(T, 5.0), f(5.0));
	420	try expectEqual(@as(T, -10.0), f(-10.0));
	421	try expectEqual(@as(T, 0.0), f(0.0));
335	422
336	423	// Very large numbers return unchanged
337		const large: f64 = 9007199254740992.0; // 2^53
338		try expectEqual(large, rint(large));
339		try expectEqual(-large, rint(-large));
	424	const large: T = 9007199254740992.0; // 2^53
	425	try expectEqual(large, f(large));
	426	try expectEqual(-large, f(-large));
340	427
341	428	// Small positive numbers round to zero
342		const pos_result = rint(0.3);
343		try expectEqual(@as(f64, 0.0), pos_result);
344		try expect(@as(u64, @bitCast(pos_result)) == 0);
	429	const pos_result = f(0.3);
	430	try expect(math.isPositiveZero(pos_result));
345	431
346	432	// Small negative numbers round to negative zero
347		const neg_result = rint(-0.3);
348		try expectEqual(@as(f64, 0.0), neg_result);
349		const bits: u64 = @bitCast(neg_result);
350		try expect((bits >> 63) == 1);
	433	const neg_result = f(-0.3);
	434	try expect(math.isNegativeZero(neg_result));
351	435
352	436	// Exact half rounds to nearest even (banker's rounding)
353		try expectEqual(@as(f64, 2.0), rint(2.5));
354		try expectEqual(@as(f64, 4.0), rint(3.5));
	437	try expectEqual(@as(T, 2.0), f(2.5));
	438	try expectEqual(@as(T, 4.0), f(3.5));
	439	}
	440
	441	test "rint" {
	442	try testRint(f32);
	443	try testRint(f64);
355	444	}
356	445
357	446	fn tanh(x: f64) callconv(.c) f64 {

lib/compiler_rt/cos.zig+141-52

...	...	@@ -1,19 +1,30 @@
	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/cosf.c
	5	//! https://git.musl-libc.org/cgit/musl/tree/src/math/cos.c
	6	//! https://git.musl-libc.org/cgit/musl/tree/src/math/cosl.c
	7
1	8	const std = @import("std");
2	9	const math = std.math;
3	10	const mem = std.mem;
4	11	const expect = std.testing.expect;
	12	const expectApproxEqAbs = std.testing.expectApproxEqAbs;
5	13
6	14	const compiler_rt = @import("../compiler_rt.zig");
7	15	const symbol = @import("../compiler_rt.zig").symbol;
8	16	const trig = @import("trig.zig");
9	17	const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
10	18	const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
	19	const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
	20	const ld = @import("long_double.zig");
11	21
12	22	comptime {
13		symbol(&__cosh, "__cosh");
	23	symbol(&cosh, "__cosh");
	24	symbol(&cosl, "__cosl");
14	25	symbol(&cosf, "cosf");
15	26	symbol(&cos, "cos");
16		symbol(&__cosx, "__cosx");
	27	symbol(&cosx, "__cosx");
17	28	if (compiler_rt.want_ppc_abi) {
18	29	symbol(&cosq, "cosf128");
19	30	}
...	...	@@ -21,7 +32,7 @@ comptime {
21	32	symbol(&cosl, "cosl");
22	33	}
23	34
24		pub fn __cosh(a: f16) callconv(.c) f16 {
	35	pub fn cosh(a: f16) callconv(.c) f16 {
25	36	// TODO: more efficient implementation
26	37	return @floatCast(cosf(a));
27	38	}
...	...	@@ -43,27 +54,27 @@ pub fn cosf(x: f32) callconv(.c) f32 {
43	54	if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
44	55	return 1.0;
45	56	}
46		return trig.__cosdf(x);
	57	return trig.cosdf(x);
47	58	}
48	59	if (ix <= 0x407b53d1) { // \|x\| ~<= 5*pi/4
49	60	if (ix > 0x4016cbe3) { // \|x\| ~> 3*pi/4
50		return -trig.__cosdf(if (sign) x + c2pio2 else x - c2pio2);
	61	return -trig.cosdf(if (sign) x + c2pio2 else x - c2pio2);
51	62	} else {
52	63	if (sign) {
53		return trig.__sindf(x + c1pio2);
	64	return trig.sindf(x + c1pio2);
54	65	} else {
55		return trig.__sindf(c1pio2 - x);
	66	return trig.sindf(c1pio2 - x);
56	67	}
57	68	}
58	69	}
59	70	if (ix <= 0x40e231d5) { // \|x\| ~<= 9*pi/4
60	71	if (ix > 0x40afeddf) { // \|x\| ~> 7*pi/4
61		return trig.__cosdf(if (sign) x + c4pio2 else x - c4pio2);
	72	return trig.cosdf(if (sign) x + c4pio2 else x - c4pio2);
62	73	} else {
63	74	if (sign) {
64		return trig.__sindf(-x - c3pio2);
	75	return trig.sindf(-x - c3pio2);
65	76	} else {
66		return trig.__sindf(x - c3pio2);
	77	return trig.sindf(x - c3pio2);
67	78	}
68	79	}
69	80	}
...	...	@@ -76,10 +87,10 @@ pub fn cosf(x: f32) callconv(.c) f32 {
76	87	var y: f64 = undefined;
77	88	const n = rem_pio2f(x, &y);
78	89	return switch (n & 3) {
79		0 => trig.__cosdf(y),
80		1 => trig.__sindf(-y),
81		2 => -trig.__cosdf(y),
82		else => trig.__sindf(y),
	90	0 => trig.cosdf(y),
	91	1 => trig.sindf(-y),
	92	2 => -trig.cosdf(y),
	93	else => trig.sindf(y),
83	94	};
84	95	}
85	96
...	...	@@ -94,7 +105,7 @@ pub fn cos(x: f64) callconv(.c) f64 {
94	105	if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x + 0x1p120);
95	106	return 1.0;
96	107	}
97		return trig.__cos(x, 0);
	108	return trig.cos(x, 0);
98	109	}
99	110
100	111	// cos(Inf or NaN) is NaN
...	...	@@ -105,66 +116,144 @@ pub fn cos(x: f64) callconv(.c) f64 {
105	116	var y: [2]f64 = undefined;
106	117	const n = rem_pio2(x, &y);
107	118	return switch (n & 3) {
108		0 => trig.__cos(y[0], y[1]),
109		1 => -trig.__sin(y[0], y[1], 1),
110		2 => -trig.__cos(y[0], y[1]),
111		else => trig.__sin(y[0], y[1], 1),
	119	0 => trig.cos(y[0], y[1]),
	120	1 => -trig.sin(y[0], y[1], 1),
	121	2 => -trig.cos(y[0], y[1]),
	122	else => trig.sin(y[0], y[1], 1),
112	123	};
113	124	}
114	125
115		pub fn __cosx(a: f80) callconv(.c) f80 {
116		// TODO: more efficient implementation
117		return @floatCast(cosq(a));
	126	pub fn cosx(x: f80) callconv(.c) f80 {
	127	const se = ld.signExponent(x) & 0x7fff;
	128	if (se == 0x7fff) {
	129	return x - x;
	130	}
	131
	132	if (@abs(x) < trig.pi_4) {
	133	if (se < 0x3fff - math.floatMantissaBits(f80)) {
	134	// raise inexact if x!=0
	135	return 1.0 + x;
	136	}
	137	return trig.cosx(x, 0.0);
	138	}
	139
	140	var y: [2]f80 = undefined;
	141	const n = rem_pio2l(f80, x, &y);
	142	return switch (n & 3) {
	143	0 => trig.cosx(y[0], y[1]),
	144	1 => -trig.sinx(y[0], y[1], 1),
	145	2 => -trig.cosx(y[0], y[1]),
	146	else => trig.sinx(y[0], y[1], 1),
	147	};
118	148	}
119	149
120		pub fn cosq(a: f128) callconv(.c) f128 {
121		// TODO: more correct implementation
122		return cos(@floatCast(a));
	150	pub fn cosq(x: f128) callconv(.c) f128 {
	151	const se = ld.signExponent(x) & 0x7fff;
	152	if (se == 0x7fff) {
	153	return x - x;
	154	}
	155
	156	if (@abs(x) < trig.pi_4) {
	157	if (se < 0x3fff - math.floatMantissaBits(f128)) {
	158	// raise inexact if x!=0
	159	return 1.0 + x;
	160	}
	161	return trig.cosq(x, 0.0);
	162	}
	163
	164	var y: [2]f128 = undefined;
	165	const n = rem_pio2l(f128, x, &y);
	166	return switch (n & 3) {
	167	0 => trig.cosq(y[0], y[1]),
	168	1 => -trig.sinq(y[0], y[1], 1),
	169	2 => -trig.cosq(y[0], y[1]),
	170	else => trig.sinq(y[0], y[1], 1),
	171	};
123	172	}
124	173
125	174	pub fn cosl(x: c_longdouble) callconv(.c) c_longdouble {
126	175	switch (@typeInfo(c_longdouble).float.bits) {
127		16 => return __cosh(x),
	176	16 => return cosh(x),
128	177	32 => return cosf(x),
129	178	64 => return cos(x),
130		80 => return __cosx(x),
	179	80 => return cosx(x),
131	180	128 => return cosq(x),
132	181	else => @compileError("unreachable"),
133	182	}
134	183	}
135	184
136		test "cos32" {
137		const epsilon = 0.00001;
	185	fn testCosSpecial(comptime T: type) !void {
	186	const f = switch (T) {
	187	f32 => cosf,
	188	f64 => cos,
	189	f80 => cosx,
	190	f128 => cosq,
	191	else => @compileError("unimplemented"),
	192	};
138	193
139		try expect(math.approxEqAbs(f32, cosf(0.0), 1.0, epsilon));
140		try expect(math.approxEqAbs(f32, cosf(0.2), 0.980067, epsilon));
141		try expect(math.approxEqAbs(f32, cosf(0.8923), 0.627623, epsilon));
142		try expect(math.approxEqAbs(f32, cosf(1.5), 0.070737, epsilon));
143		try expect(math.approxEqAbs(f32, cosf(-1.5), 0.070737, epsilon));
144		try expect(math.approxEqAbs(f32, cosf(37.45), 0.969132, epsilon));
145		try expect(math.approxEqAbs(f32, cosf(89.123), 0.400798, epsilon));
	194	try expect(f(0.0) == 1.0);
	195	try expect(f(-0.0) == 1.0);
	196	try expect(math.isNan(f(math.inf(T))));
	197	try expect(math.isNan(f(-math.inf(T))));
	198	try expect(math.isNan(f(math.nan(T))));
146	199	}
147	200
148		test "cos64" {
149		const epsilon = 0.000001;
150
151		try expect(math.approxEqAbs(f64, cos(0.0), 1.0, epsilon));
152		try expect(math.approxEqAbs(f64, cos(0.2), 0.980067, epsilon));
153		try expect(math.approxEqAbs(f64, cos(0.8923), 0.627623, epsilon));
154		try expect(math.approxEqAbs(f64, cos(1.5), 0.070737, epsilon));
155		try expect(math.approxEqAbs(f64, cos(-1.5), 0.070737, epsilon));
156		try expect(math.approxEqAbs(f64, cos(37.45), 0.969132, epsilon));
157		try expect(math.approxEqAbs(f64, cos(89.123), 0.40080, epsilon));
	201	test "cos32.normal" {
	202	const epsilon = math.floatEps(f32);
	203	try expectApproxEqAbs(@as(f32, 1.0), cosf(0.0), epsilon);
	204	try expectApproxEqAbs(@as(f32, 0.9800666), cosf(0.2), epsilon);
	205	try expectApproxEqAbs(@as(f32, 0.6276231), cosf(0.8923), epsilon);
	206	try expectApproxEqAbs(@as(f32, 0.0707372), cosf(1.5), epsilon);
	207	try expectApproxEqAbs(@as(f32, 0.0707372), cosf(-1.5), epsilon);
	208	try expectApproxEqAbs(@as(f32, 0.96913195), cosf(37.45), epsilon);
	209	try expectApproxEqAbs(@as(f32, 0.40079966), cosf(89.123), epsilon);
158	210	}
159	211
160	212	test "cos32.special" {
161		try expect(math.isNan(cosf(math.inf(f32))));
162		try expect(math.isNan(cosf(-math.inf(f32))));
163		try expect(math.isNan(cosf(math.nan(f32))));
	213	try testCosSpecial(f32);
	214	}
	215
	216	test "cos64.normal" {
	217	const epsilon = math.floatEps(f64);
	218	try expectApproxEqAbs(@as(f64, 1.0), cos(0.0), epsilon);
	219	try expectApproxEqAbs(@as(f64, 0.9800665778412416), cos(0.2), epsilon);
	220	try expectApproxEqAbs(@as(f64, 0.6276230983360804), cos(0.8923), epsilon);
	221	try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos(1.5), epsilon);
	222	try expectApproxEqAbs(@as(f64, 0.0707372016677029), cos(-1.5), epsilon);
	223	try expectApproxEqAbs(@as(f64, 0.9691317730707778), cos(37.45), epsilon);
	224	try expectApproxEqAbs(@as(f64, 0.4008006809354791), cos(89.123), epsilon);
164	225	}
165	226
166	227	test "cos64.special" {
167		try expect(math.isNan(cos(math.inf(f64))));
168		try expect(math.isNan(cos(-math.inf(f64))));
169		try expect(math.isNan(cos(math.nan(f64))));
	228	try testCosSpecial(f64);
	229	}
	230
	231	test "cos80.normal" {
	232	const epsilon = math.floatEps(f80);
	233	try expectApproxEqAbs(@as(f80, 1.0), cosx(0.0), epsilon);
	234	try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), cosx(0.2), epsilon);
	235	try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), cosx(0.8923), epsilon);
	236	try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cosx(1.5), epsilon);
	237	try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), cosx(-1.5), epsilon);
	238	try expectApproxEqAbs(@as(f80, 0.9691317730707771246), cosx(37.45), epsilon);
	239	try expectApproxEqAbs(@as(f80, 0.4008006809354834001), cosx(89.123), epsilon);
	240	}
	241
	242	test "cos80.special" {
	243	try testCosSpecial(f80);
	244	}
	245
	246	test "cos128.normal" {
	247	const epsilon = math.floatEps(f128);
	248	try expectApproxEqAbs(@as(f128, 1.0), cosq(0.0), epsilon);
	249	try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), cosq(0.2), epsilon);
	250	try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), cosq(0.8923), epsilon);
	251	try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cosq(1.5), epsilon);
	252	try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), cosq(-1.5), epsilon);
	253	try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), cosq(37.45), epsilon);
	254	try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), cosq(89.123), epsilon);
	255	}
	256
	257	test "cos128.special" {
	258	try testCosSpecial(f128);
170	259	}

lib/compiler_rt/long_double.zig created+37

...	...	@@ -0,0 +1,37 @@
	1	//! Utilities for dealing with the `long double` type (`f80` or `f128`)
	2
	3	const std = @import("std");
	4
	5	pub const U80 = std.meta.Int(.unsigned, 80);
	6
	7	/// Returns the sign + exponent bits of a `long double`
	8	pub fn signExponent(x: anytype) u16 {
	9	const T = @TypeOf(x);
	10	switch (T) {
	11	f80 => {
	12	const bits: U80 = @bitCast(x);
	13	return @intCast(bits >> 64);
	14	},
	15	f128 => {
	16	const bits: u128 = @bitCast(x);
	17	return @intCast(bits >> 112);
	18	},
	19	else => @compileError("`signExponent` supports only `f80` and `f128`, got: " ++ @typeName(T)),
	20	}
	21	}
	22
	23	/// Takes the top 16 bits of a `long double`'s mantissa
	24	pub fn mantissaTop(x: anytype) u16 {
	25	const T = @TypeOf(x);
	26	switch (T) {
	27	f80 => {
	28	const bits: U80 = @bitCast(x);
	29	return @intCast((bits >> 48) & 0xFFFF);
	30	},
	31	f128 => {
	32	const bits: u128 = @bitCast(x);
	33	return @intCast((bits >> 96) & 0xFFFF);
	34	},
	35	else => @compileError("`mantissaTop` supports only `f80` and `f128`, got: " ++ @typeName(T)),
	36	}
	37	}

lib/compiler_rt/rem_pio2l.zig created+173

...	...	@@ -0,0 +1,173 @@
	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/__rem_pio2l.c
	5
	6	const std = @import("std");
	7	const math = std.math;
	8
	9	const ld = @import("long_double.zig");
	10	const rem_pio2_large = @import("rem_pio2_large.zig").rem_pio2_large;
	11
	12	pub fn rem_pio2l(comptime T: type, x: T, y: *[2]T) i32 {
	13	const impl = switch (T) {
	14	f80 => struct {
	15	const round1: i8 = 22;
	16	const round2: i8 = 61;
	17	const nx: i8 = 3;
	18	const ny: i8 = 2;
	19
	20	const pio4: T = 0x1.921fb54442d1846ap-1;
	21	// 64 bits of 2/pi
	22	const invpio2: T = 6.36619772367581343076e-01; // 0xa2f9836e4e44152a.0p-64
	23	// first 39 bits of pi/2
	24	const pio2_1: f64 = 1.57079632679597125389e+00; // 0x3FF921FB, 0x54444000
	25	// pi/2 - pio2_1
	26	const pio2_1t: T = -1.07463465549719416346e-12; // -0x973dcb3b399d747f.0p-103
	27	// second 39 bits of pi/2
	28	const pio2_2: f64 = -1.07463465549783099519e-12; // -0x12e7b967674000.0p-92
	29	// pi/2 - (pio2_1+pio2_2)
	30	const pio2_2t: T = 6.36831716351095013979e-25; // 0xc51701b839a25205.0p-144
	31	// pi/2 - (pio2_1+pio2_2+pio2_3)
	32	const pio2_3t: T = -2.75299651904407171810e-37; // -0xbb5bf6c7ddd660ce.0p-185
	33	// third 39 bits of pi/2
	34	const pio2_3: f64 = 6.36831716351370313614e-25; // 0x18a2e037074000.0p-133
	35
	36	fn small(x_val: T) bool {
	37	const se = ld.signExponent(x_val);
	38	const top = ld.mantissaTop(x_val);
	39	const lhs = (@as(u32, se & 0x7fff) << 16) \| top;
	40	const rhs: u32 = ((0x3fff + 25) << 16) \| 0x921f >> 1 \| 0x8000;
	41	return lhs < rhs;
	42	}
	43
	44	fn quobits(v: T) i32 {
	45	const q: i32 = @intFromFloat(v);
	46	return @intCast(@as(u32, @bitCast(q)) & 0x7fffffff);
	47	}
	48	},
	49	f128 => struct {
	50	const round1: i8 = 51;
	51	const round2: i8 = 119;
	52	const nx: i8 = 5;
	53	const ny: i8 = 3;
	54
	55	const pio4: T = 0x1.921fb54442d18469898cc51701b8p-1;
	56	const invpio2: T = 6.3661977236758134307553505349005747e-01;
	57	const pio2_1: T = 1.5707963267948966192292994253909555e+00;
	58	const pio2_1t: T = 2.0222662487959507323996846200947577e-21;
	59	const pio2_2: T = 2.0222662487959507323994779168837751e-21;
	60	const pio2_2t: T = 2.0670321098263988236496903051604844e-43;
	61	const pio2_3: T = 2.0670321098263988236499468110329591e-43;
	62	const pio2_3t: T = -2.5650587247459238361625433492959285e-65;
	63
	64	fn small(x_val: T) bool {
	65	const se = ld.signExponent(x_val);
	66	const top = ld.mantissaTop(x_val);
	67	const lhs = (@as(u32, se & 0x7fff) << 16) \| top;
	68	const rhs: u32 = ((0x3fff + 45) << 16) \| 0x921f;
	69	return lhs < rhs;
	70	}
	71
	72	fn quobits(fn_val: T) i32 {
	73	const q: i64 = @intFromFloat(fn_val);
	74	return @intCast(@as(u64, @bitCast(q)) & 0x7fffffff);
	75	}
	76	},
	77	else => @compileError("rem_pio2l supports only f80 and f128, got: " ++ @typeName(T)),
	78	};
	79
	80	const x_se = ld.signExponent(x);
	81	const ex: i32 = @intCast(x_se & 0x7fff);
	82
	83	if (impl.small(x)) {
	84	// rint(x/(pi/2))
	85	const toint: T = 1.5 / math.floatEps(T);
	86	var fn_ = x * impl.invpio2 + toint - toint;
	87	var n = impl.quobits(fn_);
	88	var r = x - fn_ * @as(T, impl.pio2_1);
	89	var w = fn_ * impl.pio2_1t; // 1st round good to 102/180 bits
	90
	91	// Matters with directed rounding.
	92	if (r - w < -impl.pio4) {
	93	@branchHint(.unlikely);
	94	n -= 1;
	95	fn_ -= 1;
	96	r = x - fn_ * @as(T, impl.pio2_1);
	97	w = fn_ * impl.pio2_1t;
	98	} else if (r - w > impl.pio4) {
	99	@branchHint(.unlikely);
	100	n += 1;
	101	fn_ += 1;
	102	r = x - fn_ * @as(T, impl.pio2_1);
	103	w = fn_ * impl.pio2_1t;
	104	}
	105
	106	y[0] = r - w;
	107
	108	const ey: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff);
	109	if (ex - ey > impl.round1) {
	110	var t = r;
	111	w = fn_ * impl.pio2_2;
	112	r = t - w;
	113	w = fn_ * impl.pio2_2t - ((t - r) - w);
	114	y[0] = r - w;
	115	const ey2: i32 = @intCast(ld.signExponent(y[0]) & 0x7fff);
	116	if (ex - ey2 > impl.round2) {
	117	t = r;
	118	w = fn_ * impl.pio2_3;
	119	r = t - w;
	120	w = fn_ * impl.pio2_3t - ((t - r) - w);
	121	y[0] = r - w;
	122	}
	123	}
	124	y[1] = (r - y[0]) - w;
	125	return n;
	126	}
	127
	128	// all other (large) arguments
	129	if (ex == 0x7fff) { // x is inf or NaN
	130	y[0] = x - x;
	131	y[1] = y[0];
	132	return 0;
	133	}
	134
	135	var z: T = math.scalbn(@abs(x), -math.ilogb(x) + 23);
	136	var tx: [impl.nx]f64 = undefined;
	137	var ty: [impl.ny]f64 = undefined;
	138	var i: usize = 0;
	139
	140	while (i < impl.nx - 1) : (i += 1) {
	141	tx[i] = @floatFromInt(@as(i32, @intFromFloat(z)));
	142	z = (z - @as(T, tx[i])) * 0x1p24;
	143	}
	144
	145	tx[i] = @floatCast(z);
	146	while (tx[i] == 0.0) {
	147	i -= 1;
	148	}
	149
	150	const n = rem_pio2_large(
	151	tx[0..(i + 1)],
	152	ty[0..impl.ny],
	153	ex - 0x3fff - 23,
	154	@intCast(i + 1),
	155	impl.ny,
	156	);
	157	var w: f64 = ty[1];
	158	if (impl.ny == 3) {
	159	w += ty[2];
	160	}
	161	const r = ty[0] + w;
	162	w -= r - ty[0];
	163
	164	if (x_se >> 15 != 0) {
	165	y[0] = -@as(T, r);
	166	y[1] = -@as(T, w);
	167	return -n;
	168	}
	169
	170	y[0] = @as(T, r);
	171	y[1] = @as(T, w);
	172	return n;
	173	}

lib/compiler_rt/sin.zig+140-55

...	...	@@ -3,23 +3,28 @@
3	3	//!
4	4	//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinf.c
5	5	//! https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c
	6	//! https://git.musl-libc.org/cgit/musl/tree/src/math/sinl.c
6	7
7	8	const std = @import("std");
8	9	const math = std.math;
9	10	const mem = std.mem;
10	11	const expect = std.testing.expect;
	12	const expectApproxEqAbs = std.testing.expectApproxEqAbs;
11	13
12	14	const compiler_rt = @import("../compiler_rt.zig");
13	15	const symbol = @import("../compiler_rt.zig").symbol;
14	16	const trig = @import("trig.zig");
15	17	const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
16	18	const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
	19	const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
	20	const ld = @import("long_double.zig");
17	21
18	22	comptime {
19		symbol(&__sinh, "__sinh");
	23	symbol(&sinh, "__sinh");
	24	symbol(&sinl, "__sinl");
20	25	symbol(&sinf, "sinf");
21	26	symbol(&sin, "sin");
22		symbol(&__sinx, "__sinx");
	27	symbol(&sinx, "__sinx");
23	28	if (compiler_rt.want_ppc_abi) {
24	29	symbol(&sinq, "sinf128");
25	30	}
...	...	@@ -27,7 +32,7 @@ comptime {
27	32	symbol(&sinl, "sinl");
28	33	}
29	34
30		pub fn __sinh(x: f16) callconv(.c) f16 {
	35	pub fn sinh(x: f16) callconv(.c) f16 {
31	36	// TODO: more efficient implementation
32	37	return @floatCast(sinf(x));
33	38	}
...	...	@@ -55,27 +60,27 @@ pub fn sinf(x: f32) callconv(.c) f32 {
55	60	}
56	61	return x;
57	62	}
58		return trig.__sindf(x);
	63	return trig.sindf(x);
59	64	}
60	65	if (ix <= 0x407b53d1) { // \|x\| ~<= 5*pi/4
61	66	if (ix <= 0x4016cbe3) { // \|x\| ~<= 3pi/4
62	67	if (sign) {
63		return -trig.__cosdf(x + s1pio2);
	68	return -trig.cosdf(x + s1pio2);
64	69	} else {
65		return trig.__cosdf(x - s1pio2);
	70	return trig.cosdf(x - s1pio2);
66	71	}
67	72	}
68		return trig.__sindf(if (sign) -(x + s2pio2) else -(x - s2pio2));
	73	return trig.sindf(if (sign) -(x + s2pio2) else -(x - s2pio2));
69	74	}
70	75	if (ix <= 0x40e231d5) { // \|x\| ~<= 9*pi/4
71	76	if (ix <= 0x40afeddf) { // \|x\| ~<= 7*pi/4
72	77	if (sign) {
73		return trig.__cosdf(x + s3pio2);
	78	return trig.cosdf(x + s3pio2);
74	79	} else {
75		return -trig.__cosdf(x - s3pio2);
	80	return -trig.cosdf(x - s3pio2);
76	81	}
77	82	}
78		return trig.__sindf(if (sign) x + s4pio2 else x - s4pio2);
	83	return trig.sindf(if (sign) x + s4pio2 else x - s4pio2);
79	84	}
80	85
81	86	// sin(Inf or NaN) is NaN
...	...	@@ -86,10 +91,10 @@ pub fn sinf(x: f32) callconv(.c) f32 {
86	91	var y: f64 = undefined;
87	92	const n = rem_pio2f(x, &y);
88	93	return switch (n & 3) {
89		0 => trig.__sindf(y),
90		1 => trig.__cosdf(y),
91		2 => trig.__sindf(-y),
92		else => -trig.__cosdf(y),
	94	0 => trig.sindf(y),
	95	1 => trig.cosdf(y),
	96	2 => trig.sindf(-y),
	97	else => -trig.cosdf(y),
93	98	};
94	99	}
95	100
...	...	@@ -110,7 +115,7 @@ pub fn sin(x: f64) callconv(.c) f64 {
110	115	}
111	116	return x;
112	117	}
113		return trig.__sin(x, 0.0, 0);
	118	return trig.sin(x, 0.0, 0);
114	119	}
115	120
116	121	// sin(Inf or NaN) is NaN
...	...	@@ -121,72 +126,152 @@ pub fn sin(x: f64) callconv(.c) f64 {
121	126	var y: [2]f64 = undefined;
122	127	const n = rem_pio2(x, &y);
123	128	return switch (n & 3) {
124		0 => trig.__sin(y[0], y[1], 1),
125		1 => trig.__cos(y[0], y[1]),
126		2 => -trig.__sin(y[0], y[1], 1),
127		else => -trig.__cos(y[0], y[1]),
	129	0 => trig.sin(y[0], y[1], 1),
	130	1 => trig.cos(y[0], y[1]),
	131	2 => -trig.sin(y[0], y[1], 1),
	132	else => -trig.cos(y[0], y[1]),
128	133	};
129	134	}
130	135
131		pub fn __sinx(x: f80) callconv(.c) f80 {
132		// TODO: more efficient implementation
133		return @floatCast(sinq(x));
	136	fn sinx(x: f80) callconv(.c) f80 {
	137	const se = ld.signExponent(x) & 0x7fff;
	138	if (se == 0x7fff) {
	139	return x - x;
	140	}
	141
	142	if (@abs(x) < trig.pi_4) {
	143	if (se < 0x3fff - (math.floatMantissaBits(f80) / 2)) {
	144	// raise inexact if x!=0 and underflow if subnormal
	145	if (compiler_rt.want_float_exceptions) {
	146	mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
	147	}
	148	return x;
	149	}
	150	return trig.sinx(x, 0.0, 0);
	151	}
	152
	153	var y: [2]f80 = undefined;
	154	const n = rem_pio2l(f80, x, &y);
	155	return switch (n & 3) {
	156	0 => trig.sinx(y[0], y[1], 1),
	157	1 => trig.cosx(y[0], y[1]),
	158	2 => -trig.sinx(y[0], y[1], 1),
	159	else => -trig.cosx(y[0], y[1]),
	160	};
134	161	}
135	162
136	163	pub fn sinq(x: f128) callconv(.c) f128 {
137		// TODO: more correct implementation
138		return sin(@floatCast(x));
	164	const se = ld.signExponent(x) & 0x7fff;
	165	if (se == 0x7fff) {
	166	return x - x;
	167	}
	168
	169	if (@abs(x) < trig.pi_4) {
	170	if (se < 0x3fff - (math.floatMantissaBits(f128) / 2)) {
	171	// raise inexact if x!=0 and underflow if subnormal
	172	if (compiler_rt.want_float_exceptions) {
	173	mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
	174	}
	175	return x;
	176	}
	177	return trig.sinq(x, 0.0, 0);
	178	}
	179
	180	var y: [2]f128 = undefined;
	181	const n = rem_pio2l(f128, x, &y);
	182	return switch (n & 3) {
	183	0 => trig.sinq(y[0], y[1], 1),
	184	1 => trig.cosq(y[0], y[1]),
	185	2 => -trig.sinq(y[0], y[1], 1),
	186	else => -trig.cosq(y[0], y[1]),
	187	};
139	188	}
140	189
141	190	pub fn sinl(x: c_longdouble) callconv(.c) c_longdouble {
142	191	switch (@typeInfo(c_longdouble).float.bits) {
143		16 => return __sinh(x),
	192	16 => return sinh(x),
144	193	32 => return sinf(x),
145	194	64 => return sin(x),
146		80 => return __sinx(x),
	195	80 => return sinx(x),
147	196	128 => return sinq(x),
148	197	else => @compileError("unreachable"),
149	198	}
150	199	}
151	200
152		test "sin32" {
153		const epsilon = 0.00001;
	201	fn testSinSpecial(comptime T: type) !void {
	202	const f = switch (T) {
	203	f32 => sinf,
	204	f64 => sin,
	205	f80 => sinx,
	206	f128 => sinq,
	207	else => @compileError("unimplemented"),
	208	};
154	209
155		try expect(math.approxEqAbs(f32, sinf(0.0), 0.0, epsilon));
156		try expect(math.approxEqAbs(f32, sinf(0.2), 0.198669, epsilon));
157		try expect(math.approxEqAbs(f32, sinf(0.8923), 0.778517, epsilon));
158		try expect(math.approxEqAbs(f32, sinf(1.5), 0.997495, epsilon));
159		try expect(math.approxEqAbs(f32, sinf(-1.5), -0.997495, epsilon));
160		try expect(math.approxEqAbs(f32, sinf(37.45), -0.246544, epsilon));
161		try expect(math.approxEqAbs(f32, sinf(89.123), 0.916166, epsilon));
	210	try expect(math.isPositiveZero(f(0.0)));
	211	try expect(math.isNegativeZero(f(-0.0)));
	212	try expect(math.isNan(f(math.inf(T))));
	213	try expect(math.isNan(f(-math.inf(T))));
	214	try expect(math.isNan(f(math.nan(T))));
162	215	}
163	216
164		test "sin64" {
165		const epsilon = 0.000001;
166
167		try expect(math.approxEqAbs(f64, sin(0.0), 0.0, epsilon));
168		try expect(math.approxEqAbs(f64, sin(0.2), 0.198669, epsilon));
169		try expect(math.approxEqAbs(f64, sin(0.8923), 0.778517, epsilon));
170		try expect(math.approxEqAbs(f64, sin(1.5), 0.997495, epsilon));
171		try expect(math.approxEqAbs(f64, sin(-1.5), -0.997495, epsilon));
172		try expect(math.approxEqAbs(f64, sin(37.45), -0.246543, epsilon));
173		try expect(math.approxEqAbs(f64, sin(89.123), 0.916166, epsilon));
	217	test "sin32.normal" {
	218	const epsilon = math.floatEps(f32);
	219	try expectApproxEqAbs(@as(f32, 0.0), sinf(0.0), epsilon);
	220	try expectApproxEqAbs(@as(f32, 0.19866933), sinf(0.2), epsilon);
	221	try expectApproxEqAbs(@as(f32, 0.77851737), sinf(0.8923), epsilon);
	222	try expectApproxEqAbs(@as(f32, 0.997495), sinf(1.5), epsilon);
	223	try expectApproxEqAbs(@as(f32, -0.997495), sinf(-1.5), epsilon);
	224	try expectApproxEqAbs(@as(f32, -0.24654257), sinf(37.45), epsilon);
	225	try expectApproxEqAbs(@as(f32, 0.9161657), sinf(89.123), epsilon);
174	226	}
175	227
176	228	test "sin32.special" {
177		try expect(sinf(0.0) == 0.0);
178		try expect(sinf(-0.0) == -0.0);
179		try expect(math.isNan(sinf(math.inf(f32))));
180		try expect(math.isNan(sinf(-math.inf(f32))));
181		try expect(math.isNan(sinf(math.nan(f32))));
	229	try testSinSpecial(f32);
	230	}
	231
	232	test "sin64.normal" {
	233	const epsilon = math.floatEps(f64);
	234	try expectApproxEqAbs(@as(f64, 0.0), sin(0.0), epsilon);
	235	try expectApproxEqAbs(@as(f64, 0.19866933079506122), sin(0.2), epsilon);
	236	try expectApproxEqAbs(@as(f64, 0.7785173385577349), sin(0.8923), epsilon);
	237	try expectApproxEqAbs(@as(f64, 0.9974949866040544), sin(1.5), epsilon);
	238	try expectApproxEqAbs(@as(f64, -0.9974949866040544), sin(-1.5), epsilon);
	239	try expectApproxEqAbs(@as(f64, -0.24654331551411082), sin(37.45), epsilon);
	240	try expectApproxEqAbs(@as(f64, 0.9161652766622714), sin(89.123), epsilon);
182	241	}
183	242
184	243	test "sin64.special" {
185		try expect(sin(0.0) == 0.0);
186		try expect(sin(-0.0) == -0.0);
187		try expect(math.isNan(sin(math.inf(f64))));
188		try expect(math.isNan(sin(-math.inf(f64))));
189		try expect(math.isNan(sin(math.nan(f64))));
	244	try testSinSpecial(f64);
	245	}
	246
	247	test "sin80.normal" {
	248	const epsilon = math.floatEps(f80);
	249	try expectApproxEqAbs(@as(f80, 0.0), sinx(0.0), epsilon);
	250	try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), sinx(0.2), epsilon);
	251	try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), sinx(0.8923), epsilon);
	252	try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), sinx(1.5), epsilon);
	253	try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), sinx(-1.5), epsilon);
	254	try expectApproxEqAbs(@as(f80, -0.24654331551411356504), sinx(37.45), epsilon);
	255	try expectApproxEqAbs(@as(f80, 0.91616527666226951006), sinx(89.123), epsilon);
	256	}
	257
	258	test "sin80.special" {
	259	try testSinSpecial(f80);
	260	}
	261
	262	test "sin128.normal" {
	263	const epsilon = math.floatEps(f128);
	264	try expectApproxEqAbs(@as(f128, 0.0), sinq(0.0), epsilon);
	265	try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), sinq(0.2), epsilon);
	266	try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), sinq(0.8923), epsilon);
	267	try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), sinq(1.5), epsilon);
	268	try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), sinq(-1.5), epsilon);
	269	try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), sinq(37.45), epsilon);
	270	try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), sinq(89.123), epsilon);
	271	}
	272
	273	test "sin128.special" {
	274	try testSinSpecial(f128);
190	275	}
191	276
192	277	test "sin32 #9901" {

lib/compiler_rt/sincos.zig+279-72

...	...	@@ -3,17 +3,21 @@ const builtin = @import("builtin");
3	3	const arch = builtin.cpu.arch;
4	4	const math = std.math;
5	5	const mem = std.mem;
	6	const expect = std.testing.expect;
	7	const expectApproxEqAbs = std.testing.expectApproxEqAbs;
6	8	const trig = @import("trig.zig");
7	9	const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
8	10	const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
	11	const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
	12	const ld = @import("long_double.zig");
9	13	const compiler_rt = @import("../compiler_rt.zig");
10	14	const symbol = compiler_rt.symbol;
11	15
12	16	comptime {
13		symbol(&__sincosh, "__sincosh");
	17	symbol(&sincosh, "__sincosh");
14	18	symbol(&sincosf, "sincosf");
15	19	symbol(&sincos, "sincos");
16		symbol(&__sincosx, "__sincosx");
	20	symbol(&sincosx, "__sincosx");
17	21	if (compiler_rt.want_ppc_abi) {
18	22	symbol(&sincosq, "sincosf128");
19	23	}
...	...	@@ -21,7 +25,7 @@ comptime {
21	25	symbol(&sincosl, "sincosl");
22	26	}
23	27
24		pub fn __sincosh(x: f16, r_sin: f16, r_cos: f16) callconv(.c) void {
	28	pub fn sincosh(x: f16, r_sin: f16, r_cos: f16) callconv(.c) void {
25	29	// TODO: more efficient implementation
26	30	var big_sin: f32 = undefined;
27	31	var big_cos: f32 = undefined;
...	...	@@ -56,8 +60,8 @@ pub fn sincosf(x: f32, r_sin: f32, r_cos: f32) callconv(.c) void {
56	60	r_cos.* = 1.0;
57	61	return;
58	62	}
59		r_sin.* = trig.__sindf(x);
60		r_cos.* = trig.__cosdf(x);
	63	r_sin.* = trig.sindf(x);
	64	r_cos.* = trig.cosdf(x);
61	65	return;
62	66	}
63	67
...	...	@@ -66,17 +70,17 @@ pub fn sincosf(x: f32, r_sin: f32, r_cos: f32) callconv(.c) void {
66	70	// \|x\| ~<= 3pi/4
67	71	if (ix <= 0x4016cbe3) {
68	72	if (sign) {
69		r_sin.* = -trig.__cosdf(x + sc1pio2);
70		r_cos.* = trig.__sindf(x + sc1pio2);
	73	r_sin.* = -trig.cosdf(x + sc1pio2);
	74	r_cos.* = trig.sindf(x + sc1pio2);
71	75	} else {
72		r_sin.* = trig.__cosdf(sc1pio2 - x);
73		r_cos.* = trig.__sindf(sc1pio2 - x);
	76	r_sin.* = trig.cosdf(sc1pio2 - x);
	77	r_cos.* = trig.sindf(sc1pio2 - x);
74	78	}
75	79	return;
76	80	}
77	81	// -sin(x+c) is not correct if x+c could be 0: -0 vs +0
78		r_sin.* = -trig.__sindf(if (sign) x + sc2pio2 else x - sc2pio2);
79		r_cos.* = -trig.__cosdf(if (sign) x + sc2pio2 else x - sc2pio2);
	82	r_sin.* = -trig.sindf(if (sign) x + sc2pio2 else x - sc2pio2);
	83	r_cos.* = -trig.cosdf(if (sign) x + sc2pio2 else x - sc2pio2);
80	84	return;
81	85	}
82	86
...	...	@@ -85,16 +89,16 @@ pub fn sincosf(x: f32, r_sin: f32, r_cos: f32) callconv(.c) void {
85	89	// \|x\| ~<= 7*pi/4
86	90	if (ix <= 0x40afeddf) {
87	91	if (sign) {
88		r_sin.* = trig.__cosdf(x + sc3pio2);
89		r_cos.* = -trig.__sindf(x + sc3pio2);
	92	r_sin.* = trig.cosdf(x + sc3pio2);
	93	r_cos.* = -trig.sindf(x + sc3pio2);
90	94	} else {
91		r_sin.* = -trig.__cosdf(x - sc3pio2);
92		r_cos.* = trig.__sindf(x - sc3pio2);
	95	r_sin.* = -trig.cosdf(x - sc3pio2);
	96	r_cos.* = trig.sindf(x - sc3pio2);
93	97	}
94	98	return;
95	99	}
96		r_sin.* = trig.__sindf(if (sign) x + sc4pio2 else x - sc4pio2);
97		r_cos.* = trig.__cosdf(if (sign) x + sc4pio2 else x - sc4pio2);
	100	r_sin.* = trig.sindf(if (sign) x + sc4pio2 else x - sc4pio2);
	101	r_cos.* = trig.cosdf(if (sign) x + sc4pio2 else x - sc4pio2);
98	102	return;
99	103	}
100	104
...	...	@@ -109,8 +113,8 @@ pub fn sincosf(x: f32, r_sin: f32, r_cos: f32) callconv(.c) void {
109	113	// general argument reduction needed
110	114	var y: f64 = undefined;
111	115	const n = rem_pio2f(x, &y);
112		const s = trig.__sindf(y);
113		const c = trig.__cosdf(y);
	116	const s = trig.sindf(y);
	117	const c = trig.cosdf(y);
114	118	switch (n & 3) {
115	119	0 => {
116	120	r_sin.* = s;
...	...	@@ -150,8 +154,8 @@ pub fn sincos(x: f64, r_sin: f64, r_cos: f64) callconv(.c) void {
150	154	r_cos.* = 1.0;
151	155	return;
152	156	}
153		r_sin.* = trig.__sin(x, 0.0, 0);
154		r_cos.* = trig.__cos(x, 0.0);
	157	r_sin.* = trig.sin(x, 0.0, 0);
	158	r_cos.* = trig.cos(x, 0.0);
155	159	return;
156	160	}
157	161
...	...	@@ -166,8 +170,8 @@ pub fn sincos(x: f64, r_sin: f64, r_cos: f64) callconv(.c) void {
166	170	// argument reduction needed
167	171	var y: [2]f64 = undefined;
168	172	const n = rem_pio2(x, &y);
169		const s = trig.__sin(y[0], y[1], 1);
170		const c = trig.__cos(y[0], y[1]);
	173	const s = trig.sin(y[0], y[1], 1);
	174	const c = trig.cos(y[0], y[1]);
171	175	switch (n & 3) {
172	176	0 => {
173	177	r_sin.* = s;
...	...	@@ -188,50 +192,57 @@ pub fn sincos(x: f64, r_sin: f64, r_cos: f64) callconv(.c) void {
188	192	}
189	193	}
190	194
191		pub fn __sincosx(x: f80, r_sin: f80, r_cos: f80) callconv(.c) void {
192		// TODO: more efficient implementation
193		//return sincos_generic(f80, x, r_sin, r_cos);
194		var big_sin: f128 = undefined;
195		var big_cos: f128 = undefined;
196		sincosq(x, &big_sin, &big_cos);
197		r_sin.* = @as(f80, @floatCast(big_sin));
198		r_cos.* = @as(f80, @floatCast(big_cos));
199		}
	195	pub fn sincosx(x: f80, r_sin: f80, r_cos: f80) callconv(.c) void {
	196	const se = ld.signExponent(x) & 0x7fff;
	197	if (se == 0x7fff) {
	198	const result = x - x;
	199	r_sin.* = result;
	200	r_cos.* = result;
	201	return;
	202	}
200	203
201		pub fn sincosq(x: f128, r_sin: f128, r_cos: f128) callconv(.c) void {
202		// TODO: more correct implementation
203		//return sincos_generic(f128, x, r_sin, r_cos);
204		var small_sin: f64 = undefined;
205		var small_cos: f64 = undefined;
206		sincos(@as(f64, @floatCast(x)), &small_sin, &small_cos);
207		r_sin.* = small_sin;
208		r_cos.* = small_cos;
209		}
	204	if (@abs(x) < trig.pi_4) {
	205	if (se < 0x3fff - math.floatMantissaBits(f80)) {
	206	// raise underflow if subnormal
	207	if (compiler_rt.want_float_exceptions and se == 0) {
	208	mem.doNotOptimizeAway(x * 0x1p-120);
	209	}
	210	r_sin.* = x;
	211	// raise inexact if x!=0
	212	r_cos.* = 1.0 + x;
	213	return;
	214	}
	215	r_sin.* = trig.sinx(x, 0.0, 0);
	216	r_cos.* = trig.cosx(x, 0.0);
	217	return;
	218	}
210	219
211		pub fn sincosl(x: c_longdouble, r_sin: c_longdouble, r_cos: c_longdouble) callconv(.c) void {
212		switch (@typeInfo(c_longdouble).float.bits) {
213		16 => return __sincosh(x, r_sin, r_cos),
214		32 => return sincosf(x, r_sin, r_cos),
215		64 => return sincos(x, r_sin, r_cos),
216		80 => return __sincosx(x, r_sin, r_cos),
217		128 => return sincosq(x, r_sin, r_cos),
218		else => @compileError("unreachable"),
	220	var y: [2]f80 = undefined;
	221	const n = rem_pio2l(f80, x, &y);
	222	const s = trig.sinx(y[0], y[1], 1);
	223	const c = trig.cosx(y[0], y[1]);
	224	switch (n & 3) {
	225	0 => {
	226	r_sin.* = s;
	227	r_cos.* = c;
	228	},
	229	1 => {
	230	r_sin.* = c;
	231	r_cos.* = -s;
	232	},
	233	2 => {
	234	r_sin.* = -s;
	235	r_cos.* = -c;
	236	},
	237	else => {
	238	r_sin.* = -c;
	239	r_cos.* = s;
	240	},
219	241	}
220	242	}
221	243
222		pub const rem_pio2_generic = @compileError("TODO");
223
224		/// Ported from musl sincosl.c. Needs the following dependencies to be complete:
225		/// * rem_pio2_generic ported from __rem_pio2l.c
226		/// * trig.sin_generic ported from __sinl.c
227		/// * trig.cos_generic ported from __cosl.c
228		inline fn sincos_generic(comptime F: type, x: F, r_sin: F, r_cos: F) void {
229		const sc1pio4: F = 1.0 * math.pi / 4.0;
230		const bits = @typeInfo(F).float.bits;
231		const I = std.meta.Int(.unsigned, bits);
232		const ix = @as(I, @bitCast(x)) & (math.maxInt(I) >> 1);
233		const se: u16 = @truncate(ix >> (bits - 16));
234
	244	pub fn sincosq(x: f128, r_sin: f128, r_cos: f128) callconv(.c) void {
	245	const se = ld.signExponent(x) & 0x7fff;
235	246	if (se == 0x7fff) {
236	247	const result = x - x;
237	248	r_sin.* = result;
...	...	@@ -239,26 +250,26 @@ inline fn sincos_generic(comptime F: type, x: F, r_sin: F, r_cos: F) void {
239	250	return;
240	251	}
241	252
242		if (@as(F, @bitCast(ix)) < sc1pio4) {
243		if (se < 0x3fff - math.floatFractionalBits(F) - 1) {
	253	if (@abs(x) < trig.pi_4) {
	254	if (se < 0x3fff - math.floatMantissaBits(f128)) {
244	255	// raise underflow if subnormal
245		if (se == 0) {
246		if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(x * 0x1p-120);
	256	if (compiler_rt.want_float_exceptions and se == 0) {
	257	mem.doNotOptimizeAway(x * 0x1p-120);
247	258	}
248	259	r_sin.* = x;
249	260	// raise inexact if x!=0
250	261	r_cos.* = 1.0 + x;
251	262	return;
252	263	}
253		r_sin.* = trig.sin_generic(F, x, 0, 0);
254		r_cos.* = trig.cos_generic(F, x, 0);
	264	r_sin.* = trig.sinq(x, 0.0, 0);
	265	r_cos.* = trig.cosq(x, 0.0);
255	266	return;
256	267	}
257	268
258		var y: [2]F = undefined;
259		const n = rem_pio2_generic(F, x, &y);
260		const s = trig.sin_generic(F, y[0], y[1], 1);
261		const c = trig.cos_generic(F, y[0], y[1]);
	269	var y: [2]f128 = undefined;
	270	const n = rem_pio2l(f128, x, &y);
	271	const s = trig.sinq(y[0], y[1], 1);
	272	const c = trig.cosq(y[0], y[1]);
262	273	switch (n & 3) {
263	274	0 => {
264	275	r_sin.* = s;
...	...	@@ -278,3 +289,199 @@ inline fn sincos_generic(comptime F: type, x: F, r_sin: F, r_cos: F) void {
278	289	},
279	290	}
280	291	}
	292
	293	pub fn sincosl(x: c_longdouble, r_sin: c_longdouble, r_cos: c_longdouble) callconv(.c) void {
	294	switch (@typeInfo(c_longdouble).float.bits) {
	295	16 => return sincosh(x, r_sin, r_cos),
	296	32 => return sincosf(x, r_sin, r_cos),
	297	64 => return sincos(x, r_sin, r_cos),
	298	80 => return sincosx(x, r_sin, r_cos),
	299	128 => return sincosq(x, r_sin, r_cos),
	300	else => @compileError("unreachable"),
	301	}
	302	}
	303
	304	fn testSincosSpecial(comptime T: type) !void {
	305	const f = switch (T) {
	306	f32 => sincosf,
	307	f64 => sincos,
	308	f80 => sincosx,
	309	f128 => sincosq,
	310	else => @compileError("unimplemented"),
	311	};
	312
	313	var s: T = undefined;
	314	var c: T = undefined;
	315
	316	f(0.0, &s, &c);
	317	try expect(math.isPositiveZero(s));
	318	try expect(c == 1.0);
	319
	320	f(-0.0, &s, &c);
	321	try expect(math.isNegativeZero(s));
	322	try expect(c == 1.0);
	323
	324	f(math.inf(T), &s, &c);
	325	try expect(math.isNan(s));
	326	try expect(math.isNan(c));
	327
	328	f(-math.inf(T), &s, &c);
	329	try expect(math.isNan(s));
	330	try expect(math.isNan(c));
	331
	332	f(math.nan(T), &s, &c);
	333	try expect(math.isNan(s));
	334	try expect(math.isNan(c));
	335	}
	336
	337	test "sincos32.normal" {
	338	const epsilon = math.floatEps(f32);
	339	var s: f32 = undefined;
	340	var c: f32 = undefined;
	341
	342	sincosf(0.0, &s, &c);
	343	try expectApproxEqAbs(@as(f32, 0.0), s, epsilon);
	344	try expectApproxEqAbs(@as(f32, 1.0), c, epsilon);
	345
	346	sincosf(0.2, &s, &c);
	347	try expectApproxEqAbs(@as(f32, 0.19866933), s, epsilon);
	348	try expectApproxEqAbs(@as(f32, 0.9800666), c, epsilon);
	349
	350	sincosf(0.8923, &s, &c);
	351	try expectApproxEqAbs(@as(f32, 0.77851737), s, epsilon);
	352	try expectApproxEqAbs(@as(f32, 0.6276231), c, epsilon);
	353
	354	sincosf(1.5, &s, &c);
	355	try expectApproxEqAbs(@as(f32, 0.997495), s, epsilon);
	356	try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
	357
	358	sincosf(-1.5, &s, &c);
	359	try expectApproxEqAbs(@as(f32, -0.997495), s, epsilon);
	360	try expectApproxEqAbs(@as(f32, 0.0707372), c, epsilon);
	361
	362	sincosf(37.45, &s, &c);
	363	try expectApproxEqAbs(@as(f32, -0.24654257), s, epsilon);
	364	try expectApproxEqAbs(@as(f32, 0.96913195), c, epsilon);
	365
	366	sincosf(89.123, &s, &c);
	367	try expectApproxEqAbs(@as(f32, 0.9161657), s, epsilon);
	368	try expectApproxEqAbs(@as(f32, 0.40079966), c, epsilon);
	369	}
	370
	371	test "sincos32.special" {
	372	try testSincosSpecial(f32);
	373	}
	374
	375	test "sincos64.normal" {
	376	const epsilon = math.floatEps(f64);
	377	var s: f64 = undefined;
	378	var c: f64 = undefined;
	379
	380	sincos(0.0, &s, &c);
	381	try expectApproxEqAbs(@as(f64, 0.0), s, epsilon);
	382	try expectApproxEqAbs(@as(f64, 1.0), c, epsilon);
	383
	384	sincos(0.2, &s, &c);
	385	try expectApproxEqAbs(@as(f64, 0.19866933079506122), s, epsilon);
	386	try expectApproxEqAbs(@as(f64, 0.9800665778412416), c, epsilon);
	387
	388	sincos(0.8923, &s, &c);
	389	try expectApproxEqAbs(@as(f64, 0.7785173385577349), s, epsilon);
	390	try expectApproxEqAbs(@as(f64, 0.6276230983360804), c, epsilon);
	391
	392	sincos(1.5, &s, &c);
	393	try expectApproxEqAbs(@as(f64, 0.9974949866040544), s, epsilon);
	394	try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
	395
	396	sincos(-1.5, &s, &c);
	397	try expectApproxEqAbs(@as(f64, -0.9974949866040544), s, epsilon);
	398	try expectApproxEqAbs(@as(f64, 0.0707372016677029), c, epsilon);
	399
	400	sincos(37.45, &s, &c);
	401	try expectApproxEqAbs(@as(f64, -0.24654331551411082), s, epsilon);
	402	try expectApproxEqAbs(@as(f64, 0.9691317730707778), c, epsilon);
	403
	404	sincos(89.123, &s, &c);
	405	try expectApproxEqAbs(@as(f64, 0.9161652766622714), s, epsilon);
	406	try expectApproxEqAbs(@as(f64, 0.4008006809354791), c, epsilon);
	407	}
	408
	409	test "sincos64.special" {
	410	try testSincosSpecial(f64);
	411	}
	412
	413	test "sincos80.normal" {
	414	const epsilon = math.floatEps(f80);
	415	var s: f80 = undefined;
	416	var c: f80 = undefined;
	417
	418	sincosx(0.0, &s, &c);
	419	try expectApproxEqAbs(@as(f80, 0.0), s, epsilon);
	420	try expectApproxEqAbs(@as(f80, 1.0), c, epsilon);
	421
	422	sincosx(0.2, &s, &c);
	423	try expectApproxEqAbs(@as(f80, 0.19866933079506121545941262711838975), s, epsilon);
	424	try expectApproxEqAbs(@as(f80, 0.98006657784124163112419651674816888), c, epsilon);
	425
	426	sincosx(0.8923, &s, &c);
	427	try expectApproxEqAbs(@as(f80, 0.77851733855773487830689285621486050), s, epsilon);
	428	try expectApproxEqAbs(@as(f80, 0.62762309833608037003563995939286067), c, epsilon);
	429
	430	sincosx(1.5, &s, &c);
	431	try expectApproxEqAbs(@as(f80, 0.99749498660405443094172337114148732), s, epsilon);
	432	try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
	433
	434	sincosx(-1.5, &s, &c);
	435	try expectApproxEqAbs(@as(f80, -0.99749498660405443094172337114148732), s, epsilon);
	436	try expectApproxEqAbs(@as(f80, 0.070737201667702910088189851434268747), c, epsilon);
	437
	438	sincosx(37.45, &s, &c);
	439	try expectApproxEqAbs(@as(f80, -0.24654331551411356504), s, epsilon);
	440	try expectApproxEqAbs(@as(f80, 0.9691317730707771246), c, epsilon);
	441
	442	sincosx(89.123, &s, &c);
	443	try expectApproxEqAbs(@as(f80, 0.91616527666226951006), s, epsilon);
	444	try expectApproxEqAbs(@as(f80, 0.4008006809354834001), c, epsilon);
	445	}
	446
	447	test "sincos80.special" {
	448	try testSincosSpecial(f80);
	449	}
	450
	451	test "sincos128.normal" {
	452	const epsilon = math.floatEps(f128);
	453	var s: f128 = undefined;
	454	var c: f128 = undefined;
	455
	456	sincosq(0.0, &s, &c);
	457	try expectApproxEqAbs(@as(f128, 0.0), s, epsilon);
	458	try expectApproxEqAbs(@as(f128, 1.0), c, epsilon);
	459
	460	sincosq(0.2, &s, &c);
	461	try expectApproxEqAbs(@as(f128, 0.19866933079506121545941262711838975), s, epsilon);
	462	try expectApproxEqAbs(@as(f128, 0.98006657784124163112419651674816888), c, epsilon);
	463
	464	sincosq(0.8923, &s, &c);
	465	try expectApproxEqAbs(@as(f128, 0.77851733855773487830689285621486050), s, epsilon);
	466	try expectApproxEqAbs(@as(f128, 0.62762309833608037003563995939286067), c, epsilon);
	467
	468	sincosq(1.5, &s, &c);
	469	try expectApproxEqAbs(@as(f128, 0.99749498660405443094172337114148732), s, epsilon);
	470	try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
	471
	472	sincosq(-1.5, &s, &c);
	473	try expectApproxEqAbs(@as(f128, -0.99749498660405443094172337114148732), s, epsilon);
	474	try expectApproxEqAbs(@as(f128, 0.070737201667702910088189851434268747), c, epsilon);
	475
	476	sincosq(37.45, &s, &c);
	477	try expectApproxEqAbs(@as(f128, -0.24654331551411356571238581321661085), s, epsilon);
	478	try expectApproxEqAbs(@as(f128, 0.96913177307077712443149563847233230), c, epsilon);
	479
	480	sincosq(89.123, &s, &c);
	481	try expectApproxEqAbs(@as(f128, 0.91616527666226951075019849560482170), s, epsilon);
	482	try expectApproxEqAbs(@as(f128, 0.40080068093548339848199454493704702), c, epsilon);
	483	}
	484
	485	test "sincos128.special" {
	486	try testSincosSpecial(f128);
	487	}

lib/compiler_rt/tan.zig+118-47

...	...	@@ -3,6 +3,7 @@
3	3	//!
4	4	//! https://git.musl-libc.org/cgit/musl/tree/src/math/tanf.c
5	5	//! https://git.musl-libc.org/cgit/musl/tree/src/math/tan.c
	6	//! https://git.musl-libc.org/cgit/musl/tree/src/math/tanl.c
6	7	//! https://golang.org/src/math/tan.go
7	8
8	9	const std = @import("std");
...	...	@@ -10,20 +11,23 @@ const builtin = @import("builtin");
10	11	const math = std.math;
11	12	const mem = std.mem;
12	13	const expect = std.testing.expect;
	14	const expectApproxEqAbs = std.testing.expectApproxEqAbs;
13	15
14	16	const kernel = @import("trig.zig");
15	17	const rem_pio2 = @import("rem_pio2.zig").rem_pio2;
16	18	const rem_pio2f = @import("rem_pio2f.zig").rem_pio2f;
	19	const rem_pio2l = @import("rem_pio2l.zig").rem_pio2l;
	20	const ld = @import("long_double.zig");
17	21
18	22	const arch = builtin.cpu.arch;
19	23	const compiler_rt = @import("../compiler_rt.zig");
20	24	const symbol = @import("../compiler_rt.zig").symbol;
21	25
22	26	comptime {
23		symbol(&__tanh, "__tanh");
	27	symbol(&tanh, "__tanh");
24	28	symbol(&tanf, "tanf");
25	29	symbol(&tan, "tan");
26		symbol(&__tanx, "__tanx");
	30	symbol(&tanx, "__tanx");
27	31	if (compiler_rt.want_ppc_abi) {
28	32	symbol(&tanq, "tanf128");
29	33	}
...	...	@@ -31,7 +35,7 @@ comptime {
31	35	symbol(&tanl, "tanl");
32	36	}
33	37
34		pub fn __tanh(x: f16) callconv(.c) f16 {
	38	pub fn tanh(x: f16) callconv(.c) f16 {
35	39	// TODO: more efficient implementation
36	40	return @floatCast(tanf(x));
37	41	}
...	...	@@ -59,20 +63,20 @@ pub fn tanf(x: f32) callconv(.c) f32 {
59	63	}
60	64	return x;
61	65	}
62		return kernel.__tandf(x, false);
	66	return kernel.tandf(x, false);
63	67	}
64	68	if (ix <= 0x407b53d1) { // \|x\| ~<= 5*pi/4
65	69	if (ix <= 0x4016cbe3) { // \|x\| ~<= 3pi/4
66		return kernel.__tandf((if (sign) x + t1pio2 else x - t1pio2), true);
	70	return kernel.tandf((if (sign) x + t1pio2 else x - t1pio2), true);
67	71	} else {
68		return kernel.__tandf((if (sign) x + t2pio2 else x - t2pio2), false);
	72	return kernel.tandf((if (sign) x + t2pio2 else x - t2pio2), false);
69	73	}
70	74	}
71	75	if (ix <= 0x40e231d5) { // \|x\| ~<= 9*pi/4
72	76	if (ix <= 0x40afeddf) { // \|x\| ~<= 7*pi/4
73		return kernel.__tandf((if (sign) x + t3pio2 else x - t3pio2), true);
	77	return kernel.tandf((if (sign) x + t3pio2 else x - t3pio2), true);
74	78	} else {
75		return kernel.__tandf((if (sign) x + t4pio2 else x - t4pio2), false);
	79	return kernel.tandf((if (sign) x + t4pio2 else x - t4pio2), false);
76	80	}
77	81	}
78	82
...	...	@@ -83,7 +87,7 @@ pub fn tanf(x: f32) callconv(.c) f32 {
83	87
84	88	var y: f64 = undefined;
85	89	const n = rem_pio2f(x, &y);
86		return kernel.__tandf(y, n & 1 != 0);
	90	return kernel.tandf(y, n & 1 != 0);
87	91	}
88	92
89	93	pub fn tan(x: f64) callconv(.c) f64 {
...	...	@@ -103,7 +107,7 @@ pub fn tan(x: f64) callconv(.c) f64 {
103	107	}
104	108	return x;
105	109	}
106		return kernel.__tan(x, 0.0, false);
	110	return kernel.tan(x, 0.0, false);
107	111	}
108	112
109	113	// tan(Inf or NaN) is NaN
...	...	@@ -113,69 +117,136 @@ pub fn tan(x: f64) callconv(.c) f64 {
113	117
114	118	var y: [2]f64 = undefined;
115	119	const n = rem_pio2(x, &y);
116		return kernel.__tan(y[0], y[1], n & 1 != 0);
	120	return kernel.tan(y[0], y[1], n & 1 != 0);
117	121	}
118	122
119		pub fn __tanx(x: f80) callconv(.c) f80 {
120		// TODO: more efficient implementation
121		return @floatCast(tanq(x));
	123	pub fn tanx(x: f80) callconv(.c) f80 {
	124	const se = ld.signExponent(x) & 0x7fff;
	125	if (se == 0x7fff) {
	126	return x - x;
	127	}
	128
	129	if (@abs(x) < kernel.pi_4) {
	130	if (se < 0x3fff - math.floatMantissaBits(f80) / 2) {
	131	if (compiler_rt.want_float_exceptions) {
	132	mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
	133	}
	134	return x;
	135	}
	136	return kernel.tanx(x, 0.0, 0);
	137	}
	138
	139	var y: [2]f80 = undefined;
	140	const n = rem_pio2l(f80, x, &y);
	141	return kernel.tanx(y[0], y[1], n & 1);
122	142	}
123	143
124	144	pub fn tanq(x: f128) callconv(.c) f128 {
125		// TODO: more correct implementation
126		return tan(@floatCast(x));
	145	const se = ld.signExponent(x) & 0x7fff;
	146	if (se == 0x7fff) {
	147	return x - x;
	148	}
	149
	150	if (@abs(x) < kernel.pi_4) {
	151	if (se < 0x3fff - math.floatMantissaBits(f128) / 2) {
	152	if (compiler_rt.want_float_exceptions) {
	153	mem.doNotOptimizeAway(if (se == 0) x * 0x1p-120 else x + 0x1p120);
	154	}
	155	return x;
	156	}
	157	return kernel.tanq(x, 0.0, 0);
	158	}
	159
	160	var y: [2]f128 = undefined;
	161	const n = rem_pio2l(f128, x, &y);
	162	return kernel.tanq(y[0], y[1], n & 1);
127	163	}
128	164
129	165	pub fn tanl(x: c_longdouble) callconv(.c) c_longdouble {
130	166	switch (@typeInfo(c_longdouble).float.bits) {
131		16 => return __tanh(x),
	167	16 => return tanh(x),
132	168	32 => return tanf(x),
133	169	64 => return tan(x),
134		80 => return __tanx(x),
	170	80 => return tanx(x),
135	171	128 => return tanq(x),
136	172	else => @compileError("unreachable"),
137	173	}
138	174	}
139	175
140		test "tan" {
141		try expect(tan(@as(f32, 0.0)) == tanf(0.0));
142		try expect(tan(@as(f64, 0.0)) == tan(0.0));
	176	fn testTanNormal(comptime T: type) !void {
	177	const f = switch (T) {
	178	f32 => tanf,
	179	f64 => tan,
	180	else => @compileError("unimplemented"),
	181	};
	182	const epsilon = 0.00001;
	183
	184	try expectApproxEqAbs(@as(T, 0.0), f(0.0), epsilon);
	185	try expectApproxEqAbs(@as(T, 0.202710), f(0.2), epsilon);
	186	try expectApproxEqAbs(@as(T, 1.240422), f(0.8923), epsilon);
	187	try expectApproxEqAbs(@as(T, 14.101420), f(1.5), epsilon);
	188	try expectApproxEqAbs(@as(T, -0.254397), f(37.45), epsilon);
	189	try expectApproxEqAbs(@as(T, 2.285837), f(89.123), epsilon);
	190	}
	191
	192	fn testTanSpecial(comptime T: type) !void {
	193	const f = switch (T) {
	194	f32 => tanf,
	195	f64 => tan,
	196	f80 => tanx,
	197	f128 => tanq,
	198	else => @compileError("unimplemented"),
	199	};
	200
	201	try expect(math.isPositiveZero(f(0.0)));
	202	try expect(math.isNegativeZero(f(-0.0)));
	203	try expect(math.isNan(f(math.inf(f32))));
	204	try expect(math.isNan(f(-math.inf(f32))));
	205	try expect(math.isNan(f(math.nan(f32))));
143	206	}
144	207
145		test "tan32" {
146		const epsilon = 0.00001;
	208	test "tan32.normal" {
	209	try testTanNormal(f32);
	210	}
	211
	212	test "tan64.normal" {
	213	try testTanNormal(f64);
	214	}
	215
	216	test "tan80.normal" {
	217	const epsilon = math.floatEps(f80);
147	218
148		try expect(math.approxEqAbs(f32, tanf(0.0), 0.0, epsilon));
149		try expect(math.approxEqAbs(f32, tanf(0.2), 0.202710, epsilon));
150		try expect(math.approxEqAbs(f32, tanf(0.8923), 1.240422, epsilon));
151		try expect(math.approxEqAbs(f32, tanf(1.5), 14.101420, epsilon));
152		try expect(math.approxEqAbs(f32, tanf(37.45), -0.254397, epsilon));
153		try expect(math.approxEqAbs(f32, tanf(89.123), 2.285852, epsilon));
	219	try expectApproxEqAbs(@as(f80, 0.0), tanx(0.0), epsilon);
	220	try expectApproxEqAbs(@as(f80, 0.2027100355086724833213582716475345), tanx(0.2), epsilon);
	221	try expectApproxEqAbs(@as(f80, 1.2404217445497097995561220131857544), tanx(0.8923), epsilon);
	222	try expectApproxEqAbs(@as(f80, 14.10141994717171938764), tanx(1.5), epsilon);
	223	try expectApproxEqAbs(@as(f80, -0.25439607116885656232), tanx(37.45), epsilon);
	224	try expectApproxEqAbs(@as(f80, 2.2858376251355320963), tanx(89.123), epsilon);
154	225	}
155	226
156		test "tan64" {
157		const epsilon = 0.000001;
	227	test "tan128.normal" {
	228	const epsilon = math.floatEps(f128);
158	229
159		try expect(math.approxEqAbs(f64, tan(0.0), 0.0, epsilon));
160		try expect(math.approxEqAbs(f64, tan(0.2), 0.202710, epsilon));
161		try expect(math.approxEqAbs(f64, tan(0.8923), 1.240422, epsilon));
162		try expect(math.approxEqAbs(f64, tan(1.5), 14.101420, epsilon));
163		try expect(math.approxEqAbs(f64, tan(37.45), -0.254397, epsilon));
164		try expect(math.approxEqAbs(f64, tan(89.123), 2.2858376, epsilon));
	230	try expectApproxEqAbs(@as(f128, 0.0), tanq(0.0), epsilon);
	231	try expectApproxEqAbs(@as(f128, 0.2027100355086724833213582716475345), tanq(0.2), epsilon);
	232	try expectApproxEqAbs(@as(f128, 1.2404217445497097995561220131857544), tanq(0.8923), epsilon);
	233	try expectApproxEqAbs(@as(f128, 14.101419947171719387646083651987755), tanq(1.5), epsilon);
	234	try expectApproxEqAbs(@as(f128, -0.2543960711688565630469573224504774), tanq(37.45), epsilon);
	235	try expectApproxEqAbs(@as(f128, 2.2858376251355321074066028114094292), tanq(89.123), epsilon);
165	236	}
166	237
167	238	test "tan32.special" {
168		try expect(tanf(0.0) == 0.0);
169		try expect(tanf(-0.0) == -0.0);
170		try expect(math.isNan(tanf(math.inf(f32))));
171		try expect(math.isNan(tanf(-math.inf(f32))));
172		try expect(math.isNan(tanf(math.nan(f32))));
	239	try testTanSpecial(f32);
173	240	}
174	241
175	242	test "tan64.special" {
176		try expect(tan(0.0) == 0.0);
177		try expect(tan(-0.0) == -0.0);
178		try expect(math.isNan(tan(math.inf(f64))));
179		try expect(math.isNan(tan(-math.inf(f64))));
180		try expect(math.isNan(tan(math.nan(f64))));
	243	try testTanSpecial(f64);
	244	}
	245
	246	test "tan80.special" {
	247	try testTanSpecial(f80);
	248	}
	249
	250	test "tan128.special" {
	251	try testTanSpecial(f128);
181	252	}

lib/compiler_rt/trig.zig+256-6

...	...	@@ -7,6 +7,13 @@
7	7	// https://git.musl-libc.org/cgit/musl/tree/src/math/__sindf.c
8	8	// https://git.musl-libc.org/cgit/musl/tree/src/math/__tand.c
9	9	// https://git.musl-libc.org/cgit/musl/tree/src/math/__tandf.c
	10	// https://git.musl-libc.org/cgit/musl/tree/src/math/__sinl.c
	11	// https://git.musl-libc.org/cgit/musl/tree/src/math/__cosl.c
	12	// https://git.musl-libc.org/cgit/musl/tree/src/math/__tanl.c
	13
	14	const std = @import("std");
	15
	16	pub const pi_4 = std.math.pi / 4.0;
10	17
11	18	/// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
12	19	/// Input x is assumed to be bounded by ~pi/4 in magnitude.
...	...	@@ -43,7 +50,7 @@
43	50	/// expression for cos(). Retention happens in all cases tested
44	51	/// under FreeBSD, so don't pessimize things by forcibly clipping
45	52	/// any extra precision in w.
46		pub fn __cos(x: f64, y: f64) f64 {
	53	pub fn cos(x: f64, y: f64) f64 {
47	54	const C1 = 4.16666666666666019037e-02; // 0x3FA55555, 0x5555554C
48	55	const C2 = -1.38888888888741095749e-03; // 0xBF56C16C, 0x16C15177
49	56	const C3 = 2.48015872894767294178e-05; // 0x3EFA01A0, 0x19CB1590
...	...	@@ -59,7 +66,7 @@ pub fn __cos(x: f64, y: f64) f64 {
59	66	return w + (((1.0 - w) - hz) + (z * r - x * y));
60	67	}
61	68
62		pub fn __cosdf(x: f64) f32 {
	69	pub fn cosdf(x: f64) f32 {
63	70	// \|cos(x) - c(x)\| < 2**-34.1 (~[-5.37e-11, 5.295e-11]).
64	71	const C0 = -0x1ffffffd0c5e81.0p-54; // -0.499999997251031003120
65	72	const C1 = 0x155553e1053a42.0p-57; // 0.0416666233237390631894
...	...	@@ -73,6 +80,46 @@ pub fn __cosdf(x: f64) f32 {
73	80	return @floatCast(((1.0 + z * C0) + w * C1) + (w * z) * r);
74	81	}
75	82
	83	pub fn cosx(x: f80, y: f80) f80 {
	84	const C1: f80 = 0.0416666666666666666136;
	85	const C2: f64 = -0.0013888888888888874;
	86	const C3: f64 = 0.000024801587301571716;
	87	const C4: f64 = -0.00000027557319215507120;
	88	const C5: f64 = 0.0000000020876754400407278;
	89	const C6: f64 = -1.1470297442401303e-11;
	90	const C7: f64 = 4.7383039476436467e-14;
	91
	92	const z = x * x;
	93	const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 +
	94	z * (C5 + z * (C6 + z * C7))))));
	95	const hz = 0.5 * z;
	96	const w = 1.0 - hz;
	97
	98	return w + (((1.0 - w) - hz) + (z * r - x * y));
	99	}
	100
	101	pub fn cosq(x: f128, y: f128) f128 {
	102	const C1: f128 = 0.04166666666666666666666666666666658424671;
	103	const C2: f128 = -0.001388888888888888888888888888863490893732;
	104	const C3: f128 = 0.00002480158730158730158730158600795304914210;
	105	const C4: f128 = -0.2755731922398589065255474947078934284324e-6;
	106	const C5: f128 = 0.2087675698786809897659225313136400793948e-8;
	107	const C6: f128 = -0.1147074559772972315817149986812031204775e-10;
	108	const C7: f128 = 0.4779477332386808976875457937252120293400e-13;
	109	const C8: f64 = -0.1561920696721507929516718307820958119868e-15;
	110	const C9: f64 = 0.4110317413744594971475941557607804508039e-18;
	111	const C10: f64 = -0.8896592467191938803288521958313920156409e-21;
	112	const C11: f64 = 0.1601061435794535138244346256065192782581e-23;
	113
	114	const z = x * x;
	115	const r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 +
	116	z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
	117	const hz = 0.5 * z;
	118	const w = 1.0 - hz;
	119
	120	return w + (((1.0 - w) - hz) + (z * r - x * y));
	121	}
	122
76	123	/// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
77	124	/// Input x is assumed to be bounded by ~pi/4 in magnitude.
78	125	/// Input y is the tail of x.
...	...	@@ -100,7 +147,7 @@ pub fn __cosdf(x: f64) f32 {
100	147	/// r = x (S2+x (S3+x (S4+x (S5+x *S6))))
101	148	/// then 3 2
102	149	/// sin(x) = x + (S1x + (x (r-y/2)+y))
103		pub fn __sin(x: f64, y: f64, iy: i32) f64 {
	150	pub fn sin(x: f64, y: f64, iy: i32) f64 {
104	151	const S1 = -1.66666666666666324348e-01; // 0xBFC55555, 0x55555549
105	152	const S2 = 8.33333333332248946124e-03; // 0x3F811111, 0x1110F8A6
106	153	const S3 = -1.98412698298579493134e-04; // 0xBF2A01A0, 0x19C161D5
...	...	@@ -119,7 +166,7 @@ pub fn __sin(x: f64, y: f64, iy: i32) f64 {
119	166	}
120	167	}
121	168
122		pub fn __sindf(x: f64) f32 {
	169	pub fn sindf(x: f64) f32 {
123	170	// \|sin(x)/x - s(x)\| < 2**-37.5 (~[-4.89e-12, 4.824e-12]).
124	171	const S1 = -0x15555554cbac77.0p-55; // -0.166666666416265235595
125	172	const S2 = 0x111110896efbb2.0p-59; // 0.0083333293858894631756
...	...	@@ -134,6 +181,52 @@ pub fn __sindf(x: f64) f32 {
134	181	return @floatCast((x + s * (S1 + z * S2)) + s * w * r);
135	182	}
136	183
	184	pub fn sinx(x: f80, y: f80, iy: i32) f80 {
	185	const S1: f80 = -0.166666666666666666671;
	186	const S2: f64 = 0.0083333333333333332;
	187	const S3: f64 = -0.00019841269841269427;
	188	const S4: f64 = 0.0000027557319223597490;
	189	const S5: f64 = -0.000000025052108218074604;
	190	const S6: f64 = 1.6059006598854211e-10;
	191	const S7: f64 = -7.6429779983024564e-13;
	192	const S8: f64 = 2.6174587166648325e-15;
	193
	194	const z = x * x;
	195	const v = z * x;
	196	const r = S2 + z * (S3 + z * (S4 + z * (S5 +
	197	z * (S6 + z * (S7 + z * S8)))));
	198
	199	if (iy == 0)
	200	return x + v * (S1 + z * r);
	201
	202	return x - ((z * (0.5 * y - v * r) - y) - v * S1);
	203	}
	204
	205	pub fn sinq(x: f128, y: f128, iy: i32) f128 {
	206	const S1: f128 = -0.16666666666666666666666666666666666606732416116558;
	207	const S2: f128 = 0.0083333333333333333333333333333331135404851288270047;
	208	const S3: f128 = -0.00019841269841269841269841269839935785325638310428717;
	209	const S4: f128 = 0.27557319223985890652557316053039946268333231205686e-5;
	210	const S5: f128 = -0.25052108385441718775048214826384312253862930064745e-7;
	211	const S6: f128 = 0.16059043836821614596571832194524392581082444805729e-9;
	212	const S7: f128 = -0.76471637318198151807063387954939213287488216303768e-12;
	213	const S8: f128 = 0.28114572543451292625024967174638477283187397621303e-14;
	214	const S9: f64 = -0.82206352458348947812512122163446202498005154296863e-17;
	215	const S10: f64 = 0.19572940011906109418080609928334380560135358385256e-19;
	216	const S11: f64 = -0.38680813379701966970673724299207480965452616911420e-22;
	217	const S12: f64 = 0.64038150078671872796678569586315881020659912139412e-25;
	218
	219	const z = x * x;
	220	const v = z * x;
	221	const r = S2 + z * (S3 + z * (S4 + z * (S5 + z * (S6 + z * (S7 + z * (S8 +
	222	z * (S9 + z * (S10 + z * (S11 + z * S12)))))))));
	223
	224	if (iy == 0)
	225	return x + v * (S1 + z * r);
	226
	227	return x - ((z * (0.5 * y - v * r) - y) - v * S1);
	228	}
	229
137	230	/// kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
138	231	/// Input x is assumed to be bounded by ~pi/4 in magnitude.
139	232	/// Input y is the tail of x.
...	...	@@ -166,7 +259,7 @@ pub fn __sindf(x: f64) f32 {
166	259	/// 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
167	260	/// tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
168	261	/// = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
169		pub fn __tan(x_: f64, y_: f64, odd: bool) f64 {
	262	pub fn tan(x_: f64, y_: f64, odd: bool) f64 {
170	263	var x = x_;
171	264	var y = y_;
172	265
...	...	@@ -239,7 +332,7 @@ pub fn __tan(x_: f64, y_: f64, odd: bool) f64 {
239	332	return a0 + a * (1.0 + a0 * w0 + a0 * v);
240	333	}
241	334
242		pub fn __tandf(x: f64, odd: bool) f32 {
	335	pub fn tandf(x: f64, odd: bool) f32 {
243	336	// \|tan(x)/x - t(x)\| < 2**-25.5 (~[-2e-08, 2e-08]).
244	337	const T = [_]f64{
245	338	0x15554d3418c99f.0p-54, // 0.333331395030791399758
...	...	@@ -271,3 +364,160 @@ pub fn __tandf(x: f64, odd: bool) f32 {
271	364	const r0 = (x + s * u) + (s * w) * (t + w * r);
272	365	return @floatCast(if (odd) -1.0 / r0 else r0);
273	366	}
	367
	368	pub fn tanx(x_: f80, y_: f80, odd: i32) f80 {
	369	const pio4: f80 = 0.785398163397448309628;
	370	const pio4lo: f80 = -1.25413940316708300586e-20;
	371
	372	const T3: f80 = 0.333333333333333333180;
	373	const T5: f80 = 0.133333333333333372290;
	374	const T7: f80 = 0.0539682539682504975744;
	375	const T9: f64 = 0.021869488536312216;
	376	const T11: f64 = 0.0088632355256619590;
	377	const T13: f64 = 0.0035921281113786528;
	378	const T15: f64 = 0.0014558334756312418;
	379	const T17: f64 = 0.00059003538700862256;
	380	const T19: f64 = 0.00023907843576635544;
	381	const T21: f64 = 0.000097154625656538905;
	382	const T23: f64 = 0.000038440165747303162;
	383	const T25: f64 = 0.000018082171885432524;
	384	const T27: f64 = 0.0000024196006108814377;
	385	const T29: f64 = 0.0000078293456938132840;
	386	const T31: f64 = -0.0000032609076735050182;
	387	const T33: f64 = 0.0000023261313142559411;
	388
	389	var x = x_;
	390	var y = y_;
	391	const big = @abs(x) >= 0.67434;
	392	var sign: i8 = 0;
	393
	394	if (big) {
	395	if (x < 0) {
	396	sign = -1;
	397	x = -x;
	398	y = -y;
	399	}
	400	x = (pio4 - x) + (pio4lo - y);
	401	y = 0.0;
	402	}
	403
	404	var z = x * x;
	405	var w = z * z;
	406
	407	var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
	408	w * (T25 + w * (T29 + w * T33))))));
	409
	410	var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
	411	w * (T27 + w * T31))))));
	412
	413	var s = z * x;
	414	r = y + z * (s * (r + v) + y) + T3 * s;
	415	w = x + r;
	416
	417	if (big) {
	418	s = @as(f80, @floatFromInt(1 - 2 * odd));
	419	v = s - 2.0 * (x + (r - w * w / (w + s)));
	420	return if (sign == -1) -v else v;
	421	}
	422
	423	if (odd == 0) {
	424	return w;
	425	}
	426
	427	// if allow error up to 2 ulp, simply return
	428	// -1.0 / (x+r) here
	429	//
	430	// compute -1.0 / (x+r) accurately
	431	z = w + 0x1p32 - 0x1p32;
	432	v = r - (z - x);
	433	const a = -1.0 / w;
	434	const t = a + 0x1p32 - 0x1p32;
	435	s = 1.0 + t * z;
	436	return t + a * (s + t * v);
	437	}
	438
	439	pub fn tanq(x_: f128, y_: f128, odd: i32) f128 {
	440	const pio4: f128 = 0x1.921fb54442d18469898cc51701b8p-1;
	441	const pio4lo: f128 = 0x1.cd129024e088a67cc74020bbea60p-116;
	442
	443	const T3: f128 = 0x1.5555555555555555555555555553p-2;
	444	const T5: f128 = 0x1.1111111111111111111111111eb5p-3;
	445	const T7: f128 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5;
	446	const T9: f128 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6;
	447	const T11: f128 = 0x1.226e355e6c23c8f5b4f5762322eep-7;
	448	const T13: f128 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9;
	449	const T15: f128 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10;
	450	const T17: f128 = 0x1.355824803674477dfcf726649efep-11;
	451	const T19: f128 = 0x1.f57d7734d1656e0aceb716f614c2p-13;
	452	const T21: f128 = 0x1.967e18afcb180ed942dfdc518d6cp-14;
	453	const T23: f128 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15;
	454	const T25: f128 = 0x1.0b132d39f055c81be49eff7afd50p-16;
	455	const T27: f128 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18;
	456	const T29: f128 = 0x1.5ef2daf21d1113df38d0fbc00267p-19;
	457	const T31: f128 = 0x1.1c77d6eac0234988cdaa04c96626p-20;
	458	const T33: f128 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22;
	459	const T35: f128 = 0x1.75c7357d0298c01a31d0a6f7d518p-23;
	460	const T37: f128 = 0x1.2f3190f4718a9a520f98f50081fcp-24;
	461	const T39: f64 = 0.000000028443389121318352;
	462	const T41: f64 = 0.000000011981013102001973;
	463	const T43: f64 = 0.0000000038303578044958070;
	464	const T45: f64 = 0.0000000034664378216909893;
	465	const T47: f64 = -0.0000000015090641701997785;
	466	const T49: f64 = 0.0000000029449552300483952;
	467	const T51: f64 = -0.0000000022006995706097711;
	468	const T53: f64 = 0.0000000015468200913196612;
	469	const T55: f64 = -0.00000000061311613386849674;
	470	const T57: f64 = 1.4912469681508012e-10;
	471
	472	var x = x_;
	473	var y = y_;
	474
	475	const big = @abs(x) >= 0.67434;
	476	var sign: i8 = 0;
	477
	478	if (big) {
	479	if (x < 0) {
	480	sign = -1;
	481	x = -x;
	482	y = -y;
	483	}
	484	x = (pio4 - x) + (pio4lo - y);
	485	y = 0.0;
	486	}
	487
	488	var z = x * x;
	489	var w = z * z;
	490
	491	var r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
	492	w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 +
	493	w * (T45 + w * (T49 + w * (T53 + w * T57))))))))))));
	494
	495	var v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
	496	w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 +
	497	w * (T47 + w * (T51 + w * T55))))))))))));
	498
	499	var s = z * x;
	500	r = y + z * (s * (r + v) + y) + T3 * s;
	501	w = x + r;
	502
	503	if (big) {
	504	s = @as(f128, @floatFromInt(1 - 2 * odd));
	505	v = s - 2.0 * (x + (r - w * w / (w + s)));
	506	return if (sign == -1) -v else v;
	507	}
	508
	509	if (odd == 0) {
	510	return w;
	511	}
	512
	513	// if allow error up to 2 ulp, simply return
	514	// -1.0 / (x+r) here
	515	//
	516	// compute -1.0 / (x+r) accurately
	517	z = w + 0x1p32 - 0x1p32;
	518	v = r - (z - x);
	519	const a = -1.0 / w;
	520	const t = a + 0x1p32 - 0x1p32;
	521	s = 1.0 + t * z;
	522	return t + a * (s + t * v);
	523	}

lib/libc/mingw/math/arm-common/sincosl.c deleted-13

...	...	@@ -1,13 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6
7		#include <math.h>
8
9		void sincosl(long double x, long double s, long double c)
10		{
11		*s = sinl(x);
12		*c = cosl(x);
13		}

lib/libc/mingw/math/arm/sincos.S deleted-30

...	...	@@ -1,30 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <_mingw_mac.h>
7
8		.file "sincos.S"
9		.text
10		.align 2
11		/* zig patch: remove sincos symbol because sincos in compiler_rt is used instead */
12		.globl __MINGW_USYMBOL(sincosl)
13		.def __MINGW_USYMBOL(sincosl); .scl 2; .type 32; .endef
14		__MINGW_USYMBOL(sincosl):
15		push {r4, r5, r11, lr}
16		add r11, sp, #8
17		vpush {d8}
18
19		mov r4, r0
20		mov r5, r1
21		vmov.f64 d8, d0
22		bl sin
23		vstr d0, [r4]
24
25		vmov.f64 d0, d8
26		bl cos
27		vstr d0, [r5]
28
29		vpop {d8}
30		pop {r4, r5, r11, pc}

lib/libc/mingw/math/arm64/rint.c deleted-12

...	...	@@ -1,12 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <math.h>
7
8		double rint (double x) {
9		double retval = 0.0;
10		__asm__ __volatile__ ("frintx %d0, %d1\n\t" : "=w" (retval) : "w" (x));
11		return retval;
12		}

lib/libc/mingw/math/arm64/rintf.c deleted-12

...	...	@@ -1,12 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <math.h>
7
8		float rintf (float x) {
9		float retval = 0.0F;
10		__asm__ __volatile__ ("frintx %s0, %s1\n\t" : "=w" (retval) : "w" (x));
11		return retval;
12		}

lib/libc/mingw/math/arm64/sincos.S deleted-32

...	...	@@ -1,32 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <_mingw_mac.h>
7
8		.file "sincos.S"
9		.text
10		.align 2
11		/* zig patch: remove sincos symbol because sincos in compiler_rt is used instead */
12		.globl __MINGW_USYMBOL(sincosl)
13		.def __MINGW_USYMBOL(sincosl); .scl 2; .type 32; .endef
14		__MINGW_USYMBOL(sincosl):
15		str d8, [sp, #-32]!
16		str x30, [sp, #8]
17		stp x19, x20, [sp, #16]
18
19		mov x19, x0
20		mov x20, x1
21		fmov d8, d0
22		bl sin
23		str d0, [x19]
24
25		fmov d0, d8
26		bl cos
27		str d0, [x20]
28
29		ldp x19, x20, [sp, #16]
30		ldr x30, [sp, #8]
31		ldr d8, [sp], #32
32		ret

lib/libc/mingw/math/frexpf.c deleted-13

...	...	@@ -1,13 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		extern double __cdecl frexp(double _X,int *_Y);
7
8		float frexpf (float, int *);
9		float frexpf (float x, int *expn)
10		{
11		return (float)frexp(x, expn);
12		}
13

lib/libc/mingw/math/frexpl.c deleted-71

...	...	@@ -1,71 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		long double frexpl(long double value, int* exp);
7
8		#if __SIZEOF_LONG_DOUBLE__ == __SIZEOF_DOUBLE__
9
10		double frexp(double value, int* exp);
11
12		/* On ARM `long double` is 64 bits. */
13		long double frexpl(long double value, int* exp)
14		{
15		return frexp(value, exp);
16		}
17
18		#elif defined(_AMD64_) \|\| defined(__x86_64__) \|\| defined(_X86_) \|\| defined(__i386__)
19
20		#include <stdint.h>
21
22		/* https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format */
23		typedef union x87reg_ {
24		struct __attribute__((__packed__)) {
25		uint64_t f64;
26		uint16_t exp : 15;
27		uint16_t sgn : 1;
28		};
29		long double f;
30		} x87reg;
31
32		long double frexpl(long double value, int* exp)
33		{
34		int n;
35		x87reg reg;
36		reg.f = value;
37		if(reg.exp == 0x7FFF) {
38		/* The value is an infinity or NaN.
39		* Store zero in `exp`. Return the value as is. /
40		*exp = 0;
41		return reg.f;
42		}
43		if(reg.exp != 0) {
44		/* The value is normalized.
45		* Extract and zero out the exponent. */
46		*exp = reg.exp - 0x3FFE;
47		reg.exp = 0x3FFE;
48		return reg.f;
49		}
50		if(reg.f64 == 0) {
51		/* The value is zero.
52		* Store zero in `*exp`. Return the value as is.
53		* Note the signness. */
54		*exp = 0;
55		return reg.f;
56		}
57		/* The value is denormalized.
58		* Extract the exponent, normalize the value, then zero out
59		* the exponent. Note that x87 uses an explicit leading bit. */
60		n = __builtin_clzll(reg.f64);
61		reg.f64 <<= n;
62		*exp = 1 - 0x3FFE - n;
63		reg.exp = 0x3FFE;
64		return reg.f;
65		}
66
67		#else
68
69		#error Please add `frexpl()` implementation for this platform.
70
71		#endif

lib/libc/mingw/math/x86/cos.def.h deleted-65

...	...	@@ -1,65 +0,0 @@
1		/*
2		This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4		Copyright (c) 2009, 2010 by the mingw-w64 project
5
6		See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8		This license has been certified as open source. It has also been designated
9		as GPL compatible by the Free Software Foundation (FSF).
10
11		Redistribution and use in source and binary forms, with or without
12		modification, are permitted provided that the following conditions are met:
13
14		1. Redistributions in source code must retain the accompanying copyright
15		notice, this list of conditions, and the following disclaimer.
16		2. Redistributions in binary form must reproduce the accompanying
17		copyright notice, this list of conditions, and the following disclaimer
18		in the documentation and/or other materials provided with the
19		distribution.
20		3. Names of the copyright holders must not be used to endorse or promote
21		products derived from this software without prior written permission
22		from the copyright holders.
23		4. The right to distribute this software or to use it for any purpose does
24		not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25		the copyright holders. Use of them is covered by separate agreement
26		with the copyright holders.
27		5. If any files are modified, you must cause the modified files to carry
28		prominent notices stating that you changed the files and the date of
29		any change.
30
31		Disclaimer
32
33		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34		OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35		OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36		EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37		INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38		LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39		OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40		LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41		NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42		EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43		*/
44
45		#include "../complex/complex_internal.h"
46		#include <errno.h>
47
48		extern long double __cosl_internal (long double);
49
50		__FLT_TYPE
51		__FLT_ABI(cos) (__FLT_TYPE x)
52		{
53		int x_class = fpclassify (x);
54		if (x_class == FP_NAN)
55		{
56		__FLT_RPT_DOMAIN ("cos", x, 0.0, x);
57		return x;
58		}
59		else if (x_class == FP_INFINITE)
60		{
61		__FLT_RPT_DOMAIN ("cos", x, 0.0, __FLT_NAN);
62		return __FLT_NAN;
63		}
64		return (__FLT_TYPE) __cosl_internal ((long double) x);
65		}

lib/libc/mingw/math/x86/cosl.c deleted-46

...	...	@@ -1,46 +0,0 @@
1		/*
2		This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4		Copyright (c) 2009, 2010 by the mingw-w64 project
5
6		See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8		This license has been certified as open source. It has also been designated
9		as GPL compatible by the Free Software Foundation (FSF).
10
11		Redistribution and use in source and binary forms, with or without
12		modification, are permitted provided that the following conditions are met:
13
14		1. Redistributions in source code must retain the accompanying copyright
15		notice, this list of conditions, and the following disclaimer.
16		2. Redistributions in binary form must reproduce the accompanying
17		copyright notice, this list of conditions, and the following disclaimer
18		in the documentation and/or other materials provided with the
19		distribution.
20		3. Names of the copyright holders must not be used to endorse or promote
21		products derived from this software without prior written permission
22		from the copyright holders.
23		4. The right to distribute this software or to use it for any purpose does
24		not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25		the copyright holders. Use of them is covered by separate agreement
26		with the copyright holders.
27		5. If any files are modified, you must cause the modified files to carry
28		prominent notices stating that you changed the files and the date of
29		any change.
30
31		Disclaimer
32
33		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34		OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35		OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36		EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37		INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38		LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39		OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40		LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41		NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42		EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43		*/
44
45		#define _NEW_COMPLEX_LDOUBLE 1
46		#include "cos.def.h"

lib/libc/mingw/math/x86/cosl_internal.S deleted-55

...	...	@@ -1,55 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <_mingw_mac.h>
7
8		.file "cosl_internal.S"
9		.text
10		#ifdef __x86_64__
11		.align 8
12		#else
13		.align 4
14		#endif
15		.globl __MINGW_USYMBOL(__cosl_internal)
16		.def __MINGW_USYMBOL(__cosl_internal); .scl 2; .type 32; .endef
17		__MINGW_USYMBOL(__cosl_internal):
18		#ifdef __x86_64__
19		fldt (%rdx)
20		fcos
21		fnstsw %ax
22		testl $0x400,%eax
23		jz 1f
24		fldpi
25		fadd %st(0)
26		fxch %st(1)
27		2: fprem1
28		fnstsw %ax
29		testl $0x400,%eax
30		jnz 2b
31		fstp %st(1)
32		fcos
33		1: movq %rcx,%rax
34		movq $0,8(%rcx)
35		fstpt (%rcx)
36		ret
37		#else
38		fldt 4(%esp)
39		fcos
40		fnstsw %ax
41		testl $0x400,%eax
42		jnz 1f
43		ret
44		1: fldpi
45		fadd %st(0)
46		fxch %st(1)
47		2: fprem1
48		fnstsw %ax
49		testl $0x400,%eax
50		jnz 2b
51		fstp %st(1)
52		fcos
53		ret
54		#endif
55

lib/libc/mingw/math/x86/sin.def.h deleted-65

...	...	@@ -1,65 +0,0 @@
1		/*
2		This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4		Copyright (c) 2009, 2010 by the mingw-w64 project
5
6		See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8		This license has been certified as open source. It has also been designated
9		as GPL compatible by the Free Software Foundation (FSF).
10
11		Redistribution and use in source and binary forms, with or without
12		modification, are permitted provided that the following conditions are met:
13
14		1. Redistributions in source code must retain the accompanying copyright
15		notice, this list of conditions, and the following disclaimer.
16		2. Redistributions in binary form must reproduce the accompanying
17		copyright notice, this list of conditions, and the following disclaimer
18		in the documentation and/or other materials provided with the
19		distribution.
20		3. Names of the copyright holders must not be used to endorse or promote
21		products derived from this software without prior written permission
22		from the copyright holders.
23		4. The right to distribute this software or to use it for any purpose does
24		not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25		the copyright holders. Use of them is covered by separate agreement
26		with the copyright holders.
27		5. If any files are modified, you must cause the modified files to carry
28		prominent notices stating that you changed the files and the date of
29		any change.
30
31		Disclaimer
32
33		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34		OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35		OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36		EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37		INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38		LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39		OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40		LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41		NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42		EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43		*/
44
45		#include "../complex/complex_internal.h"
46		#include <errno.h>
47
48		extern long double __sinl_internal (long double);
49
50		__FLT_TYPE
51		__FLT_ABI(sin) (__FLT_TYPE x)
52		{
53		int x_class = fpclassify (x);
54		if (x_class == FP_NAN)
55		{
56		__FLT_RPT_DOMAIN ("sin", x, 0.0, x);
57		return x;
58		}
59		else if (x_class == FP_INFINITE)
60		{
61		__FLT_RPT_DOMAIN ("sin", x, 0.0, __FLT_NAN);
62		return __FLT_NAN;
63		}
64		return (__FLT_TYPE) __sinl_internal ((long double) x);
65		}

lib/libc/mingw/math/x86/sinl.c deleted-46

...	...	@@ -1,46 +0,0 @@
1		/*
2		This Software is provided under the Zope Public License (ZPL) Version 2.1.
3
4		Copyright (c) 2009, 2010 by the mingw-w64 project
5
6		See the AUTHORS file for the list of contributors to the mingw-w64 project.
7
8		This license has been certified as open source. It has also been designated
9		as GPL compatible by the Free Software Foundation (FSF).
10
11		Redistribution and use in source and binary forms, with or without
12		modification, are permitted provided that the following conditions are met:
13
14		1. Redistributions in source code must retain the accompanying copyright
15		notice, this list of conditions, and the following disclaimer.
16		2. Redistributions in binary form must reproduce the accompanying
17		copyright notice, this list of conditions, and the following disclaimer
18		in the documentation and/or other materials provided with the
19		distribution.
20		3. Names of the copyright holders must not be used to endorse or promote
21		products derived from this software without prior written permission
22		from the copyright holders.
23		4. The right to distribute this software or to use it for any purpose does
24		not give you the right to use Servicemarks (sm) or Trademarks (tm) of
25		the copyright holders. Use of them is covered by separate agreement
26		with the copyright holders.
27		5. If any files are modified, you must cause the modified files to carry
28		prominent notices stating that you changed the files and the date of
29		any change.
30
31		Disclaimer
32
33		THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
34		OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
35		OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
36		EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
37		INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
38		LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39		OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40		LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41		NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
42		EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43		*/
44
45		#define _NEW_COMPLEX_LDOUBLE 1
46		#include "sin.def.h"

lib/libc/mingw/math/x86/sinl_internal.S deleted-58

...	...	@@ -1,58 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <_mingw_mac.h>
7
8		.file "sinl_internal.S"
9		.text
10		#ifdef __x86_64__
11		.align 8
12		#else
13		.align 4
14		#endif
15		.globl __MINGW_USYMBOL(__sinl_internal)
16		.def __MINGW_USYMBOL(__sinl_internal); .scl 2; .type 32; .endef
17		__MINGW_USYMBOL(__sinl_internal):
18		#ifdef __x86_64__
19		fldt (%rdx)
20		fsin
21		fnstsw %ax
22		testl $0x400,%eax
23		jnz 1f
24		movq %rcx,%rax
25		movq $0,8(%rcx)
26		fstpt (%rcx)
27		ret
28		1: fldpi
29		fadd %st(0)
30		fxch %st(1)
31		2: fprem1
32		fnstsw %ax
33		testl $0x400,%eax
34		jnz 2b
35		fstp %st(1)
36		fsin
37		movq %rcx,%rax
38		movq $0,8(%rcx)
39		fstpt (%rcx)
40		ret
41		#else
42		fldt 4(%esp)
43		fsin
44		fnstsw %ax
45		testl $0x400,%eax
46		jnz 1f
47		ret
48		1: fldpi
49		fadd %st(0)
50		fxch %st(1)
51		2: fprem1
52		fnstsw %ax
53		testl $0x400,%eax
54		jnz 2b
55		fstp %st(1)
56		fsin
57		ret
58		#endif

lib/libc/mingw/math/x86/tanl.S deleted-62

...	...	@@ -1,62 +0,0 @@
1		/**
2		* This file has no copyright assigned and is placed in the Public Domain.
3		* This file is part of the mingw-w64 runtime package.
4		* No warranty is given; refer to the file DISCLAIMER.PD within this package.
5		*/
6		#include <_mingw_mac.h>
7
8		.file "tanl.S"
9		.text
10		#ifdef __x86_64__
11		.align 8
12		#else
13		.align 4
14		#endif
15		.globl __MINGW_USYMBOL(tanl)
16		.def __MINGW_USYMBOL(tanl); .scl 2; .type 32; .endef
17		__MINGW_USYMBOL(tanl):
18		#ifdef __x86_64__
19		fldt (%rdx)
20		fptan
21		fnstsw %ax
22		testl $0x400,%eax
23		jnz 1f
24		fstp %st(0)
25		movq %rcx,%rax
26		movq $0,8(%rcx)
27		fstpt (%rcx)
28		ret
29		1: fldpi
30		fadd %st(0)
31		fxch %st(1)
32		2: fprem1
33		fstsw %ax
34		testl $0x400,%eax
35		jnz 2b
36		fstp %st(1)
37		fptan
38		fstp %st(0)
39		movq %rcx,%rax
40		movq $0,8(%rcx)
41		fstpt (%rcx)
42		ret
43		#else
44		fldt 4(%esp)
45		fptan
46		fnstsw %ax
47		testl $0x400,%eax
48		jnz 1f
49		fstp %st(0)
50		ret
51		1: fldpi
52		fadd %st(0)
53		fxch %st(1)
54		2: fprem1
55		fstsw %ax
56		testl $0x400,%eax
57		jnz 2b
58		fstp %st(1)
59		fptan
60		fstp %st(0)
61		ret
62		#endif

lib/libc/musl/src/math/__cosl.c deleted-96

...	...	@@ -1,96 +0,0 @@
1		/* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */
2		/* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */
3		/*
4		* ====================================================
5		* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6		* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
7		*
8		* Developed at SunSoft, a Sun Microsystems, Inc. business.
9		* Permission to use, copy, modify, and distribute this
10		* software is freely granted, provided that this notice
11		* is preserved.
12		* ====================================================
13		*/
14
15
16		#include "libm.h"
17
18		#if (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
19		#if LDBL_MANT_DIG == 64
20		/*
21		* ld80 version of __cos.c. See __cos.c for most comments.
22		*/
23		/*
24		* Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]:
25		* \|cos(x) - c(x)\| < 2**-75.1
26		*
27		* The coefficients of c(x) were generated by a pari-gp script using
28		* a Remez algorithm that searches for the best higher coefficients
29		* after rounding leading coefficients to a specified precision.
30		*
31		* Simpler methods like Chebyshev or basic Remez barely suffice for
32		* cos() in 64-bit precision, because we want the coefficient of x^2
33		* to be precisely -0.5 so that multiplying by it is exact, and plain
34		* rounding of the coefficients of a good polynomial approximation only
35		* gives this up to about 64-bit precision. Plain rounding also gives
36		* a mediocre approximation for the coefficient of x^4, but a rounding
37		* error of 0.5 ulps for this coefficient would only contribute ~0.01
38		* ulps to the final error, so this is unimportant. Rounding errors in
39		* higher coefficients are even less important.
40		*
41		* In fact, coefficients above the x^4 one only need to have 53-bit
42		* precision, and this is more efficient. We get this optimization
43		* almost for free from the complications needed to search for the best
44		* higher coefficients.
45		*/
46		static const long double
47		C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */
48		static const double
49		C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */
50		C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */
51		C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */
52		C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */
53		C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */
54		C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */
55		#define POLY(z) (z(C1+z(C2+z(C3+z(C4+z(C5+z(C6+z*C7)))))))
56		#elif LDBL_MANT_DIG == 113
57		/*
58		* ld128 version of __cos.c. See __cos.c for most comments.
59		*/
60		/*
61		* Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]:
62		* \|cos(x) - c(x))\| < 2**-122.0
63		*
64		* 113-bit precision requires more care than 64-bit precision, since
65		* simple methods give a minimax polynomial with coefficient for x^2
66		* that is 1 ulp below 0.5, but we want it to be precisely 0.5. See
67		* above for more details.
68		*/
69		static const long double
70		C1 = 0.04166666666666666666666666666666658424671L,
71		C2 = -0.001388888888888888888888888888863490893732L,
72		C3 = 0.00002480158730158730158730158600795304914210L,
73		C4 = -0.2755731922398589065255474947078934284324e-6L,
74		C5 = 0.2087675698786809897659225313136400793948e-8L,
75		C6 = -0.1147074559772972315817149986812031204775e-10L,
76		C7 = 0.4779477332386808976875457937252120293400e-13L;
77		static const double
78		C8 = -0.1561920696721507929516718307820958119868e-15,
79		C9 = 0.4110317413744594971475941557607804508039e-18,
80		C10 = -0.8896592467191938803288521958313920156409e-21,
81		C11 = 0.1601061435794535138244346256065192782581e-23;
82		#define POLY(z) (z(C1+z(C2+z(C3+z(C4+z(C5+z(C6+z*(C7+ \
83		z(C8+z(C9+z(C10+zC11)))))))))))
84		#endif
85
86		long double __cosl(long double x, long double y)
87		{
88		long double hz,z,r,w;
89
90		z = x*x;
91		r = POLY(z);
92		hz = 0.5*z;
93		w = 1.0-hz;
94		return w + (((1.0-w)-hz) + (zr-xy));
95		}
96		#endif

lib/libc/musl/src/math/__sinl.c deleted-78

...	...	@@ -1,78 +0,0 @@
1		/* origin: FreeBSD /usr/src/lib/msun/ld80/k_sinl.c */
2		/* origin: FreeBSD /usr/src/lib/msun/ld128/k_sinl.c */
3		/*
4		* ====================================================
5		* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6		* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
7		*
8		* Developed at SunSoft, a Sun Microsystems, Inc. business.
9		* Permission to use, copy, modify, and distribute this
10		* software is freely granted, provided that this notice
11		* is preserved.
12		* ====================================================
13		*/
14
15		#include "libm.h"
16
17		#if (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
18		#if LDBL_MANT_DIG == 64
19		/*
20		* ld80 version of __sin.c. See __sin.c for most comments.
21		*/
22		/*
23		* Domain [-0.7854, 0.7854], range ~[-1.89e-22, 1.915e-22]
24		* \|sin(x)/x - s(x)\| < 2**-72.1
25		*
26		* See __cosl.c for more details about the polynomial.
27		*/
28		static const long double
29		S1 = -0.166666666666666666671L; /* -0xaaaaaaaaaaaaaaab.0p-66 */
30		static const double
31		S2 = 0.0083333333333333332, /* 0x11111111111111.0p-59 */
32		S3 = -0.00019841269841269427, /* -0x1a01a01a019f81.0p-65 */
33		S4 = 0.0000027557319223597490, /* 0x171de3a55560f7.0p-71 */
34		S5 = -0.000000025052108218074604, /* -0x1ae64564f16cad.0p-78 */
35		S6 = 1.6059006598854211e-10, /* 0x161242b90243b5.0p-85 */
36		S7 = -7.6429779983024564e-13, /* -0x1ae42ebd1b2e00.0p-93 */
37		S8 = 2.6174587166648325e-15; /* 0x179372ea0b3f64.0p-101 */
38		#define POLY(z) (S2+z(S3+z(S4+z(S5+z(S6+z(S7+zS8))))))
39		#elif LDBL_MANT_DIG == 113
40		/*
41		* ld128 version of __sin.c. See __sin.c for most comments.
42		*/
43		/*
44		* Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
45		* \|sin(x)/x - s(x)\| < 2**-122.1
46		*
47		* See __cosl.c for more details about the polynomial.
48		*/
49		static const long double
50		S1 = -0.16666666666666666666666666666666666606732416116558L,
51		S2 = 0.0083333333333333333333333333333331135404851288270047L,
52		S3 = -0.00019841269841269841269841269839935785325638310428717L,
53		S4 = 0.27557319223985890652557316053039946268333231205686e-5L,
54		S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
55		S6 = 0.16059043836821614596571832194524392581082444805729e-9L,
56		S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
57		S8 = 0.28114572543451292625024967174638477283187397621303e-14L;
58		static const double
59		S9 = -0.82206352458348947812512122163446202498005154296863e-17,
60		S10 = 0.19572940011906109418080609928334380560135358385256e-19,
61		S11 = -0.38680813379701966970673724299207480965452616911420e-22,
62		S12 = 0.64038150078671872796678569586315881020659912139412e-25;
63		#define POLY(z) (S2+z(S3+z(S4+z(S5+z(S6+z(S7+z(S8+ \
64		z(S9+z(S10+z(S11+zS12))))))))))
65		#endif
66
67		long double __sinl(long double x, long double y, int iy)
68		{
69		long double z,r,v;
70
71		z = x*x;
72		v = z*x;
73		r = POLY(z);
74		if (iy == 0)
75		return x+v(S1+zr);
76		return x-((z(0.5y-vr)-y)-vS1);
77		}
78		#endif

lib/libc/musl/src/math/__tanl.c deleted-143

...	...	@@ -1,143 +0,0 @@
1		/* origin: FreeBSD /usr/src/lib/msun/ld80/k_tanl.c */
2		/* origin: FreeBSD /usr/src/lib/msun/ld128/k_tanl.c */
3		/*
4		* ====================================================
5		* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
6		* Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
7		*
8		* Permission to use, copy, modify, and distribute this
9		* software is freely granted, provided that this notice
10		* is preserved.
11		* ====================================================
12		*/
13
14		#include "libm.h"
15
16		#if (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
17		#if LDBL_MANT_DIG == 64
18		/*
19		* ld80 version of __tan.c. See __tan.c for most comments.
20		*/
21		/*
22		* Domain [-0.67434, 0.67434], range ~[-2.25e-22, 1.921e-22]
23		* \|tan(x)/x - t(x)\| < 2**-71.9
24		*
25		* See __cosl.c for more details about the polynomial.
26		*/
27		static const long double
28		T3 = 0.333333333333333333180L, /* 0xaaaaaaaaaaaaaaa5.0p-65 */
29		T5 = 0.133333333333333372290L, /* 0x88888888888893c3.0p-66 */
30		T7 = 0.0539682539682504975744L, /* 0xdd0dd0dd0dc13ba2.0p-68 */
31		pio4 = 0.785398163397448309628L, /* 0xc90fdaa22168c235.0p-64 */
32		pio4lo = -1.25413940316708300586e-20L; /* -0xece675d1fc8f8cbb.0p-130 */
33		static const double
34		T9 = 0.021869488536312216, /* 0x1664f4882cc1c2.0p-58 */
35		T11 = 0.0088632355256619590, /* 0x1226e355c17612.0p-59 */
36		T13 = 0.0035921281113786528, /* 0x1d6d3d185d7ff8.0p-61 */
37		T15 = 0.0014558334756312418, /* 0x17da354aa3f96b.0p-62 */
38		T17 = 0.00059003538700862256, /* 0x13559358685b83.0p-63 */
39		T19 = 0.00023907843576635544, /* 0x1f56242026b5be.0p-65 */
40		T21 = 0.000097154625656538905, /* 0x1977efc26806f4.0p-66 */
41		T23 = 0.000038440165747303162, /* 0x14275a09b3ceac.0p-67 */
42		T25 = 0.000018082171885432524, /* 0x12f5e563e5487e.0p-68 */
43		T27 = 0.0000024196006108814377, /* 0x144c0d80cc6896.0p-71 */
44		T29 = 0.0000078293456938132840, /* 0x106b59141a6cb3.0p-69 */
45		T31 = -0.0000032609076735050182, /* -0x1b5abef3ba4b59.0p-71 */
46		T33 = 0.0000023261313142559411; /* 0x13835436c0c87f.0p-71 */
47		#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
48		w * (T25 + w * (T29 + w * T33)))))))
49		#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
50		w * (T27 + w * T31))))))
51		#elif LDBL_MANT_DIG == 113
52		/*
53		* ld128 version of __tan.c. See __tan.c for most comments.
54		*/
55		/*
56		* Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37]
57		* \|tan(x)/x - t(x)\| < 2**-117.8 (XXX should be ~1e-37)
58		*
59		* See __cosl.c for more details about the polynomial.
60		*/
61		static const long double
62		T3 = 0x1.5555555555555555555555555553p-2L,
63		T5 = 0x1.1111111111111111111111111eb5p-3L,
64		T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L,
65		T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L,
66		T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L,
67		T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L,
68		T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L,
69		T17 = 0x1.355824803674477dfcf726649efep-11L,
70		T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L,
71		T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L,
72		T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L,
73		T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L,
74		T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L,
75		T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L,
76		T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L,
77		T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L,
78		T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L,
79		T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L,
80		pio4 = 0x1.921fb54442d18469898cc51701b8p-1L,
81		pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L;
82		static const double
83		T39 = 0.000000028443389121318352, /* 0x1e8a7592977938.0p-78 */
84		T41 = 0.000000011981013102001973, /* 0x19baa1b1223219.0p-79 */
85		T43 = 0.0000000038303578044958070, /* 0x107385dfb24529.0p-80 */
86		T45 = 0.0000000034664378216909893, /* 0x1dc6c702a05262.0p-81 */
87		T47 = -0.0000000015090641701997785, /* -0x19ecef3569ebb6.0p-82 */
88		T49 = 0.0000000029449552300483952, /* 0x194c0668da786a.0p-81 */
89		T51 = -0.0000000022006995706097711, /* -0x12e763b8845268.0p-81 */
90		T53 = 0.0000000015468200913196612, /* 0x1a92fc98c29554.0p-82 */
91		T55 = -0.00000000061311613386849674, /* -0x151106cbc779a9.0p-83 */
92		T57 = 1.4912469681508012e-10; /* 0x147edbdba6f43a.0p-85 */
93		#define RPOLY(w) (T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 + \
94		w * (T25 + w * (T29 + w * (T33 + w * (T37 + w * (T41 + \
95		w * (T45 + w * (T49 + w * (T53 + w * T57)))))))))))))
96		#define VPOLY(w) (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 + \
97		w * (T27 + w * (T31 + w * (T35 + w * (T39 + w * (T43 + \
98		w * (T47 + w * (T51 + w * T55))))))))))))
99		#endif
100
101		long double __tanl(long double x, long double y, int odd) {
102		long double z, r, v, w, s, a, t;
103		int big, sign;
104
105		big = fabsl(x) >= 0.67434;
106		if (big) {
107		sign = 0;
108		if (x < 0) {
109		sign = 1;
110		x = -x;
111		y = -y;
112		}
113		x = (pio4 - x) + (pio4lo - y);
114		y = 0.0;
115		}
116		z = x * x;
117		w = z * z;
118		r = RPOLY(w);
119		v = z * VPOLY(w);
120		s = z * x;
121		r = y + z * (s * (r + v) + y) + T3 * s;
122		w = x + r;
123		if (big) {
124		s = 1 - 2*odd;
125		v = s - 2.0 * (x + (r - w * w / (w + s)));
126		return sign ? -v : v;
127		}
128		if (!odd)
129		return w;
130		/*
131		* if allow error up to 2 ulp, simply return
132		* -1.0 / (x+r) here
133		*/
134		/* compute -1.0 / (x+r) accurately */
135		z = w;
136		z = z + 0x1p32 - 0x1p32;
137		v = r - (z - x); /* z+v = r+x */
138		t = a = -1.0 / w; /* a = -1.0/w */
139		t = t + 0x1p32 - 0x1p32;
140		s = 1.0 + t * z;
141		return t + a * (s + t * v);
142		}
143		#endif

lib/libc/musl/src/math/aarch64/lrint.c deleted-10

...	...	@@ -1,10 +0,0 @@
1		#include <math.h>
2
3		long lrint(double x)
4		{
5		long n;
6		__asm__ (
7		"frintx %d1, %d1\n"
8		"fcvtzs %x0, %d1\n" : "=r"(n), "+w"(x));
9		return n;
10		}

lib/libc/musl/src/math/aarch64/lrintf.c deleted-10

...	...	@@ -1,10 +0,0 @@
1		#include <math.h>
2
3		long lrintf(float x)
4		{
5		long n;
6		__asm__ (
7		"frintx %s1, %s1\n"
8		"fcvtzs %x0, %s1\n" : "=r"(n), "+w"(x));
9		return n;
10		}

lib/libc/musl/src/math/aarch64/rintf.c deleted-7

...	...	@@ -1,7 +0,0 @@
1		#include <math.h>
2
3		float rintf(float x)
4		{
5		__asm__ ("frintx %s0, %s1" : "=w"(x) : "w"(x));
6		return x;
7		}

lib/libc/musl/src/math/cosl.c deleted-39

...	...	@@ -1,39 +0,0 @@
1		#include "libm.h"
2
3		#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4		long double cosl(long double x) {
5		return cos(x);
6		}
7		#elif (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
8		long double cosl(long double x)
9		{
10		union ldshape u = {x};
11		unsigned n;
12		long double y[2], hi, lo;
13
14		u.i.se &= 0x7fff;
15		if (u.i.se == 0x7fff)
16		return x - x;
17		x = u.f;
18		if (x < M_PI_4) {
19		if (u.i.se < 0x3fff - LDBL_MANT_DIG)
20		/* raise inexact if x!=0 */
21		return 1.0 + x;
22		return __cosl(x, 0);
23		}
24		n = __rem_pio2l(x, y);
25		hi = y[0];
26		lo = y[1];
27		switch (n & 3) {
28		case 0:
29		return __cosl(hi, lo);
30		case 1:
31		return -__sinl(hi, lo, 1);
32		case 2:
33		return -__cosl(hi, lo);
34		case 3:
35		default:
36		return __sinl(hi, lo, 1);
37		}
38		}
39		#endif

lib/libc/musl/src/math/finite.c deleted-7

...	...	@@ -1,7 +0,0 @@
1		#define _GNU_SOURCE
2		#include <math.h>
3
4		int finite(double x)
5		{
6		return isfinite(x);
7		}

lib/libc/musl/src/math/finitef.c deleted-7

...	...	@@ -1,7 +0,0 @@
1		#define _GNU_SOURCE
2		#include <math.h>
3
4		int finitef(float x)
5		{
6		return isfinite(x);
7		}

lib/libc/musl/src/math/frexp.c deleted-23

...	...	@@ -1,23 +0,0 @@
1		#include <math.h>
2		#include <stdint.h>
3
4		double frexp(double x, int *e)
5		{
6		union { double d; uint64_t i; } y = { x };
7		int ee = y.i>>52 & 0x7ff;
8
9		if (!ee) {
10		if (x) {
11		x = frexp(x*0x1p64, e);
12		*e -= 64;
13		} else *e = 0;
14		return x;
15		} else if (ee == 0x7ff) {
16		return x;
17		}
18
19		*e = ee - 0x3fe;
20		y.i &= 0x800fffffffffffffull;
21		y.i \|= 0x3fe0000000000000ull;
22		return y.d;
23		}

lib/libc/musl/src/math/frexpf.c deleted-23

...	...	@@ -1,23 +0,0 @@
1		#include <math.h>
2		#include <stdint.h>
3
4		float frexpf(float x, int *e)
5		{
6		union { float f; uint32_t i; } y = { x };
7		int ee = y.i>>23 & 0xff;
8
9		if (!ee) {
10		if (x) {
11		x = frexpf(x*0x1p64, e);
12		*e -= 64;
13		} else *e = 0;
14		return x;
15		} else if (ee == 0xff) {
16		return x;
17		}
18
19		*e = ee - 0x7e;
20		y.i &= 0x807ffffful;
21		y.i \|= 0x3f000000ul;
22		return y.f;
23		}

lib/libc/musl/src/math/frexpl.c deleted-29

...	...	@@ -1,29 +0,0 @@
1		#include "libm.h"
2
3		#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4		long double frexpl(long double x, int *e)
5		{
6		return frexp(x, e);
7		}
8		#elif (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
9		long double frexpl(long double x, int *e)
10		{
11		union ldshape u = {x};
12		int ee = u.i.se & 0x7fff;
13
14		if (!ee) {
15		if (x) {
16		x = frexpl(x*0x1p120, e);
17		*e -= 120;
18		} else *e = 0;
19		return x;
20		} else if (ee == 0x7fff) {
21		return x;
22		}
23
24		*e = ee - 0x3ffe;
25		u.i.se &= 0x8000;
26		u.i.se \|= 0x3ffe;
27		return u.f;
28		}
29		#endif

lib/libc/musl/src/math/i386/lrint.c deleted-8

...	...	@@ -1,8 +0,0 @@
1		#include <math.h>
2
3		long lrint(double x)
4		{
5		long r;
6		__asm__ ("fistpl %0" : "=m"(r) : "t"(x) : "st");
7		return r;
8		}

lib/libc/musl/src/math/i386/lrintf.c deleted-8

...	...	@@ -1,8 +0,0 @@
1		#include <math.h>
2
3		long lrintf(float x)
4		{
5		long r;
6		__asm__ ("fistpl %0" : "=m"(r) : "t"(x) : "st");
7		return r;
8		}

lib/libc/musl/src/math/i386/rintf.c deleted-7

...	...	@@ -1,7 +0,0 @@
1		#include <math.h>
2
3		float rintf(float x)
4		{
5		__asm__ ("frndint" : "+t"(x));
6		return x;
7		}

lib/libc/musl/src/math/lrint.c deleted-72

...	...	@@ -1,72 +0,0 @@
1		#include <limits.h>
2		#include <fenv.h>
3		#include <math.h>
4		#include "libm.h"
5
6		/*
7		If the result cannot be represented (overflow, nan), then
8		lrint raises the invalid exception.
9
10		Otherwise if the input was not an integer then the inexact
11		exception is raised.
12
13		C99 is a bit vague about whether inexact exception is
14		allowed to be raised when invalid is raised.
15		(F.9 explicitly allows spurious inexact exceptions, F.9.6.5
16		does not make it clear if that rule applies to lrint, but
17		IEEE 754r 7.8 seems to forbid spurious inexact exception in
18		the ineger conversion functions)
19
20		So we try to make sure that no spurious inexact exception is
21		raised in case of an overflow.
22
23		If the bit size of long > precision of double, then there
24		cannot be inexact rounding in case the result overflows,
25		otherwise LONG_MAX and LONG_MIN can be represented exactly
26		as a double.
27		*/
28
29		#if LONG_MAX < 1U<<53 && defined(FE_INEXACT)
30		#include <float.h>
31		#include <stdint.h>
32		#if FLT_EVAL_METHOD==0 \|\| FLT_EVAL_METHOD==1
33		#define EPS DBL_EPSILON
34		#elif FLT_EVAL_METHOD==2
35		#define EPS LDBL_EPSILON
36		#endif
37		#ifdef __GNUC__
38		/* avoid stack frame in lrint */
39		__attribute__((noinline))
40		#endif
41		static long lrint_slow(double x)
42		{
43		#pragma STDC FENV_ACCESS ON
44		int e;
45
46		e = fetestexcept(FE_INEXACT);
47		x = rint(x);
48		if (!e && (x > LONG_MAX \|\| x < LONG_MIN))
49		feclearexcept(FE_INEXACT);
50		/* conversion */
51		return x;
52		}
53
54		long lrint(double x)
55		{
56		uint32_t abstop = asuint64(x)>>32 & 0x7fffffff;
57		uint64_t sign = asuint64(x) & (1ULL << 63);
58
59		if (abstop < 0x41dfffff) {
60		/* \|x\| < 0x7ffffc00, no overflow */
61		double_t toint = asdouble(asuint64(1/EPS) \| sign);
62		double_t y = x + toint - toint;
63		return (long)y;
64		}
65		return lrint_slow(x);
66		}
67		#else
68		long lrint(double x)
69		{
70		return rint(x);
71		}
72		#endif

lib/libc/musl/src/math/lrintf.c deleted-8

...	...	@@ -1,8 +0,0 @@
1		#include <math.h>
2
3		/* uses LONG_MAX > 2^24, see comments in lrint.c */
4
5		long lrintf(float x)
6		{
7		return rintf(x);
8		}

lib/libc/musl/src/math/powerpc64/lrint.c deleted-16

...	...	@@ -1,16 +0,0 @@
1		#include <math.h>
2
3		#ifdef _ARCH_PWR5X
4
5		long lrint(double x)
6		{
7		long n;
8		__asm__ ("fctid %0, %1" : "=d"(n) : "d"(x));
9		return n;
10		}
11
12		#else
13
14		#include "../lrint.c"
15
16		#endif

lib/libc/musl/src/math/powerpc64/lrintf.c deleted-16

...	...	@@ -1,16 +0,0 @@
1		#include <math.h>
2
3		#ifdef _ARCH_PWR5X
4
5		long lrintf(float x)
6		{
7		long n;
8		__asm__ ("fctid %0, %1" : "=d"(n) : "f"(x));
9		return n;
10		}
11
12		#else
13
14		#include "../lrintf.c"
15
16		#endif

lib/libc/musl/src/math/rintf.c deleted-30

...	...	@@ -1,30 +0,0 @@
1		#include <float.h>
2		#include <math.h>
3		#include <stdint.h>
4
5		#if FLT_EVAL_METHOD==0
6		#define EPS FLT_EPSILON
7		#elif FLT_EVAL_METHOD==1
8		#define EPS DBL_EPSILON
9		#elif FLT_EVAL_METHOD==2
10		#define EPS LDBL_EPSILON
11		#endif
12		static const float_t toint = 1/EPS;
13
14		float rintf(float x)
15		{
16		union {float f; uint32_t i;} u = {x};
17		int e = u.i>>23 & 0xff;
18		int s = u.i>>31;
19		float_t y;
20
21		if (e >= 0x7f+23)
22		return x;
23		if (s)
24		y = x - toint + toint;
25		else
26		y = x + toint - toint;
27		if (y == 0)
28		return s ? -0.0f : 0.0f;
29		return y;
30		}

lib/libc/musl/src/math/s390x/rintf.c deleted-15

...	...	@@ -1,15 +0,0 @@
1		#include <math.h>
2
3		#if defined(__HTM__) \|\| __ARCH__ >= 9
4
5		float rintf(float x)
6		{
7		__asm__ ("fiebr %0, 0, %1" : "=f"(x) : "f"(x));
8		return x;
9		}
10
11		#else
12
13		#include "../rintf.c"
14
15		#endif

lib/libc/musl/src/math/sincosl.c deleted-60

...	...	@@ -1,60 +0,0 @@
1		#define _GNU_SOURCE
2		#include "libm.h"
3
4		#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
5		void sincosl(long double x, long double sin, long double cos)
6		{
7		double sind, cosd;
8		sincos(x, &sind, &cosd);
9		*sin = sind;
10		*cos = cosd;
11		}
12		#elif (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
13		void sincosl(long double x, long double sin, long double cos)
14		{
15		union ldshape u = {x};
16		unsigned n;
17		long double y[2], s, c;
18
19		u.i.se &= 0x7fff;
20		if (u.i.se == 0x7fff) {
21		sin = cos = x - x;
22		return;
23		}
24		if (u.f < M_PI_4) {
25		if (u.i.se < 0x3fff - LDBL_MANT_DIG) {
26		/* raise underflow if subnormal */
27		if (u.i.se == 0) FORCE_EVAL(x*0x1p-120f);
28		*sin = x;
29		/* raise inexact if x!=0 */
30		*cos = 1.0 + x;
31		return;
32		}
33		*sin = __sinl(x, 0, 0);
34		*cos = __cosl(x, 0);
35		return;
36		}
37		n = __rem_pio2l(x, y);
38		s = __sinl(y[0], y[1], 1);
39		c = __cosl(y[0], y[1]);
40		switch (n & 3) {
41		case 0:
42		*sin = s;
43		*cos = c;
44		break;
45		case 1:
46		*sin = c;
47		*cos = -s;
48		break;
49		case 2:
50		*sin = -s;
51		*cos = -c;
52		break;
53		case 3:
54		default:
55		*sin = -c;
56		*cos = s;
57		break;
58		}
59		}
60		#endif

lib/libc/musl/src/math/sinl.c deleted-41

...	...	@@ -1,41 +0,0 @@
1		#include "libm.h"
2
3		#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4		long double sinl(long double x)
5		{
6		return sin(x);
7		}
8		#elif (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
9		long double sinl(long double x)
10		{
11		union ldshape u = {x};
12		unsigned n;
13		long double y[2], hi, lo;
14
15		u.i.se &= 0x7fff;
16		if (u.i.se == 0x7fff)
17		return x - x;
18		if (u.f < M_PI_4) {
19		if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
20		/* raise inexact if x!=0 and underflow if subnormal */
21		FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
22		return x;
23		}
24		return __sinl(x, 0.0, 0);
25		}
26		n = __rem_pio2l(x, y);
27		hi = y[0];
28		lo = y[1];
29		switch (n & 3) {
30		case 0:
31		return __sinl(hi, lo, 1);
32		case 1:
33		return __cosl(hi, lo);
34		case 2:
35		return -__sinl(hi, lo, 1);
36		case 3:
37		default:
38		return -__cosl(hi, lo);
39		}
40		}
41		#endif

lib/libc/musl/src/math/tanl.c deleted-29

...	...	@@ -1,29 +0,0 @@
1		#include "libm.h"
2
3		#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
4		long double tanl(long double x)
5		{
6		return tan(x);
7		}
8		#elif (LDBL_MANT_DIG == 64 \|\| LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
9		long double tanl(long double x)
10		{
11		union ldshape u = {x};
12		long double y[2];
13		unsigned n;
14
15		u.i.se &= 0x7fff;
16		if (u.i.se == 0x7fff)
17		return x - x;
18		if (u.f < M_PI_4) {
19		if (u.i.se < 0x3fff - LDBL_MANT_DIG/2) {
20		/* raise inexact if x!=0 and underflow if subnormal */
21		FORCE_EVAL(u.i.se == 0 ? x*0x1p-120f : x+0x1p120f);
22		return x;
23		}
24		return __tanl(x, 0, 0);
25		}
26		n = __rem_pio2l(x, y);
27		return __tanl(y[0], y[1], n&1);
28		}
29		#endif

lib/libc/musl/src/math/x32/lrint.s deleted-5

...	...	@@ -1,5 +0,0 @@
1		.global lrint
2		.type lrint,@function
3		lrint:
4		cvtsd2si %xmm0,%rax
5		ret

lib/libc/musl/src/math/x32/lrintf.s deleted-5

...	...	@@ -1,5 +0,0 @@
1		.global lrintf
2		.type lrintf,@function
3		lrintf:
4		cvtss2si %xmm0,%rax
5		ret

lib/libc/musl/src/math/x86_64/lrint.c deleted-8

...	...	@@ -1,8 +0,0 @@
1		#include <math.h>
2
3		long lrint(double x)
4		{
5		long r;
6		__asm__ ("cvtsd2si %1, %0" : "=r"(r) : "x"(x));
7		return r;
8		}

lib/libc/musl/src/math/x86_64/lrintf.c deleted-8

...	...	@@ -1,8 +0,0 @@
1		#include <math.h>
2
3		long lrintf(float x)
4		{
5		long r;
6		__asm__ ("cvtss2si %1, %0" : "=r"(r) : "x"(x));
7		return r;
8		}

src/libs/mingw.zig+1-15

...	...	@@ -613,8 +613,6 @@ const mingw32_generic_src = [_][]const u8{
613	613	"math" ++ path.sep_str ++ "fpclassify.c",
614	614	"math" ++ path.sep_str ++ "fpclassifyf.c",
615	615	"math" ++ path.sep_str ++ "fpclassifyl.c",
616		"math" ++ path.sep_str ++ "frexpf.c",
617		"math" ++ path.sep_str ++ "frexpl.c",
618	616	"math" ++ path.sep_str ++ "ldexpf.c",
619	617	"math" ++ path.sep_str ++ "lgamma.c",
620	618	"math" ++ path.sep_str ++ "lgammaf.c",
...	...	@@ -942,8 +940,6 @@ const mingw32_x86_src = [_][]const u8{
942	940	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atan2l.c",
943	941	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atanhl.c",
944	942	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "atanl.c",
945		"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cosl.c",
946		"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cosl_internal.S",
947	943	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "cossinl.c",
948	944	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "exp2l.S",
949	945	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "expl.c",
...	...	@@ -965,9 +961,6 @@ const mingw32_x86_src = [_][]const u8{
965	961	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbn.S",
966	962	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnf.S",
967	963	"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "scalbnl.S",
968		"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl.c",
969		"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "sinl_internal.S",
970		"math" ++ path.sep_str ++ "x86" ++ path.sep_str ++ "tanl.S",
971	964	// ucrtbase
972	965	"math" ++ path.sep_str ++ "nextafterl.c",
973	966	"math" ++ path.sep_str ++ "nexttoward.c",
...	...	@@ -987,22 +980,15 @@ const mingw32_x86_32_src = [_][]const u8{
987	980	const mingw32_arm_src = [_][]const u8{
988	981	// mingwex
989	982	"math" ++ path.sep_str ++ "arm-common" ++ path.sep_str ++ "ldexpl.c",
990		"math" ++ path.sep_str ++ "arm-common" ++ path.sep_str ++ "sincosl.c",
991	983	};
992	984
993	985	const mingw32_arm32_src = [_][]const u8{
994	986	// mingwex
995	987	"math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "s_rint.c",
996	988	"math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "s_rintf.c",
997		"math" ++ path.sep_str ++ "arm" ++ path.sep_str ++ "sincos.S",
998	989	};
999	990
1000		const mingw32_arm64_src = [_][]const u8{
1001		// mingwex
1002		"math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "rint.c",
1003		"math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "rintf.c",
1004		"math" ++ path.sep_str ++ "arm64" ++ path.sep_str ++ "sincos.S",
1005		};
	991	const mingw32_arm64_src = [_][]const u8{};
1006	992
1007	993	const mingw32_winpthreads_src = [_][]const u8{
1008	994	// winpthreads

src/libs/musl.zig-28

...	...	@@ -787,13 +787,10 @@ const src_files = [_][]const u8{
787	787	"musl/src/math/aarch64/llrintf.c",
788	788	"musl/src/math/aarch64/llround.c",
789	789	"musl/src/math/aarch64/llroundf.c",
790		"musl/src/math/aarch64/lrint.c",
791		"musl/src/math/aarch64/lrintf.c",
792	790	"musl/src/math/aarch64/lround.c",
793	791	"musl/src/math/aarch64/lroundf.c",
794	792	"musl/src/math/aarch64/nearbyint.c",
795	793	"musl/src/math/aarch64/nearbyintf.c",
796		"musl/src/math/aarch64/rintf.c",
797	794	"musl/src/math/acosh.c",
798	795	"musl/src/math/acoshl.c",
799	796	"musl/src/math/acosl.c",
...	...	@@ -814,8 +811,6 @@ const src_files = [_][]const u8{
814	811	"musl/src/math/__cos.c",
815	812	"musl/src/math/__cosdf.c",
816	813	"musl/src/math/coshl.c",
817		"musl/src/math/__cosl.c",
818		"musl/src/math/cosl.c",
819	814	"musl/src/math/erf.c",
820	815	"musl/src/math/erff.c",
821	816	"musl/src/math/erfl.c",
...	...	@@ -831,17 +826,12 @@ const src_files = [_][]const u8{
831	826	"musl/src/math/__expo2f.c",
832	827	"musl/src/math/fdimf.c",
833	828	"musl/src/math/fdiml.c",
834		"musl/src/math/finite.c",
835		"musl/src/math/finitef.c",
836	829	"musl/src/math/fma.c",
837	830	"musl/src/math/fmaf.c",
838	831	"musl/src/math/fmal.c",
839	832	"musl/src/math/__fpclassify.c",
840	833	"musl/src/math/__fpclassifyf.c",
841	834	"musl/src/math/__fpclassifyl.c",
842		"musl/src/math/frexp.c",
843		"musl/src/math/frexpf.c",
844		"musl/src/math/frexpl.c",
845	835	"musl/src/math/i386/acosl.s",
846	836	"musl/src/math/i386/asinf.s",
847	837	"musl/src/math/i386/asinl.s",
...	...	@@ -865,8 +855,6 @@ const src_files = [_][]const u8{
865	855	"musl/src/math/i386/log1p.s",
866	856	"musl/src/math/i386/log2l.s",
867	857	"musl/src/math/i386/logl.s",
868		"musl/src/math/i386/lrint.c",
869		"musl/src/math/i386/lrintf.c",
870	858	"musl/src/math/i386/lrintl.c",
871	859	"musl/src/math/i386/remainder.c",
872	860	"musl/src/math/i386/remainderf.c",
...	...	@@ -874,7 +862,6 @@ const src_files = [_][]const u8{
874	862	"musl/src/math/i386/remquof.s",
875	863	"musl/src/math/i386/remquol.s",
876	864	"musl/src/math/i386/remquo.s",
877		"musl/src/math/i386/rintf.c",
878	865	"musl/src/math/i386/rintl.c",
879	866	"musl/src/math/i386/scalblnf.s",
880	867	"musl/src/math/i386/scalblnl.s",
...	...	@@ -915,8 +902,6 @@ const src_files = [_][]const u8{
915	902	"musl/src/math/logbf.c",
916	903	"musl/src/math/logbl.c",
917	904	"musl/src/math/logl.c",
918		"musl/src/math/lrint.c",
919		"musl/src/math/lrintf.c",
920	905	"musl/src/math/lrintl.c",
921	906	"musl/src/math/lround.c",
922	907	"musl/src/math/lroundf.c",
...	...	@@ -945,8 +930,6 @@ const src_files = [_][]const u8{
945	930	"musl/src/math/pow_data.c",
946	931	"musl/src/math/powerpc64/fma.c",
947	932	"musl/src/math/powerpc64/fmaf.c",
948		"musl/src/math/powerpc64/lrint.c",
949		"musl/src/math/powerpc64/lrintf.c",
950	933	"musl/src/math/powerpc64/lround.c",
951	934	"musl/src/math/powerpc64/lroundf.c",
952	935	"musl/src/math/powerpc/fma.c",
...	...	@@ -964,7 +947,6 @@ const src_files = [_][]const u8{
964	947	"musl/src/math/remquo.c",
965	948	"musl/src/math/remquof.c",
966	949	"musl/src/math/remquol.c",
967		"musl/src/math/rintf.c",
968	950	"musl/src/math/rintl.c",
969	951	"musl/src/math/riscv32/fma.c",
970	952	"musl/src/math/riscv32/fmaf.c",
...	...	@@ -975,7 +957,6 @@ const src_files = [_][]const u8{
975	957	"musl/src/math/s390x/nearbyint.c",
976	958	"musl/src/math/s390x/nearbyintf.c",
977	959	"musl/src/math/s390x/nearbyintl.c",
978		"musl/src/math/s390x/rintf.c",
979	960	"musl/src/math/s390x/rintl.c",
980	961	"musl/src/math/scalb.c",
981	962	"musl/src/math/scalbf.c",
...	...	@@ -992,18 +973,13 @@ const src_files = [_][]const u8{
992	973	"musl/src/math/significand.c",
993	974	"musl/src/math/significandf.c",
994	975	"musl/src/math/__sin.c",
995		"musl/src/math/sincosl.c",
996	976	"musl/src/math/__sindf.c",
997	977	"musl/src/math/sinh.c",
998	978	"musl/src/math/sinhf.c",
999	979	"musl/src/math/sinhl.c",
1000		"musl/src/math/__sinl.c",
1001		"musl/src/math/sinl.c",
1002	980	"musl/src/math/__tan.c",
1003	981	"musl/src/math/__tandf.c",
1004	982	"musl/src/math/tanhl.c",
1005		"musl/src/math/__tanl.c",
1006		"musl/src/math/tanl.c",
1007	983	"musl/src/math/tgamma.c",
1008	984	"musl/src/math/tgammaf.c",
1009	985	"musl/src/math/tgammal.c",
...	...	@@ -1023,9 +999,7 @@ const src_files = [_][]const u8{
1023	999	"musl/src/math/x32/log1pl.s",
1024	1000	"musl/src/math/x32/log2l.s",
1025	1001	"musl/src/math/x32/logl.s",
1026		"musl/src/math/x32/lrintf.s",
1027	1002	"musl/src/math/x32/lrintl.s",
1028		"musl/src/math/x32/lrint.s",
1029	1003	"musl/src/math/x32/remainderl.s",
1030	1004	"musl/src/math/x32/rintl.s",
1031	1005	"musl/src/math/x86_64/acosl.s",
...	...	@@ -1044,8 +1018,6 @@ const src_files = [_][]const u8{
1044	1018	"musl/src/math/x86_64/log1pl.s",
1045	1019	"musl/src/math/x86_64/log2l.s",
1046	1020	"musl/src/math/x86_64/logl.s",
1047		"musl/src/math/x86_64/lrint.c",
1048		"musl/src/math/x86_64/lrintf.c",
1049	1021	"musl/src/math/x86_64/lrintl.c",
1050	1022	"musl/src/math/x86_64/remainderl.c",
1051	1023	"musl/src/math/x86_64/remquol.c",

src/libs/wasi_libc.zig-14

...	...	@@ -682,8 +682,6 @@ const libc_top_half_src_files = [_][]const u8{
682	682	"musl/src/math/__cos.c",
683	683	"musl/src/math/__cosdf.c",
684	684	"musl/src/math/coshl.c",
685		"musl/src/math/__cosl.c",
686		"musl/src/math/cosl.c",
687	685	"musl/src/math/erf.c",
688	686	"musl/src/math/erff.c",
689	687	"musl/src/math/erfl.c",
...	...	@@ -697,13 +695,8 @@ const libc_top_half_src_files = [_][]const u8{
697	695	"musl/src/math/expm1l.c",
698	696	"musl/src/math/fdimf.c",
699	697	"musl/src/math/fdiml.c",
700		"musl/src/math/finite.c",
701		"musl/src/math/finitef.c",
702	698	"musl/src/math/fma.c",
703	699	"musl/src/math/fmaf.c",
704		"musl/src/math/frexp.c",
705		"musl/src/math/frexpf.c",
706		"musl/src/math/frexpl.c",
707	700	"musl/src/math/ilogb.c",
708	701	"musl/src/math/ilogbf.c",
709	702	"musl/src/math/ilogbl.c",
...	...	@@ -737,8 +730,6 @@ const libc_top_half_src_files = [_][]const u8{
737	730	"musl/src/math/logbf.c",
738	731	"musl/src/math/logbl.c",
739	732	"musl/src/math/logl.c",
740		"musl/src/math/lrint.c",
741		"musl/src/math/lrintf.c",
742	733	"musl/src/math/lrintl.c",
743	734	"musl/src/math/lround.c",
744	735	"musl/src/math/lroundf.c",
...	...	@@ -785,16 +776,11 @@ const libc_top_half_src_files = [_][]const u8{
785	776	"musl/src/math/significand.c",
786	777	"musl/src/math/significandf.c",
787	778	"musl/src/math/__sin.c",
788		"musl/src/math/sincosl.c",
789	779	"musl/src/math/__sindf.c",
790	780	"musl/src/math/sinhl.c",
791		"musl/src/math/__sinl.c",
792		"musl/src/math/sinl.c",
793	781	"musl/src/math/__tan.c",
794	782	"musl/src/math/__tandf.c",
795	783	"musl/src/math/tanhl.c",
796		"musl/src/math/__tanl.c",
797		"musl/src/math/tanl.c",
798	784	"musl/src/math/tgamma.c",
799	785	"musl/src/math/tgammaf.c",
800	786	"musl/src/math/tgammal.c",

test/libc.zig+1-1

...	...	@@ -293,7 +293,7 @@ pub fn addCases(cases: *tests.LibcContext) void {
293	293	cases.addLibcTestCase("math/remquo.c", true, .{});
294	294	cases.addLibcTestCase("math/remquof.c", true, .{});
295	295	cases.addLibcTestCase("math/remquol.c", true, .{});
296		// cases.addLibcTestCase("math/rint.c", true, .{});
	296	cases.addLibcTestCase("math/rint.c", true, .{});
297	297	cases.addLibcTestCase("math/rintf.c", true, .{});
298	298	// cases.addLibcTestCase("math/rintl.c", true, .{});
299	299	cases.addLibcTestCase("math/round.c", true, .{});