Commit 70b6b98e91fb9d1d5575da4d59c5c523864a8f0e | zig

author	alice <hi@dnydxn.moe> 2022-05-17 21:55:22+01:00
committer	alice <hi@dnydxn.moe> 2022-05-17 21:55:22+01:00
log	70b6b98e91fb9d1d5575da4d59c5c523864a8f0e
tree	a728394816c426f6e4b2c1af33fc4ee3e49add14
parent	ceeec8d19f36bf0a2cb132d35c0e119c4fec476d
signature	Commit is signed but in an unrecognized format.

Simplify `Copysign`

9 files changed, 51 insertions(+), 118 deletions(-)

lib
- compiler_rt
  - fma.zig
- std
  - math
    atanh.zig
    complex
    cosh.zig
    proj.zig
    sinh.zig
    sqrt.zig
    tanh.zig
    copysign.zig
    pow.zig

lib/compiler_rt/fma.zig+2-2

...	...	@@ -57,7 +57,7 @@ pub fn fma(x: f64, y: f64, z: f64) callconv(.C) f64 {
57	57	if (spread <= 53 * 2) {
58	58	zs = math.scalbn(zs, -spread);
59	59	} else {
60		zs = math.copysign(f64, math.floatMin(f64), zs);
	60	zs = math.copysign(math.floatMin(f64), zs);
61	61	}
62	62
63	63	const xy = dd_mul(xs, ys);
...	...	@@ -116,7 +116,7 @@ pub fn fmaq(x: f128, y: f128, z: f128) callconv(.C) f128 {
116	116	if (spread <= 113 * 2) {
117	117	zs = math.scalbn(zs, -spread);
118	118	} else {
119		zs = math.copysign(f128, math.floatMin(f128), zs);
	119	zs = math.copysign(math.floatMin(f128), zs);
120	120	}
121	121
122	122	const xy = dd_mul128(xs, ys);

lib/std/math/atanh.zig+2-2

...	...	@@ -33,7 +33,7 @@ fn atanh_32(x: f32) f32 {
33	33	var y = @bitCast(f32, i); // \|x\|
34	34
35	35	if (y == 1.0) {
36		return math.copysign(f32, math.inf(f32), x);
	36	return math.copysign(math.inf(f32), x);
37	37	}
38	38
39	39	if (u < 0x3F800000 - (1 << 23)) {
...	...	@@ -62,7 +62,7 @@ fn atanh_64(x: f64) f64 {
62	62	var y = @bitCast(f64, u & (maxInt(u64) >> 1)); // \|x\|
63	63
64	64	if (y == 1.0) {
65		return math.copysign(f64, math.inf(f64), x);
	65	return math.copysign(math.inf(f64), x);
66	66	}
67	67
68	68	if (e < 0x3FF - 1) {

lib/std/math/complex/cosh.zig+10-10

...	...	@@ -45,13 +45,13 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
45	45	if (ix < 0x42b17218) {
46	46	// x < 88.7: exp(\|x\|) won't overflow
47	47	const h = @exp(@fabs(x)) * 0.5;
48		return Complex(f32).init(math.copysign(f32, h, x) * @cos(y), h * @sin(y));
	48	return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
49	49	}
50	50	// x < 192.7: scale to avoid overflow
51	51	else if (ix < 0x4340b1e7) {
52	52	const v = Complex(f32).init(@fabs(x), y);
53	53	const r = ldexp_cexp(v, -1);
54		return Complex(f32).init(r.re, r.im * math.copysign(f32, 1, x));
	54	return Complex(f32).init(r.re, r.im * math.copysign(@as(f32, 1.0), x));
55	55	}
56	56	// x >= 192.7: result always overflows
57	57	else {
...	...	@@ -61,14 +61,14 @@ fn cosh32(z: Complex(f32)) Complex(f32) {
61	61	}
62	62
63	63	if (ix == 0 and iy >= 0x7f800000) {
64		return Complex(f32).init(y - y, math.copysign(f32, 0, x * (y - y)));
	64	return Complex(f32).init(y - y, math.copysign(@as(f32, 0.0), x * (y - y)));
65	65	}
66	66
67	67	if (iy == 0 and ix >= 0x7f800000) {
68	68	if (hx & 0x7fffff == 0) {
69		return Complex(f32).init(x * x, math.copysign(f32, 0, x) * y);
	69	return Complex(f32).init(x * x, math.copysign(@as(f32, 0.0), x) * y);
70	70	}
71		return Complex(f32).init(x, math.copysign(f32, 0, (x + x) * y));
	71	return Complex(f32).init(x, math.copysign(@as(f32, 0.0), (x + x) * y));
72	72	}
73	73
74	74	if (ix < 0x7f800000 and iy >= 0x7f800000) {
...	...	@@ -113,13 +113,13 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
113	113	if (ix < 0x40862e42) {
114	114	// x < 710: exp(\|x\|) won't overflow
115	115	const h = @exp(@fabs(x)) * 0.5;
116		return Complex(f64).init(h * @cos(y), math.copysign(f64, h, x) * @sin(y));
	116	return Complex(f64).init(h * @cos(y), math.copysign(h, x) * @sin(y));
117	117	}
118	118	// x < 1455: scale to avoid overflow
119	119	else if (ix < 0x4096bbaa) {
120	120	const v = Complex(f64).init(@fabs(x), y);
121	121	const r = ldexp_cexp(v, -1);
122		return Complex(f64).init(r.re, r.im * math.copysign(f64, 1, x));
	122	return Complex(f64).init(r.re, r.im * math.copysign(@as(f64, 1.0), x));
123	123	}
124	124	// x >= 1455: result always overflows
125	125	else {
...	...	@@ -129,14 +129,14 @@ fn cosh64(z: Complex(f64)) Complex(f64) {
129	129	}
130	130
131	131	if (ix \| lx == 0 and iy >= 0x7ff00000) {
132		return Complex(f64).init(y - y, math.copysign(f64, 0, x * (y - y)));
	132	return Complex(f64).init(y - y, math.copysign(@as(f64, 0.0), x * (y - y)));
133	133	}
134	134
135	135	if (iy \| ly == 0 and ix >= 0x7ff00000) {
136	136	if ((hx & 0xfffff) \| lx == 0) {
137		return Complex(f64).init(x * x, math.copysign(f64, 0, x) * y);
	137	return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), x) * y);
138	138	}
139		return Complex(f64).init(x * x, math.copysign(f64, 0, (x + x) * y));
	139	return Complex(f64).init(x * x, math.copysign(@as(f64, 0.0), (x + x) * y));
140	140	}
141	141
142	142	if (ix < 0x7ff00000 and iy >= 0x7ff00000) {

lib/std/math/complex/proj.zig+1-1

...	...	@@ -9,7 +9,7 @@ pub fn proj(z: anytype) Complex(@TypeOf(z.re)) {
9	9	const T = @TypeOf(z.re);
10	10
11	11	if (math.isInf(z.re) or math.isInf(z.im)) {
12		return Complex(T).init(math.inf(T), math.copysign(T, 0, z.re));
	12	return Complex(T).init(math.inf(T), math.copysign(@as(T, 0.0), z.re));
13	13	}
14	14
15	15	return Complex(T).init(z.re, z.im);

lib/std/math/complex/sinh.zig+8-8

...	...	@@ -45,13 +45,13 @@ fn sinh32(z: Complex(f32)) Complex(f32) {
45	45	if (ix < 0x42b17218) {
46	46	// x < 88.7: exp(\|x\|) won't overflow
47	47	const h = @exp(@fabs(x)) * 0.5;
48		return Complex(f32).init(math.copysign(f32, h, x) * @cos(y), h * @sin(y));
	48	return Complex(f32).init(math.copysign(h, x) * @cos(y), h * @sin(y));
49	49	}
50	50	// x < 192.7: scale to avoid overflow
51	51	else if (ix < 0x4340b1e7) {
52	52	const v = Complex(f32).init(@fabs(x), y);
53	53	const r = ldexp_cexp(v, -1);
54		return Complex(f32).init(r.re * math.copysign(f32, 1, x), r.im);
	54	return Complex(f32).init(r.re * math.copysign(@as(f32, 1.0), x), r.im);
55	55	}
56	56	// x >= 192.7: result always overflows
57	57	else {
...	...	@@ -61,14 +61,14 @@ fn sinh32(z: Complex(f32)) Complex(f32) {
61	61	}
62	62
63	63	if (ix == 0 and iy >= 0x7f800000) {
64		return Complex(f32).init(math.copysign(f32, 0, x * (y - y)), y - y);
	64	return Complex(f32).init(math.copysign(@as(f32, 0.0), x * (y - y)), y - y);
65	65	}
66	66
67	67	if (iy == 0 and ix >= 0x7f800000) {
68	68	if (hx & 0x7fffff == 0) {
69	69	return Complex(f32).init(x, y);
70	70	}
71		return Complex(f32).init(x, math.copysign(f32, 0, y));
	71	return Complex(f32).init(x, math.copysign(@as(f32, 0.0), y));
72	72	}
73	73
74	74	if (ix < 0x7f800000 and iy >= 0x7f800000) {
...	...	@@ -112,13 +112,13 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
112	112	if (ix < 0x40862e42) {
113	113	// x < 710: exp(\|x\|) won't overflow
114	114	const h = @exp(@fabs(x)) * 0.5;
115		return Complex(f64).init(math.copysign(f64, h, x) * @cos(y), h * @sin(y));
	115	return Complex(f64).init(math.copysign(h, x) * @cos(y), h * @sin(y));
116	116	}
117	117	// x < 1455: scale to avoid overflow
118	118	else if (ix < 0x4096bbaa) {
119	119	const v = Complex(f64).init(@fabs(x), y);
120	120	const r = ldexp_cexp(v, -1);
121		return Complex(f64).init(r.re * math.copysign(f64, 1, x), r.im);
	121	return Complex(f64).init(r.re * math.copysign(@as(f64, 1.0), x), r.im);
122	122	}
123	123	// x >= 1455: result always overflows
124	124	else {
...	...	@@ -128,14 +128,14 @@ fn sinh64(z: Complex(f64)) Complex(f64) {
128	128	}
129	129
130	130	if (ix \| lx == 0 and iy >= 0x7ff00000) {
131		return Complex(f64).init(math.copysign(f64, 0, x * (y - y)), y - y);
	131	return Complex(f64).init(math.copysign(@as(f64, 0.0), x * (y - y)), y - y);
132	132	}
133	133
134	134	if (iy \| ly == 0 and ix >= 0x7ff00000) {
135	135	if ((hx & 0xfffff) \| lx == 0) {
136	136	return Complex(f64).init(x, y);
137	137	}
138		return Complex(f64).init(x, math.copysign(f64, 0, y));
	138	return Complex(f64).init(x, math.copysign(@as(f64, 0.0), y));
139	139	}
140	140
141	141	if (ix < 0x7ff00000 and iy >= 0x7ff00000) {

lib/std/math/complex/sqrt.zig+6-6

...	...	@@ -43,9 +43,9 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
43	43	// sqrt(-inf + i nan) = nan +- inf i
44	44	// sqrt(-inf + iy) = 0 + inf i
45	45	if (math.signbit(x)) {
46		return Complex(f32).init(@fabs(x - y), math.copysign(f32, x, y));
	46	return Complex(f32).init(@fabs(x - y), math.copysign(x, y));
47	47	} else {
48		return Complex(f32).init(x, math.copysign(f32, y - y, y));
	48	return Complex(f32).init(x, math.copysign(y - y, y));
49	49	}
50	50	}
51	51
...	...	@@ -65,7 +65,7 @@ fn sqrt32(z: Complex(f32)) Complex(f32) {
65	65	const t = @sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5);
66	66	return Complex(f32).init(
67	67	@floatCast(f32, @fabs(y) / (2.0 * t)),
68		@floatCast(f32, math.copysign(f64, t, y)),
	68	@floatCast(f32, math.copysign(t, y)),
69	69	);
70	70	}
71	71	}
...	...	@@ -94,9 +94,9 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
94	94	// sqrt(-inf + i nan) = nan +- inf i
95	95	// sqrt(-inf + iy) = 0 + inf i
96	96	if (math.signbit(x)) {
97		return Complex(f64).init(@fabs(x - y), math.copysign(f64, x, y));
	97	return Complex(f64).init(@fabs(x - y), math.copysign(x, y));
98	98	} else {
99		return Complex(f64).init(x, math.copysign(f64, y - y, y));
	99	return Complex(f64).init(x, math.copysign(y - y, y));
100	100	}
101	101	}
102	102
...	...	@@ -116,7 +116,7 @@ fn sqrt64(z: Complex(f64)) Complex(f64) {
116	116	result = Complex(f64).init(t, y / (2.0 * t));
117	117	} else {
118	118	const t = @sqrt((-x + math.hypot(f64, x, y)) * 0.5);
119		result = Complex(f64).init(@fabs(y) / (2.0 * t), math.copysign(f64, t, y));
	119	result = Complex(f64).init(@fabs(y) / (2.0 * t), math.copysign(t, y));
120	120	}
121	121
122	122	if (scale) {

lib/std/math/complex/tanh.zig+4-4

...	...	@@ -34,7 +34,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) {
34	34	}
35	35	const xx = @bitCast(f32, hx - 0x40000000);
36	36	const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
37		return Complex(f32).init(xx, math.copysign(f32, 0, r));
	37	return Complex(f32).init(xx, math.copysign(@as(f32, 0.0), r));
38	38	}
39	39
40	40	if (!math.isFinite(y)) {
...	...	@@ -45,7 +45,7 @@ fn tanh32(z: Complex(f32)) Complex(f32) {
45	45	// x >= 11
46	46	if (ix >= 0x41300000) {
47	47	const exp_mx = @exp(-@fabs(x));
48		return Complex(f32).init(math.copysign(f32, 1, x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
	48	return Complex(f32).init(math.copysign(@as(f32, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
49	49	}
50	50
51	51	// Kahan's algorithm
...	...	@@ -77,7 +77,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
77	77
78	78	const xx = @bitCast(f64, (@as(u64, hx - 0x40000000) << 32) \| lx);
79	79	const r = if (math.isInf(y)) y else @sin(y) * @cos(y);
80		return Complex(f64).init(xx, math.copysign(f64, 0, r));
	80	return Complex(f64).init(xx, math.copysign(@as(f64, 0.0), r));
81	81	}
82	82
83	83	if (!math.isFinite(y)) {
...	...	@@ -88,7 +88,7 @@ fn tanh64(z: Complex(f64)) Complex(f64) {
88	88	// x >= 22
89	89	if (ix >= 0x40360000) {
90	90	const exp_mx = @exp(-@fabs(x));
91		return Complex(f64).init(math.copysign(f64, 1, x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
	91	return Complex(f64).init(math.copysign(@as(f64, 1.0), x), 4 * @sin(y) * @cos(y) * exp_mx * exp_mx);
92	92	}
93	93
94	94	// Kahan's algorithm

lib/std/math/copysign.zig+17-84

...	...	@@ -1,92 +1,25 @@
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/copysignf.c
5		// https://git.musl-libc.org/cgit/musl/tree/src/math/copysign.c
6
7	1	const std = @import("../std.zig");
8	2	const math = std.math;
9	3	const expect = std.testing.expect;
10		const maxInt = std.math.maxInt;
11
12		/// Returns a value with the magnitude of x and the sign of y.
13		pub fn copysign(comptime T: type, x: T, y: T) T {
14		return switch (T) {
15		f16 => copysign16(x, y),
16		f32 => copysign32(x, y),
17		f64 => copysign64(x, y),
18		f128 => copysign128(x, y),
19		else => @compileError("copysign not implemented for " ++ @typeName(T)),
20		};
21		}
22	4
23		fn copysign16(x: f16, y: f16) f16 {
24		const ux = @bitCast(u16, x);
25		const uy = @bitCast(u16, y);
26
27		const h1 = ux & (maxInt(u16) / 2);
28		const h2 = uy & (@as(u16, 1) << 15);
29		return @bitCast(f16, h1 \| h2);
30		}
31
32		fn copysign32(x: f32, y: f32) f32 {
33		const ux = @bitCast(u32, x);
34		const uy = @bitCast(u32, y);
35
36		const h1 = ux & (maxInt(u32) / 2);
37		const h2 = uy & (@as(u32, 1) << 31);
38		return @bitCast(f32, h1 \| h2);
39		}
40
41		fn copysign64(x: f64, y: f64) f64 {
42		const ux = @bitCast(u64, x);
43		const uy = @bitCast(u64, y);
44
45		const h1 = ux & (maxInt(u64) / 2);
46		const h2 = uy & (@as(u64, 1) << 63);
47		return @bitCast(f64, h1 \| h2);
48		}
49
50		fn copysign128(x: f128, y: f128) f128 {
51		const ux = @bitCast(u128, x);
52		const uy = @bitCast(u128, y);
53
54		const h1 = ux & (maxInt(u128) / 2);
55		const h2 = uy & (@as(u128, 1) << 127);
56		return @bitCast(f128, h1 \| h2);
	5	/// Returns a value with the magnitude of `magnitude` and the sign of `sign`.
	6	pub fn copysign(magnitude: anytype, sign: @TypeOf(magnitude)) @TypeOf(magnitude) {
	7	const T = @TypeOf(magnitude);
	8	const TBits = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
	9	const sign_bit_mask = @as(TBits, 1) << (@bitSizeOf(T) - 1);
	10	const mag = @bitCast(TBits, magnitude) & ~sign_bit_mask;
	11	const sgn = @bitCast(TBits, sign) & sign_bit_mask;
	12	return @bitCast(T, mag \| sgn);
57	13	}
58	14
59	15	test "math.copysign" {
60		try expect(copysign(f16, 1.0, 1.0) == copysign16(1.0, 1.0));
61		try expect(copysign(f32, 1.0, 1.0) == copysign32(1.0, 1.0));
62		try expect(copysign(f64, 1.0, 1.0) == copysign64(1.0, 1.0));
63		try expect(copysign(f128, 1.0, 1.0) == copysign128(1.0, 1.0));
64		}
65
66		test "math.copysign16" {
67		try expect(copysign16(5.0, 1.0) == 5.0);
68		try expect(copysign16(5.0, -1.0) == -5.0);
69		try expect(copysign16(-5.0, -1.0) == -5.0);
70		try expect(copysign16(-5.0, 1.0) == 5.0);
71		}
72
73		test "math.copysign32" {
74		try expect(copysign32(5.0, 1.0) == 5.0);
75		try expect(copysign32(5.0, -1.0) == -5.0);
76		try expect(copysign32(-5.0, -1.0) == -5.0);
77		try expect(copysign32(-5.0, 1.0) == 5.0);
78		}
79
80		test "math.copysign64" {
81		try expect(copysign64(5.0, 1.0) == 5.0);
82		try expect(copysign64(5.0, -1.0) == -5.0);
83		try expect(copysign64(-5.0, -1.0) == -5.0);
84		try expect(copysign64(-5.0, 1.0) == 5.0);
85		}
86
87		test "math.copysign128" {
88		try expect(copysign128(5.0, 1.0) == 5.0);
89		try expect(copysign128(5.0, -1.0) == -5.0);
90		try expect(copysign128(-5.0, -1.0) == -5.0);
91		try expect(copysign128(-5.0, 1.0) == 5.0);
	16	inline for ([_]type{ f16, f32, f64, f80, f128 }) \|T\| {
	17	try expect(copysign(@as(T, 1.0), @as(T, 1.0)) == 1.0);
	18	try expect(copysign(@as(T, 2.0), @as(T, -2.0)) == -2.0);
	19	try expect(copysign(@as(T, -3.0), @as(T, 3.0)) == 3.0);
	20	try expect(copysign(@as(T, -4.0), @as(T, -4.0)) == -4.0);
	21	try expect(copysign(@as(T, 5.0), @as(T, -500.0)) == -5.0);
	22	try expect(copysign(math.inf(T), @as(T, -0.0)) == -math.inf(T));
	23	try expect(copysign(@as(T, 6.0), -math.nan(T)) == -6.0);
	24	}
92	25	}

lib/std/math/pow.zig+1-1

...	...	@@ -60,7 +60,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {
60	60	if (y < 0) {
61	61	// pow(+-0, y) = +- 0 for y an odd integer
62	62	if (isOddInteger(y)) {
63		return math.copysign(T, math.inf(T), x);
	63	return math.copysign(math.inf(T), x);
64	64	}
65	65	// pow(+-0, y) = +inf for y an even integer
66	66	else {