diff --git a/lib/std/math.zig b/lib/std/math.zig index 0f3a22d2b1..4d2c1395f6 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -241,6 +241,7 @@ pub const isNegativeInf = @import("math/isinf.zig").isNegativeInf; pub const isNormal = @import("math/isnormal.zig").isNormal; pub const signbit = @import("math/signbit.zig").signbit; pub const scalbn = @import("math/scalbn.zig").scalbn; +pub const ldexp = @import("math/ldexp.zig").ldexp; pub const pow = @import("math/pow.zig").pow; pub const powi = @import("math/powi.zig").powi; pub const sqrt = @import("math/sqrt.zig").sqrt; diff --git a/lib/std/math/ldexp.zig b/lib/std/math/ldexp.zig new file mode 100644 index 0000000000..228bd7dd39 --- /dev/null +++ b/lib/std/math/ldexp.zig @@ -0,0 +1,92 @@ +// 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/ldexpf.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/ldexp.c + +const std = @import("std"); +const math = std.math; +const assert = std.debug.assert; +const expect = std.testing.expect; + +/// Returns x * 2^n. +pub fn ldexp(x: anytype, n: i32) @TypeOf(x) { + var base = x; + var shift = n; + + const T = @TypeOf(base); + const IntT = std.meta.Int(.unsigned, @bitSizeOf(T)); + if (@typeInfo(T) != .Float) { + @compileError("ldexp not implemented for " ++ @typeName(T)); + } + + const mantissa_bits = math.floatMantissaBits(T); + const exponent_bits = math.floatExponentBits(T); + const exponent_bias = (1 << (exponent_bits - 1)) - 1; + const exponent_min = 1 - exponent_bias; + const exponent_max = exponent_bias; + + // fix double rounding errors in subnormal ranges + // https://git.musl-libc.org/cgit/musl/commit/src/math/ldexp.c?id=8c44a060243f04283ca68dad199aab90336141db + const scale_min_expo = exponent_min + mantissa_bits + 1; + const scale_min = @bitCast(T, @as(IntT, scale_min_expo + exponent_bias) << mantissa_bits); + const scale_max = @bitCast(T, @intCast(IntT, exponent_max + exponent_bias) << mantissa_bits); + + // scale `shift` within floating point limits, if possible + // second pass is possible due to subnormal range + // third pass always results in +/-0.0 or +/-inf + if (shift > exponent_max) { + base *= scale_max; + shift -= exponent_max; + if (shift > exponent_max) { + base *= scale_max; + shift -= exponent_max; + if (shift > exponent_max) shift = exponent_max; + } + } else if (shift < exponent_min) { + base *= scale_min; + shift -= scale_min_expo; + if (shift < exponent_min) { + base *= scale_min; + shift -= scale_min_expo; + if (shift < exponent_min) shift = exponent_min; + } + } + + return base * @bitCast(T, @intCast(IntT, shift + exponent_bias) << mantissa_bits); +} + +test "math.ldexp" { + // basic usage + try expect(ldexp(@as(f16, 1.5), 4) == 24.0); + try expect(ldexp(@as(f32, 1.5), 4) == 24.0); + try expect(ldexp(@as(f64, 1.5), 4) == 24.0); + try expect(ldexp(@as(f128, 1.5), 4) == 24.0); + + // subnormals + try expect(math.isNormal(ldexp(@as(f16, 1.0), -14))); + try expect(!math.isNormal(ldexp(@as(f16, 1.0), -15))); + try expect(math.isNormal(ldexp(@as(f32, 1.0), -126))); + try expect(!math.isNormal(ldexp(@as(f32, 1.0), -127))); + try expect(math.isNormal(ldexp(@as(f64, 1.0), -1022))); + try expect(!math.isNormal(ldexp(@as(f64, 1.0), -1023))); + try expect(math.isNormal(ldexp(@as(f128, 1.0), -16382))); + try expect(!math.isNormal(ldexp(@as(f128, 1.0), -16383))); + // unreliable due to lack of native f16 support, see talk on PR #8733 + // try expect(ldexp(@as(f16, 0x1.1FFp-1), -14 - 9) == math.f16_true_min); + try expect(ldexp(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.f32_true_min); + try expect(ldexp(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.f64_true_min); + try expect(ldexp(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.f128_true_min); + + // float limits + try expect(ldexp(@as(f32, math.f32_max), -128 - 149) > 0.0); + try expect(ldexp(@as(f32, math.f32_max), -128 - 149 - 1) == 0.0); + try expect(!math.isPositiveInf(ldexp(@as(f16, math.f16_true_min), 15 + 24))); + try expect(math.isPositiveInf(ldexp(@as(f16, math.f16_true_min), 15 + 24 + 1))); + try expect(!math.isPositiveInf(ldexp(@as(f32, math.f32_true_min), 127 + 149))); + try expect(math.isPositiveInf(ldexp(@as(f32, math.f32_true_min), 127 + 149 + 1))); + try expect(!math.isPositiveInf(ldexp(@as(f64, math.f64_true_min), 1023 + 1074))); + try expect(math.isPositiveInf(ldexp(@as(f64, math.f64_true_min), 1023 + 1074 + 1))); + try expect(!math.isPositiveInf(ldexp(@as(f128, math.f128_true_min), 16383 + 16494))); + try expect(math.isPositiveInf(ldexp(@as(f128, math.f128_true_min), 16383 + 16494 + 1))); +} diff --git a/lib/std/math/scalbn.zig b/lib/std/math/scalbn.zig index 30494c820a..294ee4abf0 100644 --- a/lib/std/math/scalbn.zig +++ b/lib/std/math/scalbn.zig @@ -1,92 +1,15 @@ -// 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/scalbnf.c -// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c - const std = @import("std"); -const math = std.math; -const assert = std.debug.assert; const expect = std.testing.expect; -/// Returns x * 2^n. -pub fn scalbn(x: anytype, n: i32) @TypeOf(x) { - var base = x; - var shift = n; - - const T = @TypeOf(base); - const IntT = std.meta.Int(.unsigned, @bitSizeOf(T)); - if (@typeInfo(T) != .Float) { - @compileError("scalbn not implemented for " ++ @typeName(T)); - } - - const mantissa_bits = math.floatMantissaBits(T); - const exponent_bits = math.floatExponentBits(T); - const exponent_bias = (1 << (exponent_bits - 1)) - 1; - const exponent_min = 1 - exponent_bias; - const exponent_max = exponent_bias; - - // fix double rounding errors in subnormal ranges - // https://git.musl-libc.org/cgit/musl/commit/src/math/scalbn.c?id=8c44a060243f04283ca68dad199aab90336141db - const scale_min_expo = exponent_min + mantissa_bits + 1; - const scale_min = @bitCast(T, @as(IntT, scale_min_expo + exponent_bias) << mantissa_bits); - const scale_max = @bitCast(T, @intCast(IntT, exponent_max + exponent_bias) << mantissa_bits); - - // scale `shift` within floating point limits, if possible - // second pass is possible due to subnormal range - // third pass always results in +/-0.0 or +/-inf - if (shift > exponent_max) { - base *= scale_max; - shift -= exponent_max; - if (shift > exponent_max) { - base *= scale_max; - shift -= exponent_max; - if (shift > exponent_max) shift = exponent_max; - } - } else if (shift < exponent_min) { - base *= scale_min; - shift -= scale_min_expo; - if (shift < exponent_min) { - base *= scale_min; - shift -= scale_min_expo; - if (shift < exponent_min) shift = exponent_min; - } - } - - return base * @bitCast(T, @intCast(IntT, shift + exponent_bias) << mantissa_bits); -} +/// Returns a * FLT_RADIX ^ exp. +/// +/// Zig only supports binary radix IEEE-754 floats. Hence FLT_RADIX=2, and this is an alias for ldexp. +pub const scalbn = @import("ldexp.zig").ldexp; test "math.scalbn" { - // basic usage + // Verify we are using radix 2. try expect(scalbn(@as(f16, 1.5), 4) == 24.0); try expect(scalbn(@as(f32, 1.5), 4) == 24.0); try expect(scalbn(@as(f64, 1.5), 4) == 24.0); try expect(scalbn(@as(f128, 1.5), 4) == 24.0); - - // subnormals - try expect(math.isNormal(scalbn(@as(f16, 1.0), -14))); - try expect(!math.isNormal(scalbn(@as(f16, 1.0), -15))); - try expect(math.isNormal(scalbn(@as(f32, 1.0), -126))); - try expect(!math.isNormal(scalbn(@as(f32, 1.0), -127))); - try expect(math.isNormal(scalbn(@as(f64, 1.0), -1022))); - try expect(!math.isNormal(scalbn(@as(f64, 1.0), -1023))); - try expect(math.isNormal(scalbn(@as(f128, 1.0), -16382))); - try expect(!math.isNormal(scalbn(@as(f128, 1.0), -16383))); - // unreliable due to lack of native f16 support, see talk on PR #8733 - // try expect(scalbn(@as(f16, 0x1.1FFp-1), -14 - 9) == math.f16_true_min); - try expect(scalbn(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.f32_true_min); - try expect(scalbn(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.f64_true_min); - try expect(scalbn(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.f128_true_min); - - // float limits - try expect(scalbn(@as(f32, math.f32_max), -128 - 149) > 0.0); - try expect(scalbn(@as(f32, math.f32_max), -128 - 149 - 1) == 0.0); - try expect(!math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24))); - try expect(math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24 + 1))); - try expect(!math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149))); - try expect(math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149 + 1))); - try expect(!math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074))); - try expect(math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074 + 1))); - try expect(!math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494))); - try expect(math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494 + 1))); }