authorgravatar for lewis.gaul@gmail.comLewis Gaul <lewis.gaul@gmail.com> 2024-05-25 22:50:56+01:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2025-07-17 05:04:59+02:00
log7abb170f59ed959c4004ec37258c8905903f36b1
tree2c629e39b88558f38499fbed73c522ce5481f14d
parentf34b26231ee5c6d4810b333617940bbd682e1a36

Add tests for math.expm1(), fixing bug in 32-bit function

1 files changed, 68 insertions(+), 29 deletions(-)

lib/std/math/expm1.zig+68-29
...@@ -10,6 +10,7 @@ const std = @import("../std.zig");...@@ -10,6 +10,7 @@ const std = @import("../std.zig");
10const math = std.math;10const math = std.math;
11const mem = std.mem;11const mem = std.mem;
12const expect = std.testing.expect;12const expect = std.testing.expect;
13const expectEqual = std.testing.expectEqual;
1314
14/// Returns e raised to the power of x, minus 1 (e^x - 1). This is more accurate than exp(e, x) - 115/// Returns e raised to the power of x, minus 1 (e^x - 1). This is more accurate than exp(e, x) - 1
15/// when x is near 0.16/// when x is near 0.
...@@ -39,9 +40,9 @@ fn expm1_32(x_: f32) f32 {...@@ -39,9 +40,9 @@ fn expm1_32(x_: f32) f32 {
39 const Q2: f32 = 1.5807170421e-3;40 const Q2: f32 = 1.5807170421e-3;
4041
41 var x = x_;42 var x = x_;
42 const ux = @as(u32, @bitCast(x));43 const ux: u32 = @bitCast(x);
43 const hx = ux & 0x7FFFFFFF;44 const hx = ux & 0x7FFFFFFF;
44 const sign = hx >> 31;45 const sign = ux >> 31;
4546
46 // TODO: Shouldn't need this check explicitly.47 // TODO: Shouldn't need this check explicitly.
47 if (math.isNegativeInf(x)) {48 if (math.isNegativeInf(x)) {
...@@ -147,7 +148,7 @@ fn expm1_32(x_: f32) f32 {...@@ -147,7 +148,7 @@ fn expm1_32(x_: f32) f32 {
147 return y - 1.0;148 return y - 1.0;
148 }149 }
149150
150 const uf = @as(f32, @bitCast(@as(u32, @intCast(0x7F -% k)) << 23));151 const uf: f32 = @bitCast(@as(u32, @intCast(0x7F -% k)) << 23);
151 if (k < 23) {152 if (k < 23) {
152 return (x - e + (1 - uf)) * twopk;153 return (x - e + (1 - uf)) * twopk;
153 } else {154 } else {
...@@ -286,39 +287,77 @@ fn expm1_64(x_: f64) f64 {...@@ -286,39 +287,77 @@ fn expm1_64(x_: f64) f64 {
286 }287 }
287}288}
288289
289test expm1 {290test "expm1_32() special" {
290 try expect(expm1(@as(f32, 0.0)) == expm1_32(0.0));291 try expectEqual(expm1_32(0.0), 0.0);
291 try expect(expm1(@as(f64, 0.0)) == expm1_64(0.0));292 try expectEqual(expm1_32(-0.0), 0.0);
293 try expectEqual(expm1_32(math.ln2), 1.0);
294 try expectEqual(expm1_32(math.inf(f32)), math.inf(f32));
295 try expectEqual(expm1_32(-math.inf(f32)), -1.0);
296 try expect(math.isNan(expm1_32(math.nan(f32))));
297 try expect(math.isNan(expm1_32(math.snan(f32))));
292}298}
293299
294test expm1_32 {300test "expm1_32() sanity" {
295 const epsilon = 0.000001;301 try expectEqual(expm1_32(-0x1.0223a0p+3), -0x1.ffd6e0p-1);
296302 try expectEqual(expm1_32(0x1.161868p+2), 0x1.30712ap+6);
297 try expect(math.isPositiveZero(expm1_32(0.0)));303 try expectEqual(expm1_32(-0x1.0c34b4p+3), -0x1.ffe1fap-1);
298 try expect(math.approxEqAbs(f32, expm1_32(0.0), 0.0, epsilon));304 try expectEqual(expm1_32(-0x1.a206f0p+2), -0x1.ff4116p-1);
299 try expect(math.approxEqAbs(f32, expm1_32(0.2), 0.221403, epsilon));305 try expectEqual(expm1_32(0x1.288bbcp+3), 0x1.4ab480p+13); // Disagrees with GCC in last bit
300 try expect(math.approxEqAbs(f32, expm1_32(0.8923), 1.440737, epsilon));306 try expectEqual(expm1_32(0x1.52efd0p-1), 0x1.e09536p-1);
301 try expect(math.approxEqAbs(f32, expm1_32(1.5), 3.481689, epsilon));307 try expectEqual(expm1_32(-0x1.a05cc8p-2), -0x1.561c3ep-2);
308 try expectEqual(expm1_32(0x1.1f9efap-1), 0x1.81ec4ep-1);
309 try expectEqual(expm1_32(0x1.8c5db0p-1), 0x1.2b3364p+0);
310 try expectEqual(expm1_32(-0x1.5b86eap-1), -0x1.f8951ap-2);
302}311}
303312
304test expm1_64 {313test "expm1_32() boundary" {
305 const epsilon = 0.000001;314 // TODO: The last value before inf is actually 0x1.62e300p+6 -> 0x1.ff681ep+127
315 // try expectEqual(expm1_32(0x1.62e42ep+6), 0x1.ffff08p+127); // Last value before result is inf
316 try expectEqual(expm1_32(0x1.62e430p+6), math.inf(f32)); // First value that gives inf
317 try expectEqual(expm1_32(0x1.fffffep+127), math.inf(f32)); // Max input value
318 try expectEqual(expm1_32(0x1p-149), 0x1p-149); // Min positive input value
319 try expectEqual(expm1_32(-0x1p-149), -0x1p-149); // Min negative input value
320 try expectEqual(expm1_32(0x1p-126), 0x1p-126); // First positive subnormal input
321 try expectEqual(expm1_32(-0x1p-126), -0x1p-126); // First negative subnormal input
322 try expectEqual(expm1_32(0x1.fffffep-125), 0x1.fffffep-125); // Last positive value before subnormal
323 try expectEqual(expm1_32(-0x1.fffffep-125), -0x1.fffffep-125); // Last negative value before subnormal
324 try expectEqual(expm1_32(-0x1.154244p+4), -0x1.fffffep-1); // Last value before result is -1
325 try expectEqual(expm1_32(-0x1.154246p+4), -1); // First value where result is -1
326}
306327
307 try expect(math.isPositiveZero(expm1_64(0.0)));328test "expm1_64() special" {
308 try expect(math.approxEqAbs(f64, expm1_64(0.0), 0.0, epsilon));329 try expectEqual(expm1_64(0.0), 0.0);
309 try expect(math.approxEqAbs(f64, expm1_64(0.2), 0.221403, epsilon));330 try expectEqual(expm1_64(-0.0), 0.0);
310 try expect(math.approxEqAbs(f64, expm1_64(0.8923), 1.440737, epsilon));331 try expectEqual(expm1_64(math.ln2), 1.0);
311 try expect(math.approxEqAbs(f64, expm1_64(1.5), 3.481689, epsilon));332 try expectEqual(expm1_64(math.inf(f64)), math.inf(f64));
333 try expectEqual(expm1_64(-math.inf(f64)), -1.0);
334 try expect(math.isNan(expm1_64(math.nan(f64))));
335 try expect(math.isNan(expm1_64(math.snan(f64))));
312}336}
313337
314test "expm1_32.special" {338test "expm1_64() sanity" {
315 try expect(math.isPositiveInf(expm1_32(math.inf(f32))));339 try expectEqual(expm1_64(-0x1.02239f3c6a8f1p+3), -0x1.ffd6df9b02b3ep-1);
316 try expect(expm1_32(-math.inf(f32)) == -1.0);340 try expectEqual(expm1_64(0x1.161868e18bc67p+2), 0x1.30712ed238c04p+6);
317 try expect(math.isNan(expm1_32(math.nan(f32))));341 try expectEqual(expm1_64(-0x1.0c34b3e01e6e7p+3), -0x1.ffe1f94e493e7p-1);
342 try expectEqual(expm1_64(-0x1.a206f0a19dcc4p+2), -0x1.ff4115c03f78dp-1);
343 try expectEqual(expm1_64(0x1.288bbb0d6a1e6p+3), 0x1.4ab477496e07ep+13);
344 try expectEqual(expm1_64(0x1.52efd0cd80497p-1), 0x1.e095382100a01p-1);
345 try expectEqual(expm1_64(-0x1.a05cc754481d1p-2), -0x1.561c3e0582be6p-2);
346 try expectEqual(expm1_64(0x1.1f9ef934745cbp-1), 0x1.81ec4cd4d4a8fp-1);
347 try expectEqual(expm1_64(0x1.8c5db097f7442p-1), 0x1.2b3363a944bf7p+0);
348 try expectEqual(expm1_64(-0x1.5b86ea8118a0ep-1), -0x1.f8951aebffbafp-2);
318}349}
319350
320test "expm1_64.special" {351test "expm1_64() boundary" {
321 try expect(math.isPositiveInf(expm1_64(math.inf(f64))));352 try expectEqual(expm1_64(0x1.62e42fefa39efp+9), 0x1.fffffffffff2ap+1023); // Last value before result is inf
322 try expect(expm1_64(-math.inf(f64)) == -1.0);353 try expectEqual(expm1_64(0x1.62e42fefa39f0p+9), math.inf(f64)); // First value that gives inf
323 try expect(math.isNan(expm1_64(math.nan(f64))));354 try expectEqual(expm1_64(0x1.fffffffffffffp+1023), math.inf(f64)); // Max input value
355 try expectEqual(expm1_64(0x1p-1074), 0x1p-1074); // Min positive input value
356 try expectEqual(expm1_64(-0x1p-1074), -0x1p-1074); // Min negative input value
357 try expectEqual(expm1_64(0x1p-1022), 0x1p-1022); // First positive subnormal input
358 try expectEqual(expm1_64(-0x1p-1022), -0x1p-1022); // First negative subnormal input
359 try expectEqual(expm1_64(0x1.fffffffffffffp-1021), 0x1.fffffffffffffp-1021); // Last positive value before subnormal
360 try expectEqual(expm1_64(-0x1.fffffffffffffp-1021), -0x1.fffffffffffffp-1021); // Last negative value before subnormal
361 try expectEqual(expm1_64(-0x1.2b708872320e1p+5), -0x1.fffffffffffffp-1); // Last value before result is -1
362 try expectEqual(expm1_64(-0x1.2b708872320e2p+5), -1); // First value where result is -1
324}363}