bigint: log_n for bases that fit within one DIGIT or are a power of two.

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() {

core/math/bigint/log.odin

@@ -9,20 +9,27 @@ 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]);
+		Fast path for `Int`s that fit within a single `DIGIT`.
-   //    return MP_OKAY;
+	*/
-   // }
+	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);
@@ -31,6 +38,8 @@ log_n :: proc(a: ^Int, base: int) -> (log: int, err: Error) {
 	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;
+   	}
+}

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) {
+	if radix < 2 || radix > 64 {
-   // mp_err err;
+		return -1, .Invalid_Input;
-   // mp_int a_;
+	}
-   // int b;
 
-   // /* make sure the radix is in range */
+ 	if is_zero(a) {
-   // if ((radix < 2) || (radix > 64)) {
+ 		return 1, .OK;
-   //    return MP_VAL;
+ 	}
-   // }
 
-   // if (mp_iszero(a)) {
+ 	t.sign = .Zero_or_Positive;
-   //    *size = 2;
+ 	log: int;
-   //    return MP_OKAY;
-   // }
 
-   // a_ = *a;
+ 	log, err = log_n(t, radix);
-   // a_.sign = MP_ZPOS;
+ 	if err != .OK {
-   // if ((err = mp_log_n(&a_, radix, &b)) != MP_OKAY) {
+ 		return log, err;
-   //    return 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);
+		log truncates to zero, so we need to add one more, and one for `-` if negative.
-
+ 	*/
-   return size, .OK;
+ 	if is_neg(a) {
+ 		return log + 2, .OK;
+ 	} else {
+ 		return log + 1, .OK;
+ 	}
 }