authorgravatar for git@adrianjdelgado.comAdrian Delgado <git@adrianjdelgado.com> 2023-03-07 20:41:40-05:00
committergravatar for git@adrianjdelgado.comAdrian Delgado <git@adrianjdelgado.com> 2023-03-22 12:03:11-05:00
logf34aa2d14fbacc7e4d2471c8b4897da257cffc03
tree777e30b6ba6f59b90522819f0b0286b5cc01d272
parentd61ac0db8c62f706ea65b70d2772cbb8c4efb416

std.math: port int log10 from Rust

2 files changed, 146 insertions(+), 1 deletions(-)

lib/std/math.zig+1
......@@ -250,6 +250,7 @@ pub const ln = @import("math/ln.zig").ln;
250250pub const log = @import("math/log.zig").log;
251251pub const log2 = @import("math/log2.zig").log2;
252252pub const log10 = @import("math/log10.zig").log10;
253pub const log10_int = @import("math/log10.zig").log10_int;
253254pub const log1p = @import("math/log1p.zig").log1p;
254255pub const asinh = @import("math/asinh.zig").asinh;
255256pub const acosh = @import("math/acosh.zig").acosh;
lib/std/math/log10.zig+145-1
......@@ -2,6 +2,8 @@ const std = @import("../std.zig");
22const math = std.math;
33const testing = std.testing;
44const maxInt = std.math.maxInt;
5const assert = std.debug.assert;
6const Log2Int = std.math.Log2Int;
57
68/// Returns the base-10 logarithm of x.
79///
......@@ -22,8 +24,150 @@ pub fn log10(x: anytype) @TypeOf(x) {
2224 },
2325 .Int => |IntType| switch (IntType.signedness) {
2426 .signed => @compileError("log10 not implemented for signed integers"),
25 .unsigned => return @floatToInt(T, @floor(@log10(@intToFloat(f64, x)))),
27 .unsigned => return log10_int(x),
2628 },
2729 else => @compileError("log10 not implemented for " ++ @typeName(T)),
2830 }
2931}
32
33// Based on Rust, which is licensed under the MIT license.
34// https://github.com/rust-lang/rust/blob/f63ccaf25f74151a5d8ce057904cd944074b01d2/LICENSE-MIT
35//
36// https://github.com/rust-lang/rust/blob/f63ccaf25f74151a5d8ce057904cd944074b01d2/library/core/src/num/int_log10.rs
37
38/// Return the log base 10 of integer value x, rounding down to the
39/// nearest integer.
40pub fn log10_int(x: anytype) Log2Int(@TypeOf(x)) {
41 const T = @TypeOf(x);
42 const OutT = Log2Int(T);
43 if (@typeInfo(T) != .Int or @typeInfo(T).Int.signedness != .unsigned)
44 @compileError("log10_int requires an unsigned integer, found " ++ @typeName(T));
45
46 assert(x != 0);
47
48 const bit_size = @typeInfo(T).Int.bits;
49
50 if (bit_size <= 8) {
51 return @intCast(OutT, log10_int_u8(x));
52 } else if (bit_size <= 16) {
53 return @intCast(OutT, less_than_5(x));
54 }
55
56 var val = x;
57 var log: u32 = 0;
58
59 inline for (0..11) |i| {
60 // Unnecesary branches should be removed by the compiler
61 if (bit_size > (1 << (11 - i)) * 5 * @log2(10.0) and val >= pow10((1 << (11 - i)) * 5)) {
62 const num_digits = (1 << (11 - i)) * 5;
63 val /= pow10(num_digits);
64 log += num_digits;
65 }
66 }
67
68 if (val >= pow10(5)) {
69 val /= pow10(5);
70 log += 5;
71 }
72
73 return @intCast(OutT, log + less_than_5(@intCast(u32, val)));
74}
75
76fn pow10(comptime y: comptime_int) comptime_int {
77 var result = 1;
78 for (0..y) |_| {
79 result *= 10;
80 }
81 return result;
82}
83
84inline fn log10_int_u8(x: u8) u32 {
85 // For better performance, avoid branches by assembling the solution
86 // in the bits above the low 8 bits.
87
88 // Adding c1 to val gives 10 in the top bits for val < 10, 11 for val >= 10
89 const C1: u32 = 0b11_00000000 - 10; // 758
90 // Adding c2 to val gives 01 in the top bits for val < 100, 10 for val >= 100
91 const C2: u32 = 0b10_00000000 - 100; // 412
92
93 // Value of top bits:
94 // +c1 +c2 1&2
95 // 0..=9 10 01 00 = 0
96 // 10..=99 11 01 01 = 1
97 // 100..=255 11 10 10 = 2
98 return ((x + C1) & (x + C2)) >> 8;
99}
100
101inline fn less_than_5(x: u32) u32 {
102 // Similar to log10u8, when adding one of these constants to val,
103 // we get two possible bit patterns above the low 17 bits,
104 // depending on whether val is below or above the threshold.
105 const C1: u32 = 0b011_00000000000000000 - 10; // 393206
106 const C2: u32 = 0b100_00000000000000000 - 100; // 524188
107 const C3: u32 = 0b111_00000000000000000 - 1000; // 916504
108 const C4: u32 = 0b100_00000000000000000 - 10000; // 514288
109
110 // Value of top bits:
111 // +c1 +c2 1&2 +c3 +c4 3&4 ^
112 // 0..=9 010 011 010 110 011 010 000 = 0
113 // 10..=99 011 011 011 110 011 010 001 = 1
114 // 100..=999 011 100 000 110 011 010 010 = 2
115 // 1000..=9999 011 100 000 111 011 011 011 = 3
116 // 10000..=99999 011 100 000 111 100 100 100 = 4
117 return (((x + C1) & (x + C2)) ^ ((x + C3) & (x + C4))) >> 17;
118}
119
120fn oldlog10(x: anytype) u8 {
121 return @floatToInt(u8, @log10(@intToFloat(f64, x)));
122}
123
124test "oldlog10 doesn't work" {
125 try testing.expect(14 != oldlog10(pow10(15) - 1));
126
127 // log10(10**15 -1) should indeed be 14
128 try testing.expect(14 == log10_int(@as(u64, pow10(15) - 1)));
129}
130
131test "log10_int vs old implementation" {
132 const int_types = .{ u8, u16, u32, u64, u128 };
133
134 inline for (int_types) |T| {
135 const last = @min(maxInt(T), 100_000);
136 for (1..last) |i| {
137 const x = @intCast(T, i);
138 try testing.expectEqual(oldlog10(x), log10_int(x));
139 }
140
141 const max_int: T = maxInt(T);
142 try testing.expectEqual(oldlog10(max_int), log10_int(max_int));
143 }
144}
145
146test "log10_int close to powers of 10" {
147 const int_types = .{ u8, u16, u32, u64, u128, u256, u512 };
148 const max_log_values: [7]usize = .{ 2, 4, 9, 19, 38, 77, 154 };
149
150 inline for (int_types, max_log_values) |T, expected_max_ilog| {
151 const max_val: T = maxInt(T);
152
153 try testing.expectEqual(expected_max_ilog, log10_int(max_val));
154
155 for (0..(expected_max_ilog + 1)) |idx| {
156 const i = @intCast(T, idx);
157 const p: T = try math.powi(T, 10, i);
158
159 const b = @intCast(Log2Int(T), i);
160
161 if (p >= 10) {
162 try testing.expectEqual(b - 1, log10_int(p - 9));
163 try testing.expectEqual(b - 1, log10_int(p - 1));
164 }
165
166 try testing.expectEqual(b, log10_int(p));
167 try testing.expectEqual(b, log10_int(p + 1));
168 if (p >= 10) {
169 try testing.expectEqual(b, log10_int(p + 9));
170 }
171 }
172 }
173}