authorgravatar for stra.federico@gmail.comFederico Stra <stra.federico@gmail.com> 2023-09-14 21:33:56+02:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2023-09-14 19:33:56+00:00
log4f952c7e0e36dab15f9359f55eb8714f8fe92bcf
tree6f211f6b6a022b3801e3c9a58c13e457a2e9e5e6
parent30e1883834ac630d23cdb0e28793d76219397118
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

std.math.log_int: implement integer logarithm without using float math


3 files changed, 121 insertions(+), 2 deletions(-)

lib/std/math.zig+2
......@@ -241,6 +241,7 @@ pub const log = @import("math/log.zig").log;
241241pub const log2 = @import("math/log2.zig").log2;
242242pub const log10 = @import("math/log10.zig").log10;
243243pub const log10_int = @import("math/log10.zig").log10_int;
244pub const log_int = @import("math/log_int.zig").log_int;
244245pub const log1p = @import("math/log1p.zig").log1p;
245246pub const asinh = @import("math/asinh.zig").asinh;
246247pub const acosh = @import("math/acosh.zig").acosh;
......@@ -362,6 +363,7 @@ test {
362363 _ = log2;
363364 _ = log10;
364365 _ = log10_int;
366 _ = log_int;
365367 _ = log1p;
366368 _ = asinh;
367369 _ = acosh;
lib/std/math/log.zig+5-2
......@@ -23,14 +23,17 @@ pub fn log(comptime T: type, base: T, x: T) T {
2323 .ComptimeFloat => {
2424 return @as(comptime_float, @log(@as(f64, x)) / @log(float_base));
2525 },
26
27 // TODO: implement integer log without using float math.
28 // The present implementation is incorrect, for example
29 // `log(comptime_int, 9, 59049)` should return `5` and not `4`.
2630 .ComptimeInt => {
2731 return @as(comptime_int, @floor(@log(@as(f64, x)) / @log(float_base)));
2832 },
2933
30 // TODO implement integer log without using float math
3134 .Int => |IntType| switch (IntType.signedness) {
3235 .signed => @compileError("log not implemented for signed integers"),
33 .unsigned => return @as(T, @intFromFloat(@floor(@log(@as(f64, @floatFromInt(x))) / @log(float_base)))),
36 .unsigned => return @as(T, math.log_int(T, base, x)),
3437 },
3538
3639 .Float => {
lib/std/math/log_int.zig created+114
......@@ -0,0 +1,114 @@
1const std = @import("../std.zig");
2const math = std.math;
3const testing = std.testing;
4const assert = std.debug.assert;
5const Log2Int = math.Log2Int;
6
7/// Returns the logarithm of `x` for the provided `base`, rounding down to the nearest integer.
8/// Asserts that `base > 1` and `x > 0`.
9pub fn log_int(comptime T: type, base: T, x: T) Log2Int(T) {
10 if (@typeInfo(T) != .Int or @typeInfo(T).Int.signedness != .unsigned)
11 @compileError("log_int requires an unsigned integer, found " ++ @typeName(T));
12
13 assert(base > 1 and x > 0);
14
15 // Let's denote by [y] the integer part of y.
16
17 // Throughout the iteration the following invariant is preserved:
18 // power = base ^ exponent
19
20 // Safety and termination.
21 //
22 // We never overflow inside the loop because when we enter the loop we have
23 // power <= [maxInt(T) / base]
24 // therefore
25 // power * base <= maxInt(T)
26 // is a valid multiplication for type `T` and
27 // exponent + 1 <= log(base, maxInt(T)) <= log2(maxInt(T)) <= maxInt(Log2Int(T))
28 // is a valid addition for type `Log2Int(T)`.
29 //
30 // This implies also termination because power is strictly increasing,
31 // hence it must eventually surpass [x / base] < maxInt(T) and we then exit the loop.
32
33 var exponent: Log2Int(T) = 0;
34 var power: T = 1;
35 while (power <= x / base) {
36 power *= base;
37 exponent += 1;
38 }
39
40 // If we never entered the loop we must have
41 // [x / base] < 1
42 // hence
43 // x <= [x / base] * base < base
44 // thus the result is 0. We can then return exponent, which is still 0.
45 //
46 // Otherwise, if we entered the loop at least once,
47 // when we exit the loop we have that power is exactly divisible by base and
48 // power / base <= [x / base] < power
49 // hence
50 // power <= [x / base] * base <= x < power * base
51 // This means that
52 // base^exponent <= x < base^(exponent+1)
53 // hence the result is exponent.
54
55 return exponent;
56}
57
58test "math.log_int" {
59 // Test all unsigned integers with 2, 3, ..., 64 bits.
60 // We cannot test 0 or 1 bits since base must be > 1.
61 inline for (2..64 + 1) |bits| {
62 const T = @Type(std.builtin.Type{
63 .Int = std.builtin.Type.Int{ .signedness = .unsigned, .bits = @intCast(bits) },
64 });
65
66 // for base = 2, 3, ..., min(maxInt(T),1024)
67 var base: T = 1;
68 while (base < math.maxInt(T) and base <= 1024) {
69 base += 1;
70
71 // test that `log_int(T, base, 1) == 0`
72 try testing.expectEqual(@as(Log2Int(T), 0), log_int(T, base, 1));
73
74 // For powers `pow = base^exp > 1` that fit inside T,
75 // test that `log_int` correctly detects the jump in the logarithm
76 // from `log(pow-1) == exp-1` to `log(pow) == exp`.
77 var exp: Log2Int(T) = 0;
78 var pow: T = 1;
79 while (pow <= math.maxInt(T) / base) {
80 exp += 1;
81 pow *= base;
82
83 try testing.expectEqual(exp - 1, log_int(T, base, pow - 1));
84 try testing.expectEqual(exp, log_int(T, base, pow));
85 }
86 }
87 }
88}
89
90test "math.log_int vs math.log2" {
91 const types = [_]type{ u2, u3, u4, u8, u16 };
92 inline for (types) |T| {
93 var n: T = 0;
94 while (n < math.maxInt(T)) {
95 n += 1;
96 const special = math.log2_int(T, n);
97 const general = log_int(T, 2, n);
98 try testing.expectEqual(special, general);
99 }
100 }
101}
102
103test "math.log_int vs math.log10" {
104 const types = [_]type{ u4, u5, u6, u8, u16 };
105 inline for (types) |T| {
106 var n: T = 0;
107 while (n < math.maxInt(T)) {
108 n += 1;
109 const special = math.log10_int(n);
110 const general = log_int(T, 10, n);
111 try testing.expectEqual(special, general);
112 }
113 }
114}