srctree

Sean parent 62929d59 f32723a2
Update frexp.zig (#19370)

1. Entirely rewrote frexp with generics, reducing the implementation to a single function and enabling parameters of types f80 and f16 2. Expanded upon the tests, making them more descriptive and comprehensive, and automatically generating the test bodies for each floating point type 3. Added a doctest for frexp

inline split

lib/std/math/frexp.zig added: 198, removed: 221, total 0

@@ -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)),

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.significand = @bitCast(v);

return result;

}

/// 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);

}

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);

}

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);

}

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;

r = frexp(@as(T, 0.0));

try expectEqual(0, r.exponent);

try expect(math.isPositiveZero(r.significand));

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 = frexp(math.inf(T));

try expectEqual(0, r.exponent);

try expect(math.isPositiveInf(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));

}

};

}

// 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;

// Generate tests for each floating point type

comptime {

for ([_]type{ f16, f32, f64, f80, f128 }) |T| {

_ = FrexpTests(T);

}

result.exponent = e - 0x7E;

y &= 0x807FFFFF;

y |= 0x3F000000;

result.significand = @as(f32, @bitCast(y));

return result;

}

fn frexp64(x: f64) Frexp(f64) {

var result: Frexp(f64) = undefined;

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 y = @as(u64, @bitCast(x));

const e = @as(i32, @intCast(y >> 52)) & 0x7FF;

var result: Frexp(T) = undefined;

comptime var x: T = undefined;

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;

// 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));

// frexp(+-inf) = (+-inf, 0)

if (math.isInf(x)) {

result.exponent = 0;

}

// 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));

return result;

// 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));

}

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;

}

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;

}

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 = 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 = frexp32(-math.inf(f32));

try expect(math.isNegativeInf(r.significand) and r.exponent == 0);

r = frexp32(math.nan(f32));

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));

}

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));

}