diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/std/math/frexp.zig | 405 |
1 files changed, 191 insertions, 214 deletions
diff --git a/lib/std/math/frexp.zig b/lib/std/math/frexp.zig index f295b959cb..97b578db98 100644 --- a/lib/std/math/frexp.zig +++ b/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)); + } } |
