Commit f890de62948a0134a8d703dd2f1ee37701ec4c00 | zig

author	Lewis Gaul <lewis.gaul@gmail.com> 2021-10-26 23:57:58+01:00
committer	GitHub <noreply@github.com> 2021-10-26 18:57:58-04:00
log	f890de62948a0134a8d703dd2f1ee37701ec4c00
tree	5151dadec62428299befc0d21f3f1464c4402d0a
parent	ad5b90ab1092cb50b38e4a1e6bbed463a38ef7cf
signature	Signed by PGP key 4AEE18F83AFDEB23

Fix bug in exp2() (#9999)

* Fix bug in exp2_64 to handle negative values (bad translation from C)

* Apply fix to exp2_32() as well, and modify comment on musl behaviour

* Use +%= instead of @addWithOverflow()

1 files changed, 21 insertions(+), 11 deletions(-)

lib
- std
  - math
    exp2.zig

lib/std/math/exp2.zig+21-11

...	...	@@ -78,13 +78,15 @@ fn exp2_32(x: f32) f32 {
78	78	return 1.0 + x;
79	79	}
80	80
	81	// NOTE: musl relies on unsafe behaviours which are replicated below
	82	// (addition/bit-shift overflow). Appears that this produces the
	83	// intended result but should confirm how GCC/Clang handle this to ensure.
	84
81	85	var uf = x + redux;
82	86	var i_0 = @bitCast(u32, uf);
83		i_0 += tblsiz / 2;
	87	i_0 +%= tblsiz / 2;
84	88
85	89	const k = i_0 / tblsiz;
86		// NOTE: musl relies on undefined overflow shift behaviour. Appears that this produces the
87		// intended result but should confirm how GCC/Clang handle this to ensure.
88	90	const uk = @bitCast(f64, @as(u64, 0x3FF + k) << 52);
89	91	i_0 &= tblsiz - 1;
90	92	uf -= redux;
...	...	@@ -357,7 +359,7 @@ const exp2dt = [_]f64{
357	359	};
358	360
359	361	fn exp2_64(x: f64) f64 {
360		const tblsiz = @intCast(u32, exp2dt.len / 2);
	362	const tblsiz: u32 = @intCast(u32, exp2dt.len / 2);
361	363	const redux: f64 = 0x1.8p52 / @intToFloat(f64, tblsiz);
362	364	const P1: f64 = 0x1.62e42fefa39efp-1;
363	365	const P2: f64 = 0x1.ebfbdff82c575p-3;
...	...	@@ -400,22 +402,27 @@ fn exp2_64(x: f64) f64 {
400	402	return 1.0 + x;
401	403	}
402	404
	405	// NOTE: musl relies on unsafe behaviours which are replicated below
	406	// (addition overflow, division truncation, casting). Appears that this
	407	// produces the intended result but should confirm how GCC/Clang handle this
	408	// to ensure.
	409
403	410	// reduce x
404		var uf = x + redux;
	411	var uf: f64 = x + redux;
405	412	// NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here
406		var i_0 = @truncate(u32, @bitCast(u64, uf));
407		i_0 += tblsiz / 2;
	413	var i_0: u32 = @truncate(u32, @bitCast(u64, uf));
	414	i_0 +%= tblsiz / 2;
408	415
409	416	const k: u32 = i_0 / tblsiz * tblsiz;
410		const ik = @bitCast(i32, k / tblsiz);
	417	const ik: i32 = @divTrunc(@bitCast(i32, k), tblsiz);
411	418	i_0 %= tblsiz;
412	419	uf -= redux;
413	420
414	421	// r = exp2(y) = exp2t[i_0] * p(z - eps[i])
415		var z = x - uf;
416		const t = exp2dt[@intCast(usize, 2 * i_0)];
	422	var z: f64 = x - uf;
	423	const t: f64 = exp2dt[@intCast(usize, 2 * i_0)];
417	424	z -= exp2dt[@intCast(usize, 2 * i_0 + 1)];
418		const r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
	425	const r: f64 = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
419	426
420	427	return math.scalbn(r, ik);
421	428	}
...	...	@@ -433,6 +440,7 @@ test "math.exp2_32" {
433	440	try expect(math.approxEqAbs(f32, exp2_32(0.8923), 1.856133, epsilon));
434	441	try expect(math.approxEqAbs(f32, exp2_32(1.5), 2.828427, epsilon));
435	442	try expect(math.approxEqAbs(f32, exp2_32(37.45), 187747237888, epsilon));
	443	try expect(math.approxEqAbs(f32, exp2_32(-1), 0.5, epsilon));
436	444	}
437	445
438	446	test "math.exp2_64" {
...	...	@@ -442,6 +450,8 @@ test "math.exp2_64" {
442	450	try expect(math.approxEqAbs(f64, exp2_64(0.2), 1.148698, epsilon));
443	451	try expect(math.approxEqAbs(f64, exp2_64(0.8923), 1.856133, epsilon));
444	452	try expect(math.approxEqAbs(f64, exp2_64(1.5), 2.828427, epsilon));
	453	try expect(math.approxEqAbs(f64, exp2_64(-1), 0.5, epsilon));
	454	try expect(math.approxEqAbs(f64, exp2_64(-0x1.a05cc754481d1p-2), 0x1.824056efc687cp-1, epsilon));
445	455	}
446	456
447	457	test "math.exp2_32.special" {