Commit f32723a237
Changed files (1)
lib
std
math
lib/std/math/frexp.zig
@@ -1,13 +1,8 @@
-// Ported from musl, which is MIT licensed:
-// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
-//
-// https://git.musl-libc.org/cgit/musl/tree/src/math/frexpl.c
-// 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;
+const expectEqual = std.testing.expectEqual;
+const expectApproxEqAbs = std.testing.expectApproxEqAbs;
pub fn Frexp(comptime T: type) type {
return struct {
@@ -24,228 +19,210 @@ pub fn Frexp(comptime T: type) type {
/// - frexp(+-inf) = +-inf, 0
/// - frexp(nan) = nan, undefined
pub fn frexp(x: anytype) Frexp(@TypeOf(x)) {
- const T = @TypeOf(x);
- return switch (T) {
- f32 => frexp32(x),
- f64 => frexp64(x),
- f128 => frexp128(x),
- else => @compileError("frexp not implemented for " ++ @typeName(T)),
- };
-}
-
-// TODO: unify all these implementations using generics
-
-fn frexp32(x: f32) Frexp(f32) {
- var result: Frexp(f32) = undefined;
-
- var y = @as(u32, @bitCast(x));
- const e = @as(i32, @intCast(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;
+ const T: type = @TypeOf(x);
+
+ const bits: comptime_int = @typeInfo(T).Float.bits;
+ const Int: type = std.meta.Int(.unsigned, bits);
+
+ const exp_bits: comptime_int = math.floatExponentBits(T);
+ const mant_bits: comptime_int = math.floatMantissaBits(T);
+ const frac_bits: comptime_int = math.floatFractionalBits(T);
+ const exp_min: comptime_int = math.floatExponentMin(T);
+
+ const ExpInt: type = std.meta.Int(.unsigned, exp_bits);
+ const MantInt: type = std.meta.Int(.unsigned, mant_bits);
+ const FracInt: type = std.meta.Int(.unsigned, frac_bits);
+
+ const unreal_exponent: comptime_int = (1 << exp_bits) - 1;
+ const bias: comptime_int = (1 << (exp_bits - 1)) - 2;
+ const exp_mask: comptime_int = unreal_exponent << mant_bits;
+ const zero_exponent: comptime_int = bias << mant_bits;
+ const sign_mask: comptime_int = 1 << (bits - 1);
+ const not_exp: comptime_int = ~@as(Int, exp_mask);
+ const ones_place: comptime_int = mant_bits - frac_bits;
+ const extra_denorm_shift: comptime_int = 1 - ones_place;
+
+ var result: Frexp(T) = undefined;
+ var v: Int = @bitCast(x);
+
+ const m: MantInt = @truncate(v);
+ const e: ExpInt = @truncate(v >> mant_bits);
+
+ switch (e) {
+ 0 => {
+ if (m != 0) {
+ // subnormal
+ const offset = @clz(m);
+ const shift = offset + extra_denorm_shift;
+
+ v &= sign_mask;
+ v |= zero_exponent;
+ v |= math.shl(MantInt, m, shift);
+
+ result.exponent = exp_min - @as(i32, offset) + ones_place;
+ } else {
+ // +-0 = (+-0, 0)
+ result.exponent = 0;
+ }
+ },
+ unreal_exponent => {
+ // +-nan -> {+-nan, undefined}
+ result.exponent = undefined;
+
+ // +-inf -> {+-inf, 0}
+ if (@as(FracInt, @truncate(v)) == 0)
+ result.exponent = 0;
+ },
+ else => {
+ // normal
+ v &= not_exp;
+ v |= zero_exponent;
+ result.exponent = @as(i32, e) - bias;
+ },
}
- result.exponent = e - 0x7E;
- y &= 0x807FFFFF;
- y |= 0x3F000000;
- result.significand = @as(f32, @bitCast(y));
+ result.significand = @bitCast(v);
return result;
}
-fn frexp64(x: f64) Frexp(f64) {
- var result: Frexp(f64) = undefined;
-
- var y = @as(u64, @bitCast(x));
- const e = @as(i32, @intCast(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;
+/// Generate a namespace of tests for frexp on values of the given type
+fn FrexpTests(comptime Float: type) type {
+ return struct {
+ const T = Float;
+ test "normal" {
+ const epsilon = 1e-6;
+ var r: Frexp(T) = undefined;
+
+ r = frexp(@as(T, 1.3));
+ try expectApproxEqAbs(0.65, r.significand, epsilon);
+ try expectEqual(1, r.exponent);
+
+ r = frexp(@as(T, 78.0234));
+ try expectApproxEqAbs(0.609558, r.significand, epsilon);
+ try expectEqual(7, r.exponent);
+
+ r = frexp(@as(T, -1234.5678));
+ try expectEqual(11, r.exponent);
+ try expectApproxEqAbs(-0.602816, r.significand, epsilon);
}
- 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;
+ test "max" {
+ const exponent = math.floatExponentMax(T) + 1;
+ const significand = 1.0 - math.floatEps(T) / 2;
+ const r: Frexp(T) = frexp(math.floatMax(T));
+ try expectEqual(exponent, r.exponent);
+ try expectEqual(significand, r.significand);
}
-
- return result;
- }
-
- result.exponent = e - 0x3FE;
- y &= 0x800FFFFFFFFFFFFF;
- y |= 0x3FE0000000000000;
- result.significand = @as(f64, @bitCast(y));
- return result;
-}
-
-fn frexp128(x: f128) Frexp(f128) {
- var result: Frexp(f128) = undefined;
-
- var y = @as(u128, @bitCast(x));
- const e = @as(i32, @intCast(y >> 112)) & 0x7FFF;
-
- if (e == 0) {
- if (x != 0) {
- // subnormal
- result = frexp128(x * 0x1.0p120);
- result.exponent -= 120;
- } else {
- // frexp(+-0) = (+-0, 0)
- result.significand = x;
- result.exponent = 0;
+ test "min" {
+ const exponent = math.floatExponentMin(T) + 1;
+ const r: Frexp(T) = frexp(math.floatMin(T));
+ try expectEqual(exponent, r.exponent);
+ try expectEqual(0.5, r.significand);
}
- return result;
- } else if (e == 0x7FFF) {
- // frexp(nan) = (nan, undefined)
- result.significand = x;
- result.exponent = undefined;
-
- // frexp(+-inf) = (+-inf, 0)
- if (math.isInf(x)) {
- result.exponent = 0;
+ test "subnormal" {
+ const normal_min_exponent = math.floatExponentMin(T) + 1;
+ const exponent = normal_min_exponent - math.floatFractionalBits(T);
+ const r: Frexp(T) = frexp(math.floatTrueMin(T));
+ try expectEqual(exponent, r.exponent);
+ try expectEqual(0.5, r.significand);
}
+ test "zero" {
+ var r: Frexp(T) = undefined;
- return result;
- }
-
- result.exponent = e - 0x3FFE;
- y &= 0x8000FFFFFFFFFFFFFFFFFFFFFFFFFFFF;
- y |= 0x3FFE0000000000000000000000000000;
- result.significand = @as(f128, @bitCast(y));
- return result;
-}
-
-test "type dispatch" {
- const a = frexp(@as(f32, 1.3));
- const b = frexp32(1.3);
- try expect(a.significand == b.significand and a.exponent == b.exponent);
-
- const c = frexp(@as(f64, 1.3));
- const d = frexp64(1.3);
- try expect(c.significand == d.significand and c.exponent == d.exponent);
-
- const e = frexp(@as(f128, 1.3));
- const f = frexp128(1.3);
- try expect(e.significand == f.significand and e.exponent == f.exponent);
-}
-
-test "32" {
- const epsilon = 0.000001;
- var r: Frexp(f32) = undefined;
-
- r = frexp32(1.3);
- try expect(math.approxEqAbs(f32, r.significand, 0.65, epsilon) and r.exponent == 1);
-
- r = frexp32(78.0234);
- try expect(math.approxEqAbs(f32, r.significand, 0.609558, epsilon) and r.exponent == 7);
-}
-
-test "64" {
- const epsilon = 0.000001;
- var r: Frexp(f64) = undefined;
-
- r = frexp64(1.3);
- try expect(math.approxEqAbs(f64, r.significand, 0.65, epsilon) and r.exponent == 1);
-
- r = frexp64(78.0234);
- try expect(math.approxEqAbs(f64, r.significand, 0.609558, epsilon) and r.exponent == 7);
-}
-
-test "128" {
- const epsilon = 0.000001;
- var r: Frexp(f128) = undefined;
-
- r = frexp128(1.3);
- try expect(math.approxEqAbs(f128, r.significand, 0.65, epsilon) and r.exponent == 1);
+ r = frexp(@as(T, 0.0));
+ try expectEqual(0, r.exponent);
+ try expect(math.isPositiveZero(r.significand));
- r = frexp128(78.0234);
- try expect(math.approxEqAbs(f128, r.significand, 0.609558, epsilon) and r.exponent == 7);
-}
-
-test "32 special" {
- var r: Frexp(f32) = undefined;
-
- r = frexp32(0.0);
- try expect(r.significand == 0.0 and r.exponent == 0);
-
- r = frexp32(-0.0);
- try expect(r.significand == -0.0 and r.exponent == 0);
-
- r = frexp32(math.inf(f32));
- try expect(math.isPositiveInf(r.significand) and r.exponent == 0);
+ r = frexp(@as(T, -0.0));
+ try expectEqual(0, r.exponent);
+ try expect(math.isNegativeZero(r.significand));
+ }
+ test "inf" {
+ var r: Frexp(T) = undefined;
- r = frexp32(-math.inf(f32));
- try expect(math.isNegativeInf(r.significand) and r.exponent == 0);
+ r = frexp(math.inf(T));
+ try expectEqual(0, r.exponent);
+ try expect(math.isPositiveInf(r.significand));
- r = frexp32(math.nan(f32));
- try expect(math.isNan(r.significand));
+ r = frexp(-math.inf(T));
+ try expectEqual(0, r.exponent);
+ try expect(math.isNegativeInf(r.significand));
+ }
+ test "nan" {
+ const r: Frexp(T) = frexp(math.nan(T));
+ try expect(math.isNan(r.significand));
+ }
+ };
}
-test "64 special" {
- var r: Frexp(f64) = undefined;
-
- r = frexp64(0.0);
- try expect(r.significand == 0.0 and r.exponent == 0);
-
- r = frexp64(-0.0);
- try expect(r.significand == -0.0 and r.exponent == 0);
-
- r = frexp64(math.inf(f64));
- try expect(math.isPositiveInf(r.significand) and r.exponent == 0);
-
- r = frexp64(-math.inf(f64));
- try expect(math.isNegativeInf(r.significand) and r.exponent == 0);
-
- r = frexp64(math.nan(f64));
- try expect(math.isNan(r.significand));
+// Generate tests for each floating point type
+comptime {
+ for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
+ _ = FrexpTests(T);
+ }
}
-test "128 special" {
- var r: Frexp(f128) = undefined;
-
- r = frexp128(0.0);
- try expect(r.significand == 0.0 and r.exponent == 0);
-
- r = frexp128(-0.0);
- try expect(r.significand == -0.0 and r.exponent == 0);
-
- r = frexp128(math.inf(f128));
- try expect(math.isPositiveInf(r.significand) and r.exponent == 0);
-
- r = frexp128(-math.inf(f128));
- try expect(math.isNegativeInf(r.significand) and r.exponent == 0);
-
- r = frexp128(math.nan(f128));
- try expect(math.isNan(r.significand));
+test frexp {
+ inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
+ const max_exponent = math.floatExponentMax(T) + 1;
+ const min_exponent = math.floatExponentMin(T) + 1;
+ const truemin_exponent = min_exponent - math.floatFractionalBits(T);
+
+ var result: Frexp(T) = undefined;
+ comptime var x: T = undefined;
+
+ // basic usage
+ // value -> {significand, exponent},
+ // value == significand * (2 ^ exponent)
+ x = 1234.5678;
+ result = frexp(x);
+ try expectEqual(11, result.exponent);
+ try expectApproxEqAbs(0.602816, result.significand, 1e-6);
+ try expectEqual(x, math.ldexp(result.significand, result.exponent));
+
+ // float maximum
+ x = math.floatMax(T);
+ result = frexp(x);
+ try expectEqual(max_exponent, result.exponent);
+ try expectEqual(1.0 - math.floatEps(T) / 2, result.significand);
+ try expectEqual(x, math.ldexp(result.significand, result.exponent));
+
+ // float minimum
+ x = math.floatMin(T);
+ result = frexp(x);
+ try expectEqual(min_exponent, result.exponent);
+ try expectEqual(0.5, result.significand);
+ try expectEqual(x, math.ldexp(result.significand, result.exponent));
+
+ // float true minimum
+ // subnormal -> {normal, exponent}
+ x = math.floatTrueMin(T);
+ result = frexp(x);
+ try expectEqual(truemin_exponent, result.exponent);
+ try expectEqual(0.5, result.significand);
+ try expectEqual(x, math.ldexp(result.significand, result.exponent));
+
+ // infinity -> {infinity, zero} (+)
+ result = frexp(math.inf(T));
+ try expectEqual(0, result.exponent);
+ try expect(math.isPositiveInf(result.significand));
+
+ // infinity -> {infinity, zero} (-)
+ result = frexp(-math.inf(T));
+ try expectEqual(0, result.exponent);
+ try expect(math.isNegativeInf(result.significand));
+
+ // zero -> {zero, zero} (+)
+ result = frexp(@as(T, 0.0));
+ try expectEqual(0, result.exponent);
+ try expect(math.isPositiveZero(result.significand));
+
+ // zero -> {zero, zero} (-)
+ result = frexp(@as(T, -0.0));
+ try expectEqual(0, result.exponent);
+ try expect(math.isNegativeZero(result.significand));
+
+ // nan -> {nan, undefined}
+ result = frexp(math.nan(T));
+ try expect(math.isNan(result.significand));
+ }
}