diff options
| author | Andrew Kelley <andrew@ziglang.org> | 2022-05-06 19:22:40 -0700 |
|---|---|---|
| committer | Andrew Kelley <andrew@ziglang.org> | 2022-05-06 22:41:00 -0700 |
| commit | ec95e00e28cb23f37dc097f71afd7090e947a1cd (patch) | |
| tree | a7393f13c3d2c7895eb3687d0ebd7f3205699289 /lib/compiler_rt/log2.zig | |
| parent | 3b60ab4872355f0b9a9c7d0794ca8b548ab99412 (diff) | |
| download | zig-ec95e00e28cb23f37dc097f71afd7090e947a1cd.tar.gz zig-ec95e00e28cb23f37dc097f71afd7090e947a1cd.zip | |
flatten lib/std/special and improve "pkg inside another" logic
stage2: change logic for detecting whether the main package is inside
the std package. Previously it relied on realpath() which is not portable.
This uses resolve() which is how imports already work.
* stage2: fix cleanup bug when creating Module
* flatten lib/std/special/* to lib/*
- this was motivated by making main_pkg_is_inside_std false for
compiler_rt & friends.
* rename "mini libc" to "universal libc"
Diffstat (limited to 'lib/compiler_rt/log2.zig')
| -rw-r--r-- | lib/compiler_rt/log2.zig | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/lib/compiler_rt/log2.zig b/lib/compiler_rt/log2.zig new file mode 100644 index 0000000000..2c2d620c3d --- /dev/null +++ b/lib/compiler_rt/log2.zig @@ -0,0 +1,185 @@ +// 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/log2f.c +// https://git.musl-libc.org/cgit/musl/tree/src/math/log2.c + +const std = @import("std"); +const math = std.math; +const expect = std.testing.expect; +const maxInt = std.math.maxInt; + +pub fn __log2h(a: f16) callconv(.C) f16 { + // TODO: more efficient implementation + return @floatCast(f16, log2f(a)); +} + +pub fn log2f(x_: f32) callconv(.C) f32 { + const ivln2hi: f32 = 1.4428710938e+00; + const ivln2lo: f32 = -1.7605285393e-04; + const Lg1: f32 = 0xaaaaaa.0p-24; + const Lg2: f32 = 0xccce13.0p-25; + const Lg3: f32 = 0x91e9ee.0p-25; + const Lg4: f32 = 0xf89e26.0p-26; + + var x = x_; + var u = @bitCast(u32, x); + var ix = u; + var k: i32 = 0; + + // x < 2^(-126) + if (ix < 0x00800000 or ix >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f32); + } + // log(-#) = nan + if (ix >> 31 != 0) { + return math.nan(f32); + } + + k -= 25; + x *= 0x1.0p25; + ix = @bitCast(u32, x); + } else if (ix >= 0x7F800000) { + return x; + } else if (ix == 0x3F800000) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + ix += 0x3F800000 - 0x3F3504F3; + k += @intCast(i32, ix >> 23) - 0x7F; + ix = (ix & 0x007FFFFF) + 0x3F3504F3; + x = @bitCast(f32, ix); + + const f = x - 1.0; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * Lg4); + const t2 = z * (Lg1 + w * Lg3); + const R = t2 + t1; + const hfsq = 0.5 * f * f; + + var hi = f - hfsq; + u = @bitCast(u32, hi); + u &= 0xFFFFF000; + hi = @bitCast(f32, u); + const lo = f - hi - hfsq + s * (hfsq + R); + return (lo + hi) * ivln2lo + lo * ivln2hi + hi * ivln2hi + @intToFloat(f32, k); +} + +pub fn log2(x_: f64) callconv(.C) f64 { + const ivln2hi: f64 = 1.44269504072144627571e+00; + const ivln2lo: f64 = 1.67517131648865118353e-10; + const Lg1: f64 = 6.666666666666735130e-01; + const Lg2: f64 = 3.999999999940941908e-01; + const Lg3: f64 = 2.857142874366239149e-01; + const Lg4: f64 = 2.222219843214978396e-01; + const Lg5: f64 = 1.818357216161805012e-01; + const Lg6: f64 = 1.531383769920937332e-01; + const Lg7: f64 = 1.479819860511658591e-01; + + var x = x_; + var ix = @bitCast(u64, x); + var hx = @intCast(u32, ix >> 32); + var k: i32 = 0; + + if (hx < 0x00100000 or hx >> 31 != 0) { + // log(+-0) = -inf + if (ix << 1 == 0) { + return -math.inf(f64); + } + // log(-#) = nan + if (hx >> 31 != 0) { + return math.nan(f64); + } + + // subnormal, scale x + k -= 54; + x *= 0x1.0p54; + hx = @intCast(u32, @bitCast(u64, x) >> 32); + } else if (hx >= 0x7FF00000) { + return x; + } else if (hx == 0x3FF00000 and ix << 32 == 0) { + return 0; + } + + // x into [sqrt(2) / 2, sqrt(2)] + hx += 0x3FF00000 - 0x3FE6A09E; + k += @intCast(i32, hx >> 20) - 0x3FF; + hx = (hx & 0x000FFFFF) + 0x3FE6A09E; + ix = (@as(u64, hx) << 32) | (ix & 0xFFFFFFFF); + x = @bitCast(f64, ix); + + const f = x - 1.0; + const hfsq = 0.5 * f * f; + const s = f / (2.0 + f); + const z = s * s; + const w = z * z; + const t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); + const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); + const R = t2 + t1; + + // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f) + var hi = f - hfsq; + var hii = @bitCast(u64, hi); + hii &= @as(u64, maxInt(u64)) << 32; + hi = @bitCast(f64, hii); + const lo = f - hi - hfsq + s * (hfsq + R); + + var val_hi = hi * ivln2hi; + var val_lo = (lo + hi) * ivln2lo + lo * ivln2hi; + + // spadd(val_hi, val_lo, y) + const y = @intToFloat(f64, k); + const ww = y + val_hi; + val_lo += (y - ww) + val_hi; + val_hi = ww; + + return val_lo + val_hi; +} + +pub fn __log2x(a: f80) callconv(.C) f80 { + // TODO: more efficient implementation + return @floatCast(f80, log2q(a)); +} + +pub fn log2q(a: f128) callconv(.C) f128 { + return math.log2(a); +} + +test "log2_32" { + const epsilon = 0.000001; + + try expect(math.approxEqAbs(f32, log2f(0.2), -2.321928, epsilon)); + try expect(math.approxEqAbs(f32, log2f(0.8923), -0.164399, epsilon)); + try expect(math.approxEqAbs(f32, log2f(1.5), 0.584962, epsilon)); + try expect(math.approxEqAbs(f32, log2f(37.45), 5.226894, epsilon)); + try expect(math.approxEqAbs(f32, log2f(123123.234375), 16.909744, epsilon)); +} + +test "log2_64" { + const epsilon = 0.000001; + + try expect(math.approxEqAbs(f64, log2(0.2), -2.321928, epsilon)); + try expect(math.approxEqAbs(f64, log2(0.8923), -0.164399, epsilon)); + try expect(math.approxEqAbs(f64, log2(1.5), 0.584962, epsilon)); + try expect(math.approxEqAbs(f64, log2(37.45), 5.226894, epsilon)); + try expect(math.approxEqAbs(f64, log2(123123.234375), 16.909744, epsilon)); +} + +test "log2_32.special" { + try expect(math.isPositiveInf(log2f(math.inf(f32)))); + try expect(math.isNegativeInf(log2f(0.0))); + try expect(math.isNan(log2f(-1.0))); + try expect(math.isNan(log2f(math.nan(f32)))); +} + +test "log2_64.special" { + try expect(math.isPositiveInf(log2(math.inf(f64)))); + try expect(math.isNegativeInf(log2(0.0))); + try expect(math.isNan(log2(-1.0))); + try expect(math.isNan(log2(math.nan(f64)))); +} |
