1 files changed, 178 insertions, 0 deletions
diff --git a/lib/std/math/frexp.zig b/lib/std/math/frexp.zig
new file mode 100644
index 0000000000..2759cd6492
--- /dev/null
+++ b/lib/std/math/frexp.zig
@@ -0,0 +1,178 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/frexp.c
+
+const std = @import("../std.zig");
+const math = std.math;
+const expect = std.testing.expect;
+
+fn frexp_result(comptime T: type) type {
+    return struct {
+        significand: T,
+        exponent: i32,
+    };
+}
+pub const frexp32_result = frexp_result(f32);
+pub const frexp64_result = frexp_result(f64);
+
+/// Breaks x into a normalized fraction and an integral power of two.
+/// f == frac * 2^exp, with |frac| in the interval [0.5, 1).
+///
+/// Special Cases:
+///  - frexp(+-0)   = +-0, 0
+///  - frexp(+-inf) = +-inf, 0
+///  - frexp(nan)   = nan, undefined
+pub fn frexp(x: var) frexp_result(@typeOf(x)) {
+    const T = @typeOf(x);
+    return switch (T) {
+        f32 => frexp32(x),
+        f64 => frexp64(x),
+        else => @compileError("frexp not implemented for " ++ @typeName(T)),
+    };
+}
+
+fn frexp32(x: f32) frexp32_result {
+    var result: frexp32_result = undefined;
+
+    var y = @bitCast(u32, x);
+    const e = @intCast(i32, y >> 23) & 0xFF;
+
+    if (e == 0) {
+        if (x != 0) {
+            // subnormal
+            result = frexp32(x * 0x1.0p64);
+            result.exponent -= 64;
+        } else {
+            // frexp(+-0) = (+-0, 0)
+            result.significand = x;
+            result.exponent = 0;
+        }
+        return result;
+    } else if (e == 0xFF) {
+        // frexp(nan) = (nan, undefined)
+        result.significand = x;
+        result.exponent = undefined;
+
+        // frexp(+-inf) = (+-inf, 0)
+        if (math.isInf(x)) {
+            result.exponent = 0;
+        }
+
+        return result;
+    }
+
+    result.exponent = e - 0x7E;
+    y &= 0x807FFFFF;
+    y |= 0x3F000000;
+    result.significand = @bitCast(f32, y);
+    return result;
+}
+
+fn frexp64(x: f64) frexp64_result {
+    var result: frexp64_result = undefined;
+
+    var y = @bitCast(u64, x);
+    const e = @intCast(i32, y >> 52) & 0x7FF;
+
+    if (e == 0) {
+        if (x != 0) {
+            // subnormal
+            result = frexp64(x * 0x1.0p64);
+            result.exponent -= 64;
+        } else {
+            // frexp(+-0) = (+-0, 0)
+            result.significand = x;
+            result.exponent = 0;
+        }
+        return result;
+    } else if (e == 0x7FF) {
+        // frexp(nan) = (nan, undefined)
+        result.significand = x;
+        result.exponent = undefined;
+
+        // frexp(+-inf) = (+-inf, 0)
+        if (math.isInf(x)) {
+            result.exponent = 0;
+        }
+
+        return result;
+    }
+
+    result.exponent = e - 0x3FE;
+    y &= 0x800FFFFFFFFFFFFF;
+    y |= 0x3FE0000000000000;
+    result.significand = @bitCast(f64, y);
+    return result;
+}
+
+test "math.frexp" {
+    const a = frexp(f32(1.3));
+    const b = frexp32(1.3);
+    expect(a.significand == b.significand and a.exponent == b.exponent);
+
+    const c = frexp(f64(1.3));
+    const d = frexp64(1.3);
+    expect(c.significand == d.significand and c.exponent == d.exponent);
+}
+
+test "math.frexp32" {
+    const epsilon = 0.000001;
+    var r: frexp32_result = undefined;
+
+    r = frexp32(1.3);
+    expect(math.approxEq(f32, r.significand, 0.65, epsilon) and r.exponent == 1);
+
+    r = frexp32(78.0234);
+    expect(math.approxEq(f32, r.significand, 0.609558, epsilon) and r.exponent == 7);
+}
+
+test "math.frexp64" {
+    const epsilon = 0.000001;
+    var r: frexp64_result = undefined;
+
+    r = frexp64(1.3);
+    expect(math.approxEq(f64, r.significand, 0.65, epsilon) and r.exponent == 1);
+
+    r = frexp64(78.0234);
+    expect(math.approxEq(f64, r.significand, 0.609558, epsilon) and r.exponent == 7);
+}
+
+test "math.frexp32.special" {
+    var r: frexp32_result = undefined;
+
+    r = frexp32(0.0);
+    expect(r.significand == 0.0 and r.exponent == 0);
+
+    r = frexp32(-0.0);
+    expect(r.significand == -0.0 and r.exponent == 0);
+
+    r = frexp32(math.inf(f32));
+    expect(math.isPositiveInf(r.significand) and r.exponent == 0);
+
+    r = frexp32(-math.inf(f32));
+    expect(math.isNegativeInf(r.significand) and r.exponent == 0);
+
+    r = frexp32(math.nan(f32));
+    expect(math.isNan(r.significand));
+}
+
+test "math.frexp64.special" {
+    var r: frexp64_result = undefined;
+
+    r = frexp64(0.0);
+    expect(r.significand == 0.0 and r.exponent == 0);
+
+    r = frexp64(-0.0);
+    expect(r.significand == -0.0 and r.exponent == 0);
+
+    r = frexp64(math.inf(f64));
+    expect(math.isPositiveInf(r.significand) and r.exponent == 0);
+
+    r = frexp64(-math.inf(f64));
+    expect(math.isNegativeInf(r.significand) and r.exponent == 0);
+
+    r = frexp64(math.nan(f64));
+    expect(math.isNan(r.significand));
+}