From daafc14f30605a795dd7f47dd10f785f3a4ab187 Mon Sep 17 00:00:00 2001 From: Federico Stra Date: Wed, 13 Sep 2023 13:04:10 +0200 Subject: [PATCH 1/5] std.math.log_int: implement integer logarithm without using float math --- lib/std/math.zig | 1 + lib/std/math/log.zig | 7 +++- lib/std/math/log_int.zig | 79 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 85 insertions(+), 2 deletions(-) create mode 100644 lib/std/math/log_int.zig diff --git a/lib/std/math.zig b/lib/std/math.zig index 558629b74368..16676283e119 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -241,6 +241,7 @@ pub const log = @import("math/log.zig").log; pub const log2 = @import("math/log2.zig").log2; pub const log10 = @import("math/log10.zig").log10; pub const log10_int = @import("math/log10.zig").log10_int; +pub const log_int = @import("math/log_int.zig").log_int; pub const log1p = @import("math/log1p.zig").log1p; pub const asinh = @import("math/asinh.zig").asinh; pub const acosh = @import("math/acosh.zig").acosh; diff --git a/lib/std/math/log.zig b/lib/std/math/log.zig index 9f27130ce1dc..4dcfd35d1277 100644 --- a/lib/std/math/log.zig +++ b/lib/std/math/log.zig @@ -23,14 +23,17 @@ pub fn log(comptime T: type, base: T, x: T) T { .ComptimeFloat => { return @as(comptime_float, @log(@as(f64, x)) / @log(float_base)); }, + + // TODO: implement integer log without using float math. + // The present implementation is incorrect, for example + // `log(comptime_int, 9, 59049)` should return `5` and not `4`. .ComptimeInt => { return @as(comptime_int, @floor(@log(@as(f64, x)) / @log(float_base))); }, - // TODO implement integer log without using float math .Int => |IntType| switch (IntType.signedness) { .signed => @compileError("log not implemented for signed integers"), - .unsigned => return @as(T, @intFromFloat(@floor(@log(@as(f64, @floatFromInt(x))) / @log(float_base)))), + .unsigned => return @as(T, math.log_int(T, base, x)), }, .Float => { diff --git a/lib/std/math/log_int.zig b/lib/std/math/log_int.zig new file mode 100644 index 000000000000..128edf5d3eb9 --- /dev/null +++ b/lib/std/math/log_int.zig @@ -0,0 +1,79 @@ +const std = @import("../std.zig"); +const math = std.math; +const testing = std.testing; +const assert = std.debug.assert; +const Log2Int = math.Log2Int; + +/// Returns the logarithm of `x` for the provided `base`, rounding down to the nearest integer. +/// Asserts that `base > 1` and `x != 0`. +pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) { + if (@typeInfo(T) != .Int or @typeInfo(T).Int.signedness != .unsigned) + @compileError("log_int requires an unsigned integer, found " ++ @typeName(T)); + + assert(base > 1 and x != 0); + + var exponent: Log2Int(T) = 0; + var power: T = 1; + while (power <= x / base) { + power *= base; + exponent += 1; + } + + return exponent; +} + +test "log" { + // test all unsigned integers with 2, 3, ..., 64 bits + inline for (2..64 + 1) |bits| { + const T = @Type(std.builtin.Type{ + .Int = std.builtin.Type.Int{ .signedness = .unsigned, .bits = @intCast(bits) }, + }); + + // for base = 2, 3, ..., min(maxInt(T),1024) + var base: T = 1; + while (base < math.maxInt(T) and base <= 1024) { + base += 1; + + // test `log(1) == 0` + try testing.expectEqual(@as(Log2Int(T), 0), log_int(T, base, 1)); + + // for powers `pow = base^exp` that fit inside T + var exp: Log2Int(T) = 0; + var pow: T = 1; + while (pow < math.maxInt(T) / base) { + exp += 1; + pow *= base; + + // test that `log_int` correctly detects the threshold + try testing.expectEqual(exp - 1, log_int(T, base, pow - 1)); + try testing.expectEqual(exp, log_int(T, base, pow)); + } + } + } +} + +test "log2" { + const types = [_]type{ u2, u3, u4, u8, u16, u24 }; + inline for (types) |T| { + var n: T = 0; + while (n < math.maxInt(T)) { + n += 1; + const special = math.log2_int(T, n); + const general = log_int(T, 2, n); + try testing.expectEqual(special, general); + } + } +} + +test "log10" { + const types = [_]type{ u4, u5, u6, u8, u16, u24 }; + inline for (types) |T| { + var n: T = 0; + while (n < math.maxInt(T)) { + n += 1; + const special = math.log10_int(n); + const general = log_int(T, 10, n); + try testing.expectEqual(special, general); + } + } +} From c305efdd8419d5ebb3d5b1109c52aa3c7fd2de32 Mon Sep 17 00:00:00 2001 From: Federico Stra Date: Wed, 13 Sep 2023 14:37:06 +0200 Subject: [PATCH 2/5] std.math: ensure that log_int is tested --- lib/std/math.zig | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/std/math.zig b/lib/std/math.zig index 16676283e119..dccb127c9509 100644 --- a/lib/std/math.zig +++ b/lib/std/math.zig @@ -363,6 +363,7 @@ test { _ = log2; _ = log10; _ = log10_int; + _ = log_int; _ = log1p; _ = asinh; _ = acosh; From c64b9bd800688580f009423929c0b8b05b27c74b Mon Sep 17 00:00:00 2001 From: Federico Stra Date: Wed, 13 Sep 2023 15:05:32 +0200 Subject: [PATCH 3/5] std.math.log_int: improve tests --- lib/std/math/log_int.zig | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/lib/std/math/log_int.zig b/lib/std/math/log_int.zig index 128edf5d3eb9..3114d775c5e8 100644 --- a/lib/std/math/log_int.zig +++ b/lib/std/math/log_int.zig @@ -22,8 +22,9 @@ pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) { return exponent; } -test "log" { - // test all unsigned integers with 2, 3, ..., 64 bits +test "math.log_int" { + // Test all unsigned integers with 2, 3, ..., 64 bits. + // We cannot test 0 or 1 bits since base must be > 1. inline for (2..64 + 1) |bits| { const T = @Type(std.builtin.Type{ .Int = std.builtin.Type.Int{ .signedness = .unsigned, .bits = @intCast(bits) }, @@ -34,17 +35,18 @@ test "log" { while (base < math.maxInt(T) and base <= 1024) { base += 1; - // test `log(1) == 0` + // test that `log_int(T, base, 1) == 0` try testing.expectEqual(@as(Log2Int(T), 0), log_int(T, base, 1)); - // for powers `pow = base^exp` that fit inside T + // For powers `pow = base^exp > 1` that fit inside T, + // test that `log_int` correctly detects the jump in the logarithm + // from `log(pow-1) == exp-1` to `log(pow) == exp`. var exp: Log2Int(T) = 0; var pow: T = 1; - while (pow < math.maxInt(T) / base) { + while (pow <= math.maxInt(T) / base) { exp += 1; pow *= base; - // test that `log_int` correctly detects the threshold try testing.expectEqual(exp - 1, log_int(T, base, pow - 1)); try testing.expectEqual(exp, log_int(T, base, pow)); } @@ -52,8 +54,8 @@ test "log" { } } -test "log2" { - const types = [_]type{ u2, u3, u4, u8, u16, u24 }; +test "math.log_int vs math.log2" { + const types = [_]type{ u2, u3, u4, u8, u16 }; inline for (types) |T| { var n: T = 0; while (n < math.maxInt(T)) { @@ -65,8 +67,8 @@ test "log2" { } } -test "log10" { - const types = [_]type{ u4, u5, u6, u8, u16, u24 }; +test "math.log_int vs math.log10" { + const types = [_]type{ u4, u5, u6, u8, u16 }; inline for (types) |T| { var n: T = 0; while (n < math.maxInt(T)) { From 8a298473cfa2fd023bfb138721784bf83af97a68 Mon Sep 17 00:00:00 2001 From: Federico Stra Date: Thu, 14 Sep 2023 15:06:05 +0200 Subject: [PATCH 4/5] std.math.log_int: modify the assertion --- lib/std/math/log_int.zig | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/std/math/log_int.zig b/lib/std/math/log_int.zig index 3114d775c5e8..bcdd73cb30bf 100644 --- a/lib/std/math/log_int.zig +++ b/lib/std/math/log_int.zig @@ -5,12 +5,12 @@ const assert = std.debug.assert; const Log2Int = math.Log2Int; /// Returns the logarithm of `x` for the provided `base`, rounding down to the nearest integer. -/// Asserts that `base > 1` and `x != 0`. +/// Asserts that `base > 1` and `x > 0`. pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) { if (@typeInfo(T) != .Int or @typeInfo(T).Int.signedness != .unsigned) @compileError("log_int requires an unsigned integer, found " ++ @typeName(T)); - assert(base > 1 and x != 0); + assert(base > 1 and x > 0); var exponent: Log2Int(T) = 0; var power: T = 1; From 55573d862b14ae5dd1ae96df0553a6c2721da8d8 Mon Sep 17 00:00:00 2001 From: Federico Stra Date: Thu, 14 Sep 2023 17:59:02 +0200 Subject: [PATCH 5/5] std.math.log_int: add comments to prove correctness --- lib/std/math/log_int.zig | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/lib/std/math/log_int.zig b/lib/std/math/log_int.zig index bcdd73cb30bf..d73e273d713b 100644 --- a/lib/std/math/log_int.zig +++ b/lib/std/math/log_int.zig @@ -12,6 +12,24 @@ pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) { assert(base > 1 and x > 0); + // Let's denote by [y] the integer part of y. + + // Throughout the iteration the following invariant is preserved: + // power = base ^ exponent + + // Safety and termination. + // + // We never overflow inside the loop because when we enter the loop we have + // power <= [maxInt(T) / base] + // therefore + // power * base <= maxInt(T) + // is a valid multiplication for type `T` and + // exponent + 1 <= log(base, maxInt(T)) <= log2(maxInt(T)) <= maxInt(Log2Int(T)) + // is a valid addition for type `Log2Int(T)`. + // + // This implies also termination because power is strictly increasing, + // hence it must eventually surpass [x / base] < maxInt(T) and we then exit the loop. + var exponent: Log2Int(T) = 0; var power: T = 1; while (power <= x / base) { @@ -19,6 +37,21 @@ pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) { exponent += 1; } + // If we never entered the loop we must have + // [x / base] < 1 + // hence + // x <= [x / base] * base < base + // thus the result is 0. We can then return exponent, which is still 0. + // + // Otherwise, if we entered the loop at least once, + // when we exit the loop we have that power is exactly divisible by base and + // power / base <= [x / base] < power + // hence + // power <= [x / base] * base <= x < power * base + // This means that + // base^exponent <= x < base^(exponent+1) + // hence the result is exponent. + return exponent; }