diff --git a/arbi/src/base.rs b/arbi/src/base.rs index 1d28934..31df503 100644 --- a/arbi/src/base.rs +++ b/arbi/src/base.rs @@ -1,5 +1,5 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */ @@ -159,6 +159,7 @@ impl Base { /// /// let b36 = Base::new(36).unwrap(); /// ``` + #[inline] pub fn new(value: u8) -> Result { value.try_into() } @@ -172,6 +173,7 @@ impl Base { /// let b10 = Base::new(10).unwrap(); /// assert_eq!(b10.value(), 10); /// ``` + #[inline] pub fn value(self) -> u8 { self.0 } diff --git a/arbi/src/bits.rs b/arbi/src/bits.rs index 4ce21a9..ffe5b58 100644 --- a/arbi/src/bits.rs +++ b/arbi/src/bits.rs @@ -368,6 +368,7 @@ mod tests { #[cfg(test)] mod random { use super::*; + use alloc::vec; macro_rules! test_bit_ops_for_type { ($rng:expr, $die:expr) => { diff --git a/arbi/src/builtin_int_methods/trailing_zeros.rs b/arbi/src/builtin_int_methods/trailing_zeros.rs index b1d99b9..6bfa7f5 100644 --- a/arbi/src/builtin_int_methods/trailing_zeros.rs +++ b/arbi/src/builtin_int_methods/trailing_zeros.rs @@ -33,6 +33,7 @@ mod tests { use crate::util::test::{get_seedable_rng, get_uniform_die, Distribution}; use crate::Arbi; use crate::{BitCount, DDigit, Digit}; + use alloc::vec; macro_rules! assert_trailing_zeros { ($value:expr) => { diff --git a/arbi/src/floor.rs b/arbi/src/floor.rs deleted file mode 100644 index ec7bd2e..0000000 --- a/arbi/src/floor.rs +++ /dev/null @@ -1,157 +0,0 @@ -/* -Copyright 2024 Owain Davies -SPDX-License-Identifier: Apache-2.0 OR MIT -*/ - -use core::f64; - -/// Replacement for f64::floor() in no-std. -/// -/// The motivation behind implementing this is that it was the only function -/// when converting from std to no-std that was not available in core. Floor -/// might not be needed at all in the future. -pub(crate) fn floor(x: f64) -> f64 { - const MANTISSA_MASK: u64 = 0x000fffffffffffff_u64; - const ABSOLUTE_MASK: u64 = 0x7fffffffffffffff_u64; - const IMPLICIT_ONE: i64 = 0x0010000000000000_i64; - - let mut x_i64: i64 = x.to_bits() as i64; - let unbiased_exponent: i32 = (0x7ff & (x_i64 >> 52)) as i32 - 0x3ff; - - match unbiased_exponent { - 0x400 => { - assert!(x.is_infinite() || x.is_nan()); - x + x - } // f - - exp if exp < 0 => { - if x_i64 >= 0 { - x_i64 = 0; // a - } else if (ABSOLUTE_MASK as i64 & x_i64) != 0 { - x_i64 = 0xbff0000000000000_u64 as i64; // b - } - - f64::from_bits(x_i64 as u64) - } - - exp if (0..52).contains(&exp) => { - let m = MANTISSA_MASK >> exp; - if (m as i64 & x_i64) == 0 { - x // c - } else { - if x_i64 < 0 { - x_i64 += IMPLICIT_ONE >> exp as i64; // d - } - x_i64 &= !(m as i64); - - f64::from_bits(x_i64 as u64) // e - } - } - - _ => x, // g - } -} - -#[cfg(test)] -mod tests { - #[cfg(test)] - extern crate std; - - use super::*; - use crate::util::test::{get_seedable_rng, Rng}; - use std::f64; - - #[test] - fn test_floor_marked_code_paths() { - // a - assert_eq!(floor(0.0_f64), (0.0_f64).floor()); - - // b - assert_eq!( - floor(-0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000015134722735375557_f64), - (-0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000015134722735375557_f64).floor() - ); - - // c - assert_eq!(floor(33.0_f64), (33.0_f64).floor()); - - // d - assert_eq!(floor(-2.5_f64), (-2.5_f64).floor()); - - // e - assert_eq!(floor(21.192_f64), (21.192_f64).floor()); - - // f - assert_eq!(floor(f64::INFINITY), (f64::INFINITY).floor()); - assert_eq!(floor(f64::NEG_INFINITY), (f64::NEG_INFINITY).floor()); - assert!(floor(f64::NAN).is_nan() && (f64::NAN).floor().is_nan()); - - // g - assert_eq!( - floor(9007199254740992.0_f64), - (9007199254740992.0_f64).floor() - ); - } - - #[test] - fn test_floor_against_std() { - let test_values = [ - 0.0, - -0.0, - 1.0, - -1.0, - 2.5, - -2.5, - 3.9999999999, - -3.9999999999, - 987654321.987654321, - -987654321.987654321, - f64::NAN, - f64::INFINITY, - f64::NEG_INFINITY, - f64::MIN_POSITIVE, - f64::MIN_POSITIVE / (1u64 << 52) as f64, - (1_i64 << 53) as f64, - -(1_i64 << 53) as f64, - f64::MAX, - -f64::MAX, - f64::MIN, - -f64::MIN, - ]; - - for &value in &test_values { - let expected = value.floor(); - let result = floor(value); - if expected.is_nan() { - assert!(result.is_nan(), "Failed on value: {}", value); - continue; - } - assert_eq!( - expected.to_bits(), - result.to_bits(), - "Failed on value: {}", - value - ); - } - - let (mut rng, _) = get_seedable_rng(); - - for _ in 0..i16::MAX { - let r_u64: u64 = rng.gen(); - let r = f64::from_bits(r_u64); - - let expected = r.floor(); - let result = floor(r); - if expected.is_nan() { - assert!(result.is_nan(), "Failed on value: {}", r); - continue; - } - assert_eq!( - expected.to_bits(), - result.to_bits(), - "Failed on value: {}", - r - ); - } - } -} diff --git a/arbi/src/lib.rs b/arbi/src/lib.rs index ed96871..e800a95 100644 --- a/arbi/src/lib.rs +++ b/arbi/src/lib.rs @@ -33,7 +33,6 @@ mod division; pub mod docs; mod exponentiation; mod fits; -mod floor; mod fmt; mod from_double; mod from_integral; @@ -96,8 +95,6 @@ type QDigit = u128; #[allow(dead_code)] type SQDigit = i128; -const DBL_MAX_INT: u64 = 0x20000000000000; // 2 ** 53 - /// Arbitrary Precision Integer type. /// /// # Internal Representation diff --git a/arbi/src/ops/bit/impls.rs b/arbi/src/ops/bit/impls.rs index d9550db..e1b75f9 100644 --- a/arbi/src/ops/bit/impls.rs +++ b/arbi/src/ops/bit/impls.rs @@ -377,6 +377,7 @@ mod tests { use super::BitwiseOp; use crate::util::test::{get_seedable_rng, get_uniform_die, Distribution}; use crate::{Arbi, DDigit, Digit}; + use alloc::{format, vec}; use core::ops::RangeInclusive; struct BitwiseTestCase { diff --git a/arbi/src/right_shift.rs b/arbi/src/right_shift.rs index 6196781..e54bc44 100644 --- a/arbi/src/right_shift.rs +++ b/arbi/src/right_shift.rs @@ -214,6 +214,7 @@ mod test_arithmetic_rshift { use crate::{ Arbi, Assign, BitCount, DDigit, Digit, SDDigit, SDigit, SQDigit, }; + use alloc::vec; #[test] fn test_mark_a() { diff --git a/arbi/src/to_double.rs b/arbi/src/to_double.rs index b48c093..23bbc64 100644 --- a/arbi/src/to_double.rs +++ b/arbi/src/to_double.rs @@ -1,5 +1,5 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */ @@ -61,7 +61,7 @@ mod tests { use super::*; use crate::util::test::float::*; use crate::util::test::{get_seedable_rng, get_uniform_die, Distribution}; - use crate::{DDigit, QDigit, SDDigit, SQDigit, DBL_MAX_INT}; + use crate::{DDigit, QDigit, SDDigit, SQDigit}; #[test] fn smoke() { diff --git a/arbi/src/to_string.rs b/arbi/src/to_string.rs index 78fba19..f03cd57 100644 --- a/arbi/src/to_string.rs +++ b/arbi/src/to_string.rs @@ -1,88 +1,225 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */ use crate::from_string::configs::BASE_MBS; use crate::Base; -use crate::BitCount; -use crate::{Arbi, Digit, DBL_MAX_INT}; +use crate::{Arbi, Digit}; use alloc::string::String; +use alloc::vec::Vec; use core::convert::TryInto; -/// `LOG_BASE_2[base]` gives \\( \log_{base}2 \\) (with some rounding up). -pub(crate) const LOG_BASE_2: [f64; 37] = [ - 0.0, 0.0, 1.0, 0.631, 0.5, 0.431, 0.387, 0.357, 0.334, 0.316, 0.302, 0.290, - 0.279, 0.271, 0.263, 0.256, 0.25, 0.245, 0.240, 0.236, 0.232, 0.228, 0.225, - 0.222, 0.219, 0.216, 0.213, 0.211, 0.209, 0.206, 0.204, 0.202, 0.200, - 0.199, 0.197, 0.195, 0.194, -]; +const BASE_DIGITS_LOWER_BYTES: &[u8; 36] = + b"0123456789abcdefghijklmnopqrstuvwxyz"; +const BASE_DIGITS_UPPER_BYTES: &[u8; 36] = + b"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; impl Arbi { - // In the documentation block above `base_length()`, it is shown how to - // calculate the minimum number of base-b digits, \\( n \\), required to - // represent any integer with \\( m \\) base-c digits. - // - // We can find an upper bound another way as well. - // - // Let \\( e \\) denote the largest exponent such that \\( c^{e} \\) fits - // in one base-b digit. By construction, any e-digit base-c number will fit - // in a base-b digit, but not every \\( (e + 1) \\) digit base-c number - // will. By dividing the number of base-c digits in an integer, \\( m \\), - // by \\( e \\), and rounding up the result of the division to the nearest - // integer, we obtain an upperbound on the number of base-b digits, \\( n - // \\), needed to represent any m-digit base-c integer. - - /// Return an estimate of the number of base-`base` digits required to - /// represent this integer. - /// - /// Given an n-digit integer in base b, the largest integer representable is - /// \\[ - /// b^{n} - 1 - /// \\] - /// - /// The minimum number of base b digits, n, required to represent any - /// m-digit base c integer is such that: - /// \\[ - /// \begin{align} - /// b^{n} - 1 & \geq c^{m} - 1 \\\\ - /// b^{n} & \geq c^{m} \\\\ - /// \frac{b^{n}}{c^{m}} & \geq 1 \\\\ - /// \log(b^{n}) - \log(c^{m}) & \geq 0 \\\\ - /// n & \geq m \cdot \frac{\log(c)}{\log(b)} - /// \end{align} - /// \\] - /// - /// Use - /// \\[ - /// n = \left\lceil m \cdot \frac{\log(c)}{\log(b)} \right\rceil - /// \\] - /// - /// For example, the minimum number of base 10 digits required to represent - /// any m-digit base 2 integer is: - /// \\[ - /// n = \left\lceil m * \log_{10}(2) \right\rceil - /// \\] - /// If `x` has bit length less than or equal to 2 ** 53, return an estimate - /// of the number of base-`base` digits needed to represent the absolute - /// value of this integer. Otherwise, return the exact number of base-`base` - /// digits. - fn base_length(x: &Self, base: usize) -> BitCount { - // TODO: analyze - if x == 0 { - return 1; + fn to_string_base_pow2(&self, base: Base, lowercase: bool) -> String { + debug_assert!(base.value().is_power_of_two()); + + let base_digits = if lowercase { + BASE_DIGITS_LOWER_BYTES + } else { + BASE_DIGITS_UPPER_BYTES + }; + + if self.is_zero() { + return "0".into(); } - let bitlen: BitCount = x.size_bits(); - if bitlen as BitCount > DBL_MAX_INT as BitCount { - // TODO: find some quick upperbound. - // let ilog2_base = base.ilog2(); - // (x.size_bits() - 1) / (ilog2_base as BitCount) + (1 as BitCount) - x.size_base_ref(base as u32) + + let capacity = self + .size_base_ref(base.value() as u32) + .try_into() + .unwrap_or(usize::MAX) + .saturating_add(usize::from(self.is_negative())); + let mut bytes: Vec = Vec::new(); + bytes.reserve_exact(capacity); + + #[cfg(debug_assertions)] + let initial_capacity = bytes.capacity(); + + let base_digit_bits = base.value().trailing_zeros(); + let base_digit_mask = ((1 as Digit) << base_digit_bits) - 1; + + if matches!(base.value(), 2 | 4 | 16) { + debug_assert!(Digit::BITS % base.value().trailing_zeros() == 0); + self.process_digits_base_pow2_aligned( + &mut bytes, + base_digit_bits, + base_digit_mask, + base_digits, + true, + ); } else { - // This is much more efficient than using size_base() - crate::floor::floor(bitlen as f64 * LOG_BASE_2[base] + 1.0) - as BitCount + self.process_digits_base_pow2_generic( + &mut bytes, + base_digit_bits, + base_digit_mask, + base_digits, + capacity, + ); + } + + #[cfg(debug_assertions)] + { + // Check that no reallocation occurred + debug_assert_eq!(bytes.capacity(), initial_capacity); + // Check that len is equal to requested capacity + debug_assert_eq!(bytes.len(), capacity); + } + + String::from_utf8(bytes).unwrap() + } + + fn process_digits_base_pow2_aligned( + &self, + bytes: &mut Vec, + base_digit_bits: u32, + base_digit_mask: Digit, + base_digits: &[u8; 36], + start_from_msd: bool, + ) { + let batches_per_digit = Digit::BITS / base_digit_bits; + + match start_from_msd { + true => { + // This avoids the need to reverse `bytes` at the end and + // benchmarks show that this is also a little more efficient. + if self.is_negative() { + bytes.push(b'-'); + } + + /* Handle most significant digit specially */ + let msd_idx = self.size() - 1; + let msd = self.vec[msd_idx]; + let msd_bits = Digit::BITS - msd.leading_zeros(); + let msd_bits_mod_base_bits = msd_bits % base_digit_bits; + + // Partial chunk + if msd_bits_mod_base_bits != 0 { + let shift = base_digit_bits - msd_bits_mod_base_bits; + let value = (msd >> (msd_bits - msd_bits_mod_base_bits)) + << shift + >> shift; + bytes.push(base_digits[value as usize]); + } + + // Full chunks + let mut shift = msd_bits - msd_bits_mod_base_bits; + debug_assert!(shift % base_digit_bits == 0); + while shift != 0 { + shift -= base_digit_bits; + bytes.push( + base_digits + [((msd >> shift) & base_digit_mask) as usize], + ); + } + + /* Handle remaining digits (all full chunks) */ + let first_shift = (batches_per_digit - 1) * base_digit_bits; + for &digit in self.vec[..msd_idx].iter().rev() { + bytes.push( + base_digits[((digit >> first_shift) & base_digit_mask) + as usize], + ); + + let mut shift = first_shift; + debug_assert!(shift % base_digit_bits == 0); + while shift != 0 { + shift -= base_digit_bits; + bytes.push( + base_digits + [((digit >> shift) & base_digit_mask) as usize], + ); + } + } + } + false => { + let mut j = 0; + let last_idx = self.size() - 1; + while j < last_idx { + let mut digit = self.vec[j]; + for _ in 0..batches_per_digit { + bytes.push( + base_digits[(digit & base_digit_mask) as usize], + ); + digit >>= base_digit_bits; + } + j += 1; + } + + // Handle last digit specially to avoid pushing leading zeros + let mut digit = self.vec[last_idx]; + while digit != 0 { + bytes.push(base_digits[(digit & base_digit_mask) as usize]); + digit >>= base_digit_bits; + } + + if self.is_negative() { + bytes.push(b'-'); + } + + bytes.reverse(); + } + } + } + + fn process_digits_base_pow2_generic( + &self, + bytes: &mut Vec, + base_digit_bits: u32, + base_digit_mask: Digit, + base_digits: &[u8; 36], + num_bytes: usize, + ) { + // TODO: vec![0; capacity] might be better rather than allocating + // capacity, then resizing, like we are doing now. + bytes.resize(num_bytes, 0); + + if self.is_negative() { + bytes[0] = b'-'; + } + + /* Process base_digit_bits batches */ + let mut output_idx = num_bytes - 1; + let mut j = 0; + let mut batch_stop_idx = 0; + let last_idx = self.size() - 1; + while j < last_idx { + let mut cur = self.vec[j] >> batch_stop_idx; + batch_stop_idx += base_digit_bits; + + if Digit::BITS <= batch_stop_idx { + // For < case, batch crosses digit boundaries (takes bits from + // both the current and next digit). + j += 1; + batch_stop_idx -= Digit::BITS; + + // Technically, this is not needed for the = case, but we do + // need to update counters, as above. + cur |= self.vec[j] << (base_digit_bits - batch_stop_idx); + } + + bytes[output_idx] = base_digits[(cur & base_digit_mask) as usize]; + output_idx = output_idx.wrapping_sub(1); + } + + // Handle last digit + let last_digit = self.vec[last_idx]; + let msb = Digit::BITS - last_digit.leading_zeros(); + while batch_stop_idx < msb { + let cur = (last_digit >> batch_stop_idx) & base_digit_mask; + bytes[output_idx] = base_digits[cur as usize]; + output_idx = output_idx.wrapping_sub(1); + batch_stop_idx += base_digit_bits; } + + debug_assert_eq!( + output_idx, + if self.is_negative() { 0 } else { usize::MAX } + ); } pub(crate) fn to_string_base_( @@ -93,45 +230,36 @@ impl Arbi { let base: usize = base.value() as usize; assert!((2..=36).contains(&base)); - const BASE_DIGITS_LOWER: &str = "0123456789abcdefghijklmnopqrstuvwxyz"; - const BASE_DIGITS_UPPER: &str = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; - let base_digits = if lowercase { - BASE_DIGITS_LOWER + BASE_DIGITS_LOWER_BYTES } else { - BASE_DIGITS_UPPER + BASE_DIGITS_UPPER_BYTES }; if self.size() == 0 { return "0".into(); } - let mut result = String::new(); - let true_estimate: BitCount = - Self::base_length(self, base) + BitCount::from(self.neg); - let estimate: usize = if true_estimate > usize::MAX as BitCount { - let exact = - self.size_base_ref(base as u32) + BitCount::from(self.neg); - assert!(exact > isize::MAX as BitCount); - panic!( - "Base-{} digit estimation exceeds isize::MAX bytes. Exact = {}", - base, exact - ); + /* Allocate memory for the result. This will be exactly the number of + bytes needed or higher by one. */ + let mut bytes = Vec::new(); + let capacity: usize = if base.is_power_of_two() { + // Exact number of bytes needed. + self.size_base_ref(base as u32) } else { - true_estimate.try_into().unwrap() - }; - result.reserve(estimate); - // let exact_chars = - // self.size_base(base as u32) + if self.neg { 1 } else { 0 }; - // let exact_chars_usize = if exact_chars > usize::MAX as BitCount { - // panic!( - // "More than usize::MAX bytes needed. Arbi = {}, base = {}", - // self, base - // ); - // } else { - // exact_chars.try_into().unwrap() - // }; - // result.reserve(exact_chars_usize); + // Exact number of bytes needed or one more. + Self::size_base_with_size_bits_maybe_over_by_one( + base as u32, + self.size_bits(), + ) + } + .try_into() + .unwrap_or(usize::MAX) + .saturating_add(usize::from(self.is_negative())); + bytes.reserve_exact(capacity); + + #[cfg(debug_assertions)] + let initial_capacity = bytes.capacity(); let basembs = BASE_MBS[base]; let max_batch_size = basembs.mbs; @@ -148,17 +276,30 @@ impl Arbi { let current_digit: Digit = remainder % base as Digit; remainder /= base as Digit; - result.push( - base_digits.chars().nth(current_digit as usize).unwrap(), - ); + bytes.push(base_digits[current_digit as usize]); } } - if self.neg { - result.push('-'); + if self.is_negative() { + bytes.push(b'-'); } - result.chars().rev().collect::() + #[cfg(debug_assertions)] + { + // Check that no reallocation occurred + debug_assert_eq!(bytes.capacity(), initial_capacity); + // Check that the requested capacity was exact or higher by one + debug_assert!( + bytes.len() == capacity || bytes.len() == capacity - 1, + "Capacity estimate {} should be exact or higher by one than \ + the true value {}", + capacity, + bytes.len() + ); + } + + bytes.reverse(); + String::from_utf8(bytes).unwrap() } /// Return a [`String`] containing the base-`base` representation of the @@ -174,8 +315,14 @@ impl Arbi { /// assert_eq!(Arbi::from(123456789).to_string_base(HEX), "75bcd15"); /// assert_eq!(Arbi::from(-123456789).to_string_base(HEX), "-75bcd15"); /// ``` + #[inline] pub fn to_string_base(&self, base: Base) -> String { - self.to_string_base_(base, true) + let lowercase = true; + if base.value().is_power_of_two() { + self.to_string_base_pow2(base, lowercase) + } else { + self.to_string_base_(base, lowercase) + } } /// Equivalent to [`Arbi::to_string_base()`], but panics if the base is @@ -189,16 +336,17 @@ impl Arbi { /// let s = a.to_string_radix(10); /// assert_eq!(s, "123456789"); /// ``` + #[inline] pub fn to_string_radix(&self, radix: u32) -> String { let base: Base = match radix.try_into() { Err(_) => panic!("`radix` is not an integer in [2, 36]"), Ok(b) => b, }; - self.to_string_base(base) } } +/* TODO: clean up */ #[cfg(test)] pub(crate) mod tests { use super::*; @@ -360,4 +508,85 @@ pub(crate) mod tests { test_to_string_base(base); } } + + fn create_test_number(bits: usize) -> Arbi { + let mut num = Arbi::from(1); + for _ in 0..bits - 1 { + num <<= 1; + num += 1; + } + num + } + + #[test] + fn pow2_base_tests() { + use crate::util::test::random_arbi; + + let pow2_bases = [2, 4, 8, 16, 32]; + let bit_sizes = [ + 1, 2, 3, 4, 30, 31, 32, 33, 34, 62, 63, 64, 65, 66, 1024, 4096, + ]; + + for &bits in &bit_sizes { + for _ in 0..1000 { + let random_num = random_arbi(bits); + let neg_random = -random_num.clone(); + + for &base in &pow2_bases { + let base = Base::try_from(base).unwrap(); + + // Positive + let pow2_result = + random_num.to_string_base_pow2(base, true); + let general_result = random_num.to_string_base_(base, true); + assert_eq!(pow2_result, general_result,); + + // Negative + let neg_pow2_result = + neg_random.to_string_base_pow2(base, true); + let neg_general_result = + neg_random.to_string_base_(base, true); + assert_eq!(neg_pow2_result, neg_general_result,); + } + } + + let num = create_test_number(bits); + let neg_num = -num.clone(); + for &base in &pow2_bases { + let base = Base::try_from(base).unwrap(); + + // Positive + let pow2_result = num.to_string_base_pow2(base, true); + let general_result = num.to_string_base_(base, true); + assert_eq!(pow2_result, general_result,); + + // Negative + let neg_pow2_result = neg_num.to_string_base_pow2(base, true); + let neg_general_result = neg_num.to_string_base_(base, true); + assert_eq!(neg_pow2_result, neg_general_result,); + } + } + + let cases = [ + Arbi::zero(), + Arbi::from(1), + Arbi::from(-1), + create_test_number(1), + create_test_number(63), + create_test_number(64), + create_test_number(65), + create_test_number(127), + create_test_number(128), + create_test_number(129), + ]; + + for num in cases.iter() { + for &base in &pow2_bases { + let base = Base::try_from(base).unwrap(); + let pow2_result = num.to_string_base_pow2(base, true); + let general_result = num.to_string_base_(base, true); + assert_eq!(pow2_result, general_result,); + } + } + } } diff --git a/arbi/src/uints.rs b/arbi/src/uints.rs index e19cf99..6215a91 100644 --- a/arbi/src/uints.rs +++ b/arbi/src/uints.rs @@ -1,8 +1,10 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */ +use crate::BitCount; + #[allow(dead_code)] pub(crate) trait UnsignedUtilities: Sized { fn uaddc(r: &mut Self, a: Self, b: Self, carry: &mut u8); @@ -147,6 +149,15 @@ pub(crate) const fn div_ceil_usize(x: usize, y: usize) -> usize { } } +#[allow(dead_code)] +pub(crate) const fn div_ceil_bitcount(x: BitCount, y: BitCount) -> BitCount { + if x == 0 { + 0 + } else { + 1 + (x - 1) / y + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/arbi/src/util/mod.rs b/arbi/src/util/mod.rs index eaca99f..32edd2d 100644 --- a/arbi/src/util/mod.rs +++ b/arbi/src/util/mod.rs @@ -5,6 +5,8 @@ SPDX-License-Identifier: Apache-2.0 OR MIT pub(crate) mod arbi_like_array; pub(crate) mod max_digits; +pub(crate) mod mul; +pub(crate) mod radix_info; /// Utilities for testing. #[cfg(test)] pub(crate) mod test; diff --git a/arbi/src/util/mul.rs b/arbi/src/util/mul.rs new file mode 100644 index 0000000..2b88551 --- /dev/null +++ b/arbi/src/util/mul.rs @@ -0,0 +1,203 @@ +/* +Copyright 2025 Owain Davies +SPDX-License-Identifier: Apache-2.0 OR MIT +*/ + +macro_rules! define_mul2 { + // mul2(), mul2_halves() implement + // + // a * b -> (hi, lo) + // + // where a, b, hi, lo are all of type uint. + // + // They are equivalent in behavior. + // + // uint : primitive unsigned integer type (e.g. u128). + // uint_h : primitive unsigned integer type with size in bits half that + // of uint. + ($uint:ty, $uint_h:ty) => { + +const UINT_BITS: u32 = <$uint>::BITS; +const UINT_H_BITS: u32 = UINT_BITS >> 1; +const UINT_H_BASE: $uint = 1 as $uint << UINT_H_BITS; +const UINT_H_MAX: $uint = UINT_H_BASE - 1; // MASK for low part of UINT + +#[inline(always)] +#[allow(dead_code)] +const fn split_uint(x: $uint) -> ($uint, $uint) { + let hi = x >> UINT_H_BITS; + let lo = x & UINT_H_MAX; + (hi, lo) +} + +#[inline(always)] +const fn split_uint_halves(x: $uint) -> ($uint_h, $uint_h) { + let hi = (x >> UINT_H_BITS) as $uint_h; + let lo = (x & UINT_H_MAX) as $uint_h; + (hi, lo) +} + +#[inline] +#[allow(dead_code)] +pub(crate) const fn mul2(a: $uint, b: $uint) -> ($uint, $uint) { + let (a_hi, a_lo) = split_uint(a); + let (b_hi, b_lo) = split_uint(b); + + let mut ahbh = a_hi * b_hi; + let ahbl = a_hi * b_lo; + let mut albh = a_lo * b_hi; + let albl = a_lo * b_lo; + + let (albl_hi, albl_lo) = split_uint(albl); + + albh += albl_hi; + debug_assert!(albh >= albl_hi); + + let (albh, overflow) = albh.overflowing_add(ahbl); + if overflow { + ahbh = ahbh.wrapping_add(UINT_H_BASE); + } + + ( + ahbh.wrapping_add(albh >> UINT_H_BITS), + albl_lo.wrapping_add(albh << UINT_H_BITS), + ) +} + +#[inline] +pub(crate) const fn mul2_halves(a: $uint, b: $uint) -> ($uint, $uint) { + let (a_hi, a_lo) = split_uint_halves(a); + let (b_hi, b_lo) = split_uint_halves(b); + + let mut ahbh: $uint = (a_hi as $uint) * (b_hi as $uint); + let ahbl: $uint = (a_hi as $uint) * (b_lo as $uint); + let mut albh: $uint = (a_lo as $uint) * (b_hi as $uint); + let albl: $uint = (a_lo as $uint) * (b_lo as $uint); + + let (albl_hi, albl_lo) = split_uint_halves(albl); + + albh += albl_hi as $uint; + debug_assert!(albh >= albl_hi as $uint); + + let (albh, overflow) = albh.overflowing_add(ahbl); + if overflow { + ahbh = ahbh.wrapping_add(UINT_H_BASE); + } + + ( + ahbh.wrapping_add(albh >> UINT_H_BITS), + (albl_lo as $uint).wrapping_add(albh << UINT_H_BITS), + ) +} + + }; +} + +pub(crate) mod u128_impl { + define_mul2!(u128, u64); +} + +#[allow(dead_code)] +mod u64_impl { + define_mul2!(u64, u32); +} + +/* TODO: clean up and reduce repetition */ +#[cfg(test)] +mod tests { + use super::*; + use rand::Rng; + + #[test] + fn test_u128_implementations() { + let test_values = [ + 0u128, + 1u128, + 2u128, + 3u128, + u128::MAX, + 1u128 << 127, + (1u128 << 64) - 1, + 1u128 << 64, + ]; + + for a in test_values { + for b in test_values { + let r1 = u128_impl::mul2(a, b); + let r2 = u128_impl::mul2_halves(a, b); + assert_eq!(r1, r2, + "Results differ for u128 inputs {} and {}: basic: {:?} vs half: {:?}", + a, b, r1, r2); + } + } + + use crate::Arbi; + let mut rng = rand::thread_rng(); + for _ in 0..i16::MAX { + let a: u128 = rng.gen(); + let b: u128 = rng.gen(); + + let (hi, lo) = u128_impl::mul2(a, b); + + // Verify using Arbi integers + let res = Arbi::from(a) * Arbi::from(b); + let expected_lo = res.wrapping_to_u128(); // truncates + let expected_hi = (res >> 128u32).checked_to_u128().unwrap(); + + assert_eq!((hi, lo), (expected_hi, expected_lo), + "Mismatch for inputs {} and {}\nGot: ({}, {})\nExpected: ({}, {})", + a, b, hi, lo, expected_hi, expected_lo); + + let halves_res = u128_impl::mul2_halves(a, b); + assert_eq!(halves_res, (expected_hi, expected_lo), + "Halves implementation mismatch for inputs {} and {}\nGot: ({}, {})\nExpected: ({}, {})", + a, b, halves_res.0, halves_res.1, expected_hi, expected_lo); + } + } + + #[test] + fn test_u64_implementations() { + let test_values = [ + 0u64, + 1u64, + 2u64, + 3u64, + u64::MAX, + 1u64 << 63, + (1u64 << 32) - 1, + 1u64 << 32, + ]; + + for a in test_values { + for b in test_values { + let r1 = u64_impl::mul2(a, b); + let r2 = u64_impl::mul2_halves(a, b); + assert_eq!(r1, r2, + "Results differ for u64 inputs {} and {}: basic: {:?} vs half: {:?}", + a, b, r1, r2); + } + } + + let mut rng = rand::thread_rng(); + for _ in 0..i16::MAX { + let a: u64 = rng.gen(); + let b: u64 = rng.gen(); + + let (hi, lo) = u64_impl::mul2(a, b); + + // "native" u128 multiplication + let res = (a as u128) * (b as u128); + let expected_hi = (res >> 64) as u64; + let expected_lo = res as u64; + + assert_eq!((hi, lo), (expected_hi, expected_lo), + "Mismatch for inputs {} and {}\nGot: ({}, {})\nExpected: ({}, {})", + a, b, hi, lo, expected_hi, expected_lo); + + let halves_res = u64_impl::mul2_halves(a, b); + assert_eq!(halves_res, (expected_hi, expected_lo), + "Halves implementation mismatch for inputs {} and {}\nGot: ({}, {})\nExpected: ({}, {})", + a, b, halves_res.0, halves_res.1, expected_hi, expected_lo); + } + } +} diff --git a/arbi/src/util/radix_info.rs b/arbi/src/util/radix_info.rs new file mode 100644 index 0000000..6b019f5 --- /dev/null +++ b/arbi/src/util/radix_info.rs @@ -0,0 +1,191 @@ +/* +Copyright 2025 Owain Davies +SPDX-License-Identifier: Apache-2.0 OR MIT +*/ + +use crate::util::mul::u128_impl; +use crate::{Arbi, BitCount, Digit}; + +/* +Given an n-digit integer in base b, the largest integer representable is + +\\[ + b^{n} - 1 +\\] + +The minimum number of base b digits, n, required to represent any +m-digit base c integer is such that: + +\\[ + \begin{align} + b^{n} - 1 & \geq c^{m} - 1 \\\\ + b^{n} & \geq c^{m} \\\\ + \frac{b^{n}}{c^{m}} & \geq 1 \\\\ + \log(b^{n}) - \log(c^{m}) & \geq 0 \\\\ + n & \geq m \cdot \frac{\log(c)}{\log(b)} + \end{align} +\\] +*/ + +impl Arbi { + /* + SCALED_LOG2_DIV_LOG: + + ceil( (log(2) / log(base)) * 2^(Digit::BITS) ) + + SCALED_LOG2_DIV_LOG_64 (not used): + + ceil( (log(2) / log(base)) * 2^(64) ) + + All platforms currently use 32 as Digit::BITS, but in the future, we may + choose to use 64 for some platforms, in which case the values in the latter + array will be used. + + The values in both arrays were precomputed via a Python script using + Python `decimal.Decimal` objects. + + Technically, the results in the former array can be accurately computed + simply using f64s. However, f64s won't give us accurate values for the + latter array, due to precision issues. More complicated methods with the + help of u128s can be used for the former, but higher width integers (like + Arbi integers) would be needed for the latter. + */ + /// ceil( (log(2) / log(base)) * 2^(Digit::BITS) ) + const SCALED_LOG2_DIV_LOG: [Digit; 37] = [ + 0x00000000, 0x00000000, 0x00000000, 0xa1849cc2, 0x80000000, 0x6e40d1a5, + 0x6308c91c, 0x5b3064ec, 0x55555556, 0x50c24e61, 0x4d104d43, 0x4a002708, + 0x4768ce0e, 0x452e53e4, 0x433cfffc, 0x41867712, 0x40000000, 0x3ea16afe, + 0x3d64598e, 0x3c43c231, 0x3b3b9a43, 0x3a4898f1, 0x39680b14, 0x3897b2b8, + 0x37d5aed2, 0x372068d3, 0x3676867f, 0x35d6deec, 0x354071d7, 0x34b260c6, + 0x342be987, 0x33ac61ba, 0x33333334, 0x32bfd902, 0x3251dcf7, 0x31e8d5a0, + 0x3184648e, + ]; + + /// ceil( (log(2) / log(base)) * 2^(64) ) + #[allow(dead_code)] + const SCALED_LOG2_DIV_LOG_64: [u64; 37] = [ + 0x0000000000000000, + 0x0000000000000000, + 0x0000000000000000, + 0xa1849cc1a9a9e94f, + 0x8000000000000000, + 0x6e40d1a4143dcb95, + 0x6308c91b702a7cf5, + 0x5b3064eb3aa6d389, + 0x5555555555555556, + 0x50c24e60d4d4f4a8, + 0x4d104d427de7fbcd, + 0x4a00270775914e89, + 0x4768ce0d05818e13, + 0x452e53e365907bdb, + 0x433cfffb4b5aae56, + 0x41867711b4f85356, + 0x4000000000000000, + 0x3ea16afd58b10967, + 0x3d64598d154dc4df, + 0x3c43c23018bb5564, + 0x3b3b9a42873069c8, + 0x3a4898f06cf41aca, + 0x39680b13582e7c19, + 0x3897b2b751ae561b, + 0x37d5aed131f19c99, + 0x372068d20a1ee5cb, + 0x3676867e5d60de2a, + 0x35d6deeb388df870, + 0x354071d61c77fa2f, + 0x34b260c5671b18ad, + 0x342be986572b45cd, + 0x33ac61b998fbbdf3, + 0x3333333333333334, + 0x32bfd90114c12862, + 0x3251dcf6169e45f3, + 0x31e8d59f180dc631, + 0x3184648db8153e7b, + ]; + + /// Compute the number of base-`base` digits needed to represent a + /// nonnegative integer with `size_bits` bits. The true value or one higher + /// than it will be returned. + /// + /// # Panic + /// This function panics if `base` is not in the half-open interval + /// `(2, 36]`, or if `size_bits` is `0`. + pub(crate) const fn size_base_with_size_bits_maybe_over_by_one( + base: u32, + size_bits: BitCount, + ) -> BitCount { + assert!(base > 2 && base <= 36, "base must be in (2, 36]"); + assert!(size_bits != 0); + + let multiplicand = Self::SCALED_LOG2_DIV_LOG[base as usize] as BitCount; + + if let Some(product) = size_bits.checked_mul(multiplicand) { + product >> Digit::BITS + } else { + /* Multiplication overflowed, so we do it again, this time storing + * the product in a 256-bit unsigned integer (represented as two + * u128s). */ + let (hi, lo) = u128_impl::mul2_halves(size_bits, multiplicand); + + /* Now shift the product right by Digit::BITS */ + debug_assert!(hi >> Digit::BITS == 0); + (hi << (u128::BITS - Digit::BITS)) | (lo >> Digit::BITS) + } + .wrapping_add(1) + } +} + +#[cfg(test)] +mod tests_size_base_with_size_bits { + use super::*; + use crate::util::test::{get_seedable_rng, get_uniform_die, Distribution}; + use crate::{DDigit, Digit, QDigit}; + + fn size_base(mut x: u128, base: u32) -> BitCount { + if x == 0 { + return 0; + } + let mut count = 0; + while x > 0 { + x /= base as u128; + count += 1; + } + count + } + + fn size_bits(x: u128) -> BitCount { + (u128::BITS - x.leading_zeros()).into() + } + + #[test] + fn test_random_integers() { + let (mut rng, _) = get_seedable_rng(); + let d1 = get_uniform_die(0, Digit::MAX); + let d2 = get_uniform_die(Digit::MAX as DDigit + 1, DDigit::MAX); + let d3 = get_uniform_die(DDigit::MAX as QDigit + 1, QDigit::MAX); + + for _ in 0..i16::MAX { + let nums = [ + d1.sample(&mut rng) as u128, + d2.sample(&mut rng) as u128, + d3.sample(&mut rng) as u128, + ]; + + for num in nums { + let bits = size_bits(num); + for base in 3..=36 { + let estimated = + Arbi::size_base_with_size_bits_maybe_over_by_one( + base, bits, + ); + let actual = size_base(num, base); + assert!( + estimated == actual || estimated == actual + 1, + "Failed for num={}, base={}: estimated={}, actual={}, bits={}", + num, base, estimated, actual, bits + ); + } + } + } + } +} diff --git a/arbi/src/util/test.rs b/arbi/src/util/test.rs index 675f365..baf365c 100644 --- a/arbi/src/util/test.rs +++ b/arbi/src/util/test.rs @@ -1,5 +1,5 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */ @@ -80,6 +80,7 @@ pub(crate) fn ldexp(x: f64, exp: i32) -> f64 { pub(crate) mod float { pub(crate) const ZERO: f64 = 0.0; pub(crate) const MAX_INT: f64 = 9007199254740992.0; // 2^53 + pub(crate) const DBL_MAX_INT: u64 = 1 << 53; pub(crate) const MAX_INT_NEG: f64 = -9007199254740992.0; pub(crate) const LOWEST_DOUBLE: f64 = f64::MIN; pub(crate) const MAX_DOUBLE: f64 = f64::MAX; diff --git a/arbi/src/util/to_digits.rs b/arbi/src/util/to_digits.rs index 0b8b878..94fa836 100644 --- a/arbi/src/util/to_digits.rs +++ b/arbi/src/util/to_digits.rs @@ -1,5 +1,5 @@ /* -Copyright 2024 Owain Davies +Copyright 2024-2025 Owain Davies SPDX-License-Identifier: Apache-2.0 OR MIT */