std: Fix complex ldexp implementation

Two bugs in the implementation ported from musl made all the complex
functions relying on ldexp return incorrect results in some cases.

Spotted in #9047

2 files changed, 40 insertions, 14 deletions

diff --git a/lib/std/math/complex/exp.zig b/lib/std/math/complex/exp.zig
index af85c37e1d..80fa837ab0 100644
--- a/lib/std/math/complex/exp.zig
+++ b/lib/std/math/complex/exp.zig
@@ -124,20 +124,42 @@ fn exp64(z: Complex(f64)) Complex(f64) {
     }
 }
 
-const epsilon = 0.0001;
-
 test "complex.cexp32" {
-    const a = Complex(f32).init(5, 3);
-    const c = exp(a);
+    const tolerance_f32 = math.sqrt(math.epsilon(f32));
+
+    {
+        const a = Complex(f32).init(5, 3);
+        const c = exp(a);
+
+        try testing.expectApproxEqRel(@as(f32, -1.46927917e+02), c.re, tolerance_f32);
+        try testing.expectApproxEqRel(@as(f32, 2.0944065e+01), c.im, tolerance_f32);
+    }
+
+    {
+        const a = Complex(f32).init(88.8, 0x1p-149);
+        const c = exp(a);
 
-    try testing.expect(math.approxEqAbs(f32, c.re, -146.927917, epsilon));
-    try testing.expect(math.approxEqAbs(f32, c.im, 20.944065, epsilon));
+        try testing.expectApproxEqAbs(math.inf(f32), c.re, tolerance_f32);
+        try testing.expectApproxEqAbs(@as(f32, 5.15088629e-07), c.im, tolerance_f32);
+    }
 }
 
 test "complex.cexp64" {
-    const a = Complex(f64).init(5, 3);
-    const c = exp(a);
+    const tolerance_f64 = math.sqrt(math.epsilon(f64));
 
-    try testing.expect(math.approxEqAbs(f64, c.re, -146.927917, epsilon));
-    try testing.expect(math.approxEqAbs(f64, c.im, 20.944065, epsilon));
+    {
+        const a = Complex(f64).init(5, 3);
+        const c = exp(a);
+
+        try testing.expectApproxEqRel(@as(f64, -1.469279139083189e+02), c.re, tolerance_f64);
+        try testing.expectApproxEqRel(@as(f64, 2.094406620874596e+01), c.im, tolerance_f64);
+    }
+
+    {
+        const a = Complex(f64).init(709.8, 0x1p-1074);
+        const c = exp(a);
+
+        try testing.expectApproxEqAbs(math.inf(f64), c.re, tolerance_f64);
+        try testing.expectApproxEqAbs(@as(f64, 9.036659362159884e-16), c.im, tolerance_f64);
+    }
 }
diff --git a/lib/std/math/complex/ldexp.zig b/lib/std/math/complex/ldexp.zig
index 9f31aae549..51ac105256 100644
--- a/lib/std/math/complex/ldexp.zig
+++ b/lib/std/math/complex/ldexp.zig
@@ -13,6 +13,7 @@ const std = @import("../../std.zig");
 const debug = std.debug;
 const math = std.math;
 const cmath = math.complex;
+const testing = std.testing;
 const Complex = cmath.Complex;
 
 /// Returns exp(z) scaled to avoid overflow.
@@ -48,7 +49,10 @@ fn ldexp_cexp32(z: Complex(f32), expt: i32) Complex(f32) {
     const half_expt2 = exptf - half_expt1;
     const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23);
 
-    return Complex(f32).init(math.cos(z.im) * exp_x * scale1 * scale2, math.sin(z.im) * exp_x * scale1 * scale2);
+    return Complex(f32).init(
+        math.cos(z.im) * exp_x * scale1 * scale2,
+        math.sin(z.im) * exp_x * scale1 * scale2,
+    );
 }
 
 fn frexp_exp64(x: f64, expt: *i32) f64 {
@@ -57,7 +61,7 @@ fn frexp_exp64(x: f64, expt: *i32) f64 {
 
     const exp_x = math.exp(x - kln2);
 
-    const fx = @bitCast(u64, x);
+    const fx = @bitCast(u64, exp_x);
     const hx = @intCast(u32, fx >> 32);
     const lx = @truncate(u32, fx);
 
@@ -73,10 +77,10 @@ fn ldexp_cexp64(z: Complex(f64), expt: i32) Complex(f64) {
     const exptf = @as(i64, expt + ex_expt);
 
     const half_expt1 = @divTrunc(exptf, 2);
-    const scale1 = @bitCast(f64, (0x3ff + half_expt1) << 20);
+    const scale1 = @bitCast(f64, (0x3ff + half_expt1) << (20 + 32));
 
     const half_expt2 = exptf - half_expt1;
-    const scale2 = @bitCast(f64, (0x3ff + half_expt2) << 20);
+    const scale2 = @bitCast(f64, (0x3ff + half_expt2) << (20 + 32));
 
     return Complex(f64).init(
         math.cos(z.im) * exp_x * scale1 * scale2,