1 files changed, 75 insertions, 0 deletions
diff --git a/std/math/complex/ldexp.zig b/std/math/complex/ldexp.zig
new file mode 100644
index 0000000000..4fb5a6815f
--- /dev/null
+++ b/std/math/complex/ldexp.zig
@@ -0,0 +1,75 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn ldexp_cexp(z: var, expt: i32) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+
+    return switch (T) {
+        f32 => ldexp_cexp32(z, expt),
+        f64 => ldexp_cexp64(z, expt),
+        else => unreachable,
+    };
+}
+
+fn frexp_exp32(x: f32, expt: &i32) f32 {
+    const k = 235;              // reduction constant
+    const kln2 = 162.88958740;  // k * ln2
+
+    const exp_x = math.exp(x - kln2);
+    const hx = @bitCast(u32, exp_x);
+    *expt = i32(hx >> 23) - (0x7f + 127) + k;
+    return @bitCast(f32, (hx & 0x7fffff) | ((0x7f + 127) << 23));
+}
+
+fn ldexp_cexp32(z: &const Complex(f32), expt: i32) Complex(f32) {
+    var ex_expt: i32 = undefined;
+    const exp_x = frexp_exp32(z.re, &ex_expt);
+    const exptf = expt + ex_expt;
+
+    const half_expt1 = @divTrunc(exptf, 2);
+    const scale1 = @bitCast(f32, (0x7f + half_expt1) << 23);
+
+    const half_expt2 = exptf - half_expt1;
+    const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23);
+
+    return Complex(f32).new(
+        math.cos(z.im) * exp_x * scale1 * scale2,
+        math.sin(z.im) * exp_x * scale1 * scale2,
+    );
+}
+
+fn frexp_exp64(x: f64, expt: &i32) f64 {
+    const k = 1799;                         // reduction constant
+    const kln2 = 1246.97177782734161156;    // k * ln2
+
+    const exp_x = math.exp(x - kln2);
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+
+    *expt = i32(hx >> 20) - (0x3ff + 1023) + k;
+
+    const high_word = (hx & 0xfffff) | ((0x3ff + 1023) << 20);
+    return @bitCast(f64, (u64(high_word) << 32) | lx);
+}
+
+fn ldexp_cexp64(z: &const Complex(f64), expt: i32) Complex(f64) {
+    var ex_expt: i32 = undefined;
+    const exp_x = frexp_exp64(z.re, &ex_expt);
+    const exptf = i64(expt + ex_expt);
+
+    const half_expt1 = @divTrunc(exptf, 2);
+    const scale1 = @bitCast(f64, (0x3ff + half_expt1) << 20);
+
+    const half_expt2 = exptf - half_expt1;
+    const scale2 = @bitCast(f64, (0x3ff + half_expt2) << 20);
+
+    return Complex(f64).new(
+        math.cos(z.im) * exp_x * scale1 * scale2,
+        math.sin(z.im) * exp_x * scale1 * scale2,
+    );
+}