diff --git a/core/math/bigint/example.odin b/core/math/bigint/example.odin index 2d3f5e12d..e21925f09 100644 --- a/core/math/bigint/example.odin +++ b/core/math/bigint/example.odin @@ -76,7 +76,7 @@ demo :: proc() { fmt.printf("c: %v, bits: %v\n", cs, count_bits(c)); delete(as); delete(bs); delete(cs); - fmt.println("log2:", log_n(a, 8)); + fmt.println("radix_size:", radix_size(a, 10)); } main :: proc() { diff --git a/core/math/bigint/log.odin b/core/math/bigint/log.odin index fd8116b78..347e60d99 100644 --- a/core/math/bigint/log.odin +++ b/core/math/bigint/log.odin @@ -9,28 +9,37 @@ package bigint The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks. */ -log_n :: proc(a: ^Int, base: int) -> (log: int, err: Error) { +import "core:fmt" + +log_n_int :: proc(a: ^Int, base: int) -> (log: int, err: Error) { assert_initialized(a); if is_neg(a) || is_zero(a) || base < 2 || DIGIT(base) > _DIGIT_MAX { return -1, .Invalid_Input; } + /* + Fast path for bases that are a power of two. + */ if is_power_of_two(base) { return _log_power_of_two(a, base), .OK; } - // if (MP_HAS(S_MP_LOG_D) && (a->used == 1)) { - // *c = s_mp_log_d((mp_digit)base, a->dp[0]); - // return MP_OKAY; - // } + /* + Fast path for `Int`s that fit within a single `DIGIT`. + */ + if a.used == 1 { + return log_n_digit(a.digit[0], DIGIT(base)), .OK; + } - // if (MP_HAS(S_MP_LOG)) { - // return s_mp_log(a, (mp_digit)base, c); - // } + // if (MP_HAS(S_MP_LOG)) { + // return s_mp_log(a, (mp_digit)base, c); + // } return -1, .Unimplemented; } +log_n :: proc{log_n_int, log_n_digit}; + /* Returns the log2 of an `Int`, provided `base` is a power of two. Don't call it if it isn't. @@ -44,3 +53,73 @@ _log_power_of_two :: proc(a: ^Int, base: int) -> (log: int) { } return (count_bits(a) - 1) / y; } + +/* + +*/ +small_pow :: proc(base: _WORD, exponent: _WORD) -> (result: _WORD) { + exponent := exponent; base := base; + result = _WORD(1); + + for exponent != 0 { + if exponent & 1 == 1 { + result *= base; + } + exponent >>= 1; + base *= base; + } + return result; +} + +log_n_digit :: proc(a: DIGIT, base: DIGIT) -> (log: int) { + /* + If the number is smaller than the base, it fits within a fraction. + Therefore, we return 0. + */ + if a < base { + return 0; + } + + /* + If a number equals the base, the log is 1. + */ + if a == base { + return 1; + } + + N := _WORD(a); + bracket_low := _WORD(1); + bracket_high := _WORD(base); + high := 1; + low := 0; + + for bracket_high < N { + low = high; + bracket_low = bracket_high; + high <<= 1; + bracket_high *= bracket_high; + } + + for high - low > 1 { + mid := (low + high) >> 1; + bracket_mid := bracket_low * small_pow(_WORD(base), _WORD(mid - low)); + + if N < bracket_mid { + high = mid; + bracket_high = bracket_mid; + } + if N > bracket_mid { + low = mid; + bracket_low = bracket_mid; + } + if N == bracket_mid { + return mid; + } + } + + if bracket_high == N { + return high; + } else { + return low; + } +} \ No newline at end of file diff --git a/core/math/bigint/radix.odin b/core/math/bigint/radix.odin index 5bfb6b9a0..64118f70b 100644 --- a/core/math/bigint/radix.odin +++ b/core/math/bigint/radix.odin @@ -26,8 +26,7 @@ itoa :: proc(a: ^Int, radix: int, allocator := context.allocator) -> (res: strin /* Fast path for radixes that are a power of two. */ - if radix & 1 == 0 { - + if is_power_of_two(radix) { } @@ -51,30 +50,34 @@ itoa :: proc(a: ^Int, radix: int, allocator := context.allocator) -> (res: strin int_to_string :: itoa; +/* + We size for `string`, not `cstring`. +*/ +radix_size :: proc(a: ^Int, radix: int) -> (size: int, err: Error) { + t := a; -radix_size :: proc(a: ^Int, base: int) -> (size: int, err: Error) { - // mp_err err; - // mp_int a_; - // int b; + if radix < 2 || radix > 64 { + return -1, .Invalid_Input; + } - // /* make sure the radix is in range */ - // if ((radix < 2) || (radix > 64)) { - // return MP_VAL; - // } + if is_zero(a) { + return 1, .OK; + } - // if (mp_iszero(a)) { - // *size = 2; - // return MP_OKAY; - // } + t.sign = .Zero_or_Positive; + log: int; - // a_ = *a; - // a_.sign = MP_ZPOS; - // if ((err = mp_log_n(&a_, radix, &b)) != MP_OKAY) { - // return err; - // } + log, err = log_n(t, radix); + if err != .OK { + return log, err; + } - // /* mp_ilogb truncates to zero, hence we need one extra put on top and one for `\0`. */ - // *size = (size_t)b + 2U + (mp_isneg(a) ? 1U : 0U); - - return size, .OK; + /* + log truncates to zero, so we need to add one more, and one for `-` if negative. + */ + if is_neg(a) { + return log + 2, .OK; + } else { + return log + 1, .OK; + } } \ No newline at end of file