From 8e0baecf5c57cba8d4dda9300ef929179162e7b6 Mon Sep 17 00:00:00 2001 From: Kent Overstreet Date: Thu, 14 Jul 2016 21:56:16 -0800 Subject: [PATCH] Drop some dependencies on BigDigit's size Before: test divide_0 ... bench: 1,011 ns/iter (+/- 184) test divide_1 ... bench: 18,535 ns/iter (+/- 770) test divide_2 ... bench: 990,467 ns/iter (+/- 91,980) test fac_to_string ... bench: 1,275 ns/iter (+/- 60) test factorial_100 ... bench: 6,453 ns/iter (+/- 101) test fib_100 ... bench: 1,142 ns/iter (+/- 99) test fib_1000 ... bench: 18,713 ns/iter (+/- 2,172) test fib_10000 ... bench: 1,197,965 ns/iter (+/- 21,178) test fib_to_string ... bench: 225 ns/iter (+/- 13) test from_str_radix_02 ... bench: 3,460 ns/iter (+/- 626) test from_str_radix_08 ... bench: 1,324 ns/iter (+/- 24) test from_str_radix_10 ... bench: 1,488 ns/iter (+/- 19) test from_str_radix_16 ... bench: 969 ns/iter (+/- 22) test from_str_radix_36 ... bench: 1,135 ns/iter (+/- 23) test hash ... bench: 102,126 ns/iter (+/- 1,016) test multiply_0 ... bench: 353 ns/iter (+/- 74) test multiply_1 ... bench: 31,006 ns/iter (+/- 679) test multiply_2 ... bench: 3,438,143 ns/iter (+/- 47,640) test pow_bench ... bench: 7,457,045 ns/iter (+/- 96,175) test shl ... bench: 5,627 ns/iter (+/- 121) test shr ... bench: 5,054 ns/iter (+/- 112) test to_str_radix_02 ... bench: 2,774 ns/iter (+/- 88) test to_str_radix_08 ... bench: 980 ns/iter (+/- 425) test to_str_radix_10 ... bench: 3,029 ns/iter (+/- 115) test to_str_radix_16 ... bench: 788 ns/iter (+/- 14) test to_str_radix_36 ... bench: 8,285 ns/iter (+/- 175) After: test divide_0 ... bench: 925 ns/iter (+/- 30) test divide_1 ... bench: 17,660 ns/iter (+/- 379) test divide_2 ... bench: 972,427 ns/iter (+/- 7,560) test fac_to_string ... bench: 1,260 ns/iter (+/- 36) test factorial_100 ... bench: 7,077 ns/iter (+/- 204) test fib_100 ... bench: 1,124 ns/iter (+/- 32) test fib_1000 ... bench: 18,475 ns/iter (+/- 166) test fib_10000 ... bench: 1,192,748 ns/iter (+/- 27,128) test fib_to_string ... bench: 228 ns/iter (+/- 10) test from_str_radix_02 ... bench: 3,379 ns/iter (+/- 74) test from_str_radix_08 ... bench: 1,355 ns/iter (+/- 24) test from_str_radix_10 ... bench: 1,470 ns/iter (+/- 20) test from_str_radix_16 ... bench: 958 ns/iter (+/- 239) test from_str_radix_36 ... bench: 1,137 ns/iter (+/- 19) test hash ... bench: 102,730 ns/iter (+/- 39,897) test multiply_0 ... bench: 351 ns/iter (+/- 15) test multiply_1 ... bench: 31,139 ns/iter (+/- 1,053) test multiply_2 ... bench: 3,464,509 ns/iter (+/- 124,235) test pow_bench ... bench: 7,448,428 ns/iter (+/- 326,903) test shl ... bench: 5,784 ns/iter (+/- 190) test shr ... bench: 4,820 ns/iter (+/- 63) test to_str_radix_02 ... bench: 2,757 ns/iter (+/- 33) test to_str_radix_08 ... bench: 989 ns/iter (+/- 67) test to_str_radix_10 ... bench: 3,045 ns/iter (+/- 70) test to_str_radix_16 ... bench: 787 ns/iter (+/- 24) test to_str_radix_36 ... bench: 8,257 ns/iter (+/- 117)/ --- bigint/src/lib.rs | 360 +++++++++++++++++++++++++++++----------------- 1 file changed, 227 insertions(+), 133 deletions(-) diff --git a/bigint/src/lib.rs b/bigint/src/lib.rs index d094a3b..a32bd7d 100644 --- a/bigint/src/lib.rs +++ b/bigint/src/lib.rs @@ -306,7 +306,8 @@ impl FromStr for BigUint { } } -// Read bitwise digits that evenly divide BigDigit +// Convert from a power of two radix (bits == ilog2(radix)) where bits evenly divides +// BigDigit::BITS fn from_bitwise_digits_le(v: &[u8], bits: usize) -> BigUint { debug_assert!(!v.is_empty() && bits <= 8 && big_digit::BITS % bits == 0); debug_assert!(v.iter().all(|&c| (c as BigDigit) < (1 << bits))); @@ -315,14 +316,15 @@ fn from_bitwise_digits_le(v: &[u8], bits: usize) -> BigUint { let data = v.chunks(digits_per_big_digit) .map(|chunk| { - chunk.iter().rev().fold(0u32, |acc, &c| (acc << bits) | c as BigDigit) + chunk.iter().rev().fold(0, |acc, &c| (acc << bits) | c as BigDigit) }) .collect(); BigUint::new(data) } -// Read bitwise digits that don't evenly divide BigDigit +// Convert from a power of two radix (bits == ilog2(radix)) where bits doesn't evenly divide +// BigDigit::BITS fn from_inexact_bitwise_digits_le(v: &[u8], bits: usize) -> BigUint { debug_assert!(!v.is_empty() && bits <= 8 && big_digit::BITS % bits != 0); debug_assert!(v.iter().all(|&c| (c as BigDigit) < (1 << bits))); @@ -331,15 +333,20 @@ fn from_inexact_bitwise_digits_le(v: &[u8], bits: usize) -> BigUint { let mut data = Vec::with_capacity(big_digits); let mut d = 0; - let mut dbits = 0; + let mut dbits = 0; // number of bits we currently have in d + + // walk v accumululating bits in d; whenever we accumulate big_digit::BITS in d, spit out a + // big_digit: for &c in v { - d |= (c as DoubleBigDigit) << dbits; + d |= (c as BigDigit) << dbits; dbits += bits; + if dbits >= big_digit::BITS { - let (hi, lo) = big_digit::from_doublebigdigit(d); - data.push(lo); - d = hi as DoubleBigDigit; + data.push(d); dbits -= big_digit::BITS; + // if dbits was > big_digit::BITS, we dropped some of the bits in c (they couldn't fit + // in d) - grab the bits we lost here: + d = (c as BigDigit) >> (bits - dbits); } } @@ -362,8 +369,7 @@ fn from_radix_digits_be(v: &[u8], radix: u32) -> BigUint { let mut data = Vec::with_capacity(big_digits as usize); let (base, power) = get_radix_base(radix); - debug_assert!(base < (1 << 32)); - let base = base as BigDigit; + let radix = radix as BigDigit; let r = v.len() % power; let i = if r == 0 { @@ -435,7 +441,7 @@ impl Num for BigUint { let res = if radix.is_power_of_two() { // Powers of two can use bitwise masks and shifting instead of multiplication - let bits = radix.trailing_zeros() as usize; + let bits = ilog2(radix); v.reverse(); if big_digit::BITS % bits == 0 { from_bitwise_digits_le(&v, bits) @@ -1349,6 +1355,46 @@ impl Integer for BigUint { } } +fn high_bits_to_u64(v: &BigUint) -> u64 { + match v.data.len() { + 0 => 0, + 1 => v.data[0] as u64, + _ => { + let mut bits = v.bits(); + let mut ret = 0u64; + let mut ret_bits = 0; + + for d in v.data.iter().rev() { + let digit_bits = (bits - 1) % big_digit::BITS + 1; + let bits_want = cmp::min(64 - ret_bits, digit_bits); + + if bits_want != 64 { + ret <<= bits_want; + } + ret |= *d as u64 >> (digit_bits - bits_want); + ret_bits += bits_want; + bits -= bits_want; + + if ret_bits == 64 { + break; + } + } + + ret + } + } +} + +/// Find last set bit +/// fls(0) == 0, fls(u32::MAX) == 32 +fn fls(v: T) -> usize { + std::mem::size_of::() * 8 - v.leading_zeros() as usize +} + +fn ilog2(v: T) -> usize { + fls(v) - 1 +} + impl ToPrimitive for BigUint { #[inline] fn to_i64(&self) -> Option { @@ -1362,76 +1408,53 @@ impl ToPrimitive for BigUint { }) } - // `DoubleBigDigit` size dependent #[inline] fn to_u64(&self) -> Option { - match self.data.len() { - 0 => Some(0), - 1 => Some(self.data[0] as u64), - 2 => Some(big_digit::to_doublebigdigit(self.data[1], self.data[0]) as u64), - _ => None, + let mut ret: u64 = 0; + let mut bits = 0; + + for i in self.data.iter() { + if bits >= 64 { + return None; + } + + ret += (*i as u64) << bits; + bits += big_digit::BITS; } + + Some(ret) } - // `DoubleBigDigit` size dependent #[inline] fn to_f32(&self) -> Option { - match self.data.len() { - 0 => Some(f32::zero()), - 1 => Some(self.data[0] as f32), - len => { - // this will prevent any overflow of exponent - if len > (f32::MAX_EXP as usize) / big_digit::BITS { - None - } else { - let exponent = (len - 2) * big_digit::BITS; - // we need 25 significant digits, 24 to be stored and 1 for rounding - // this gives at least 33 significant digits - let mantissa = big_digit::to_doublebigdigit(self.data[len - 1], - self.data[len - 2]); - // this cast handles rounding - let ret = (mantissa as f32) * 2.0.powi(exponent as i32); - if ret.is_infinite() { - None - } else { - Some(ret) - } - } + let mantissa = high_bits_to_u64(self); + let exponent = self.bits() - fls(mantissa); + + if exponent > f32::MAX_EXP as usize { + None + } else { + let ret = (mantissa as f32) * 2.0f32.powi(exponent as i32); + if ret.is_infinite() { + None + } else { + Some(ret) } } } - // `DoubleBigDigit` size dependent #[inline] fn to_f64(&self) -> Option { - match self.data.len() { - 0 => Some(f64::zero()), - 1 => Some(self.data[0] as f64), - 2 => Some(big_digit::to_doublebigdigit(self.data[1], self.data[0]) as f64), - len => { - // this will prevent any overflow of exponent - if len > (f64::MAX_EXP as usize) / big_digit::BITS { - None - } else { - let mut exponent = (len - 2) * big_digit::BITS; - let mut mantissa = big_digit::to_doublebigdigit(self.data[len - 1], - self.data[len - 2]); - // we need at least 54 significant bit digits, 53 to be stored and 1 for rounding - // so we take enough from the next BigDigit to make it up to 64 - let shift = mantissa.leading_zeros() as usize; - if shift > 0 { - mantissa <<= shift; - mantissa |= self.data[len - 3] as u64 >> (big_digit::BITS - shift); - exponent -= shift; - } - // this cast handles rounding - let ret = (mantissa as f64) * 2.0.powi(exponent as i32); - if ret.is_infinite() { - None - } else { - Some(ret) - } - } + let mantissa = high_bits_to_u64(self); + let exponent = self.bits() - fls(mantissa); + + if exponent > f64::MAX_EXP as usize { + None + } else { + let ret = (mantissa as f64) * 2.0f64.powi(exponent as i32); + if ret.is_infinite() { + None + } else { + Some(ret) } } } @@ -1484,14 +1507,17 @@ impl FromPrimitive for BigUint { } impl From for BigUint { - // `DoubleBigDigit` size dependent #[inline] - fn from(n: u64) -> Self { - match big_digit::from_doublebigdigit(n) { - (0, 0) => BigUint::zero(), - (0, n0) => BigUint { data: vec![n0] }, - (n1, n0) => BigUint { data: vec![n0, n1] }, + fn from(mut n: u64) -> Self { + let mut ret: BigUint = Zero::zero(); + + while n != 0 { + ret.data.push(n as BigDigit); + // don't overflow if BITS is 64: + n = (n >> 1) >> (big_digit::BITS - 1); } + + ret } } @@ -1591,28 +1617,36 @@ fn to_bitwise_digits_le(u: &BigUint, bits: usize) -> Vec { fn to_inexact_bitwise_digits_le(u: &BigUint, bits: usize) -> Vec { debug_assert!(!u.is_zero() && bits <= 8 && big_digit::BITS % bits != 0); - let last_i = u.data.len() - 1; - let mask: DoubleBigDigit = (1 << bits) - 1; + let mask: BigDigit = (1 << bits) - 1; let digits = (u.bits() + bits - 1) / bits; let mut res = Vec::with_capacity(digits); let mut r = 0; let mut rbits = 0; - for hi in u.data[..last_i].iter().cloned() { - r |= (hi as DoubleBigDigit) << rbits; + + for c in &u.data { + r |= *c << rbits; rbits += big_digit::BITS; while rbits >= bits { res.push((r & mask) as u8); r >>= bits; + + // r had more bits than it could fit - grab the bits we lost + if rbits > big_digit::BITS { + r = *c >> (big_digit::BITS - (rbits - bits)); + } + rbits -= bits; } } - r |= (u.data[last_i] as DoubleBigDigit) << rbits; - while r != 0 { - res.push((r & mask) as u8); - r >>= bits; + if rbits != 0 { + res.push(r as u8); + } + + while let Some(&0) = res.last() { + res.pop(); } res @@ -1629,8 +1663,7 @@ fn to_radix_digits_le(u: &BigUint, radix: u32) -> Vec { let mut digits = u.clone(); let (base, power) = get_radix_base(radix); - debug_assert!(base < (1 << 32)); - let base = base as BigDigit; + let radix = radix as BigDigit; while digits.data.len() > 1 { let (q, mut r) = div_rem_digit(digits, base); @@ -1659,7 +1692,7 @@ fn to_str_radix_reversed(u: &BigUint, radix: u32) -> Vec { let mut res = if radix.is_power_of_two() { // Powers of two can use bitwise masks and shifting instead of division - let bits = radix.trailing_zeros() as usize; + let bits = ilog2(radix); if big_digit::BITS % bits == 0 { to_bitwise_digits_le(u, bits) } else { @@ -1852,57 +1885,115 @@ impl serde::Deserialize for BigUint { } } -// `DoubleBigDigit` size dependent /// Returns the greatest power of the radix <= big_digit::BASE #[inline] -fn get_radix_base(radix: u32) -> (DoubleBigDigit, usize) { +fn get_radix_base(radix: u32) -> (BigDigit, usize) { + debug_assert!(2 <= radix && radix <= 36, "The radix must be within 2...36"); + debug_assert!(!radix.is_power_of_two()); + // To generate this table: - // let target = std::u32::max as u64 + 1; // for radix in 2u64..37 { - // let power = (target as f64).log(radix as f64) as u32; - // let base = radix.pow(power); + // let mut power = big_digit::BITS / fls(radix as u64); + // let mut base = radix.pow(power as u32); + // + // while let Some(b) = base.checked_mul(radix) { + // if b > big_digit::MAX { + // break; + // } + // base = b; + // power += 1; + // } + // // println!("({:10}, {:2}), // {:2}", base, power, radix); // } - const BASES: [(DoubleBigDigit, usize); 37] = [(0, 0), - (0, 0), - (4294967296, 32), // 2 - (3486784401, 20), // 3 - (4294967296, 16), // 4 - (1220703125, 13), // 5 - (2176782336, 12), // 6 - (1977326743, 11), // 7 - (1073741824, 10), // 8 - (3486784401, 10), // 9 - (1000000000, 9), // 10 - (2357947691, 9), // 11 - (429981696, 8), // 12 - (815730721, 8), // 13 - (1475789056, 8), // 14 - (2562890625, 8), // 15 - (4294967296, 8), // 16 - (410338673, 7), // 17 - (612220032, 7), // 18 - (893871739, 7), // 19 - (1280000000, 7), // 20 - (1801088541, 7), // 21 - (2494357888, 7), // 22 - (3404825447, 7), // 23 - (191102976, 6), // 24 - (244140625, 6), // 25 - (308915776, 6), // 26 - (387420489, 6), // 27 - (481890304, 6), // 28 - (594823321, 6), // 29 - (729000000, 6), // 30 - (887503681, 6), // 31 - (1073741824, 6), // 32 - (1291467969, 6), // 33 - (1544804416, 6), // 34 - (1838265625, 6), // 35 - (2176782336, 6) /* 36 */]; - assert!(2 <= radix && radix <= 36, "The radix must be within 2...36"); - BASES[radix as usize] + match big_digit::BITS { + 32 => { + const BASES: [(u32, usize); 37] = [(0, 0), (0, 0), + (0, 0), // 2 + (3486784401, 20),// 3 + (0, 0), // 4 + (1220703125, 13),// 5 + (2176782336, 12),// 6 + (1977326743, 11),// 7 + (0, 0), // 8 + (3486784401, 10),// 9 + (1000000000, 9), // 10 + (2357947691, 9), // 11 + (429981696, 8), // 12 + (815730721, 8), // 13 + (1475789056, 8), // 14 + (2562890625, 8), // 15 + (0, 0), // 16 + (410338673, 7), // 17 + (612220032, 7), // 18 + (893871739, 7), // 19 + (1280000000, 7), // 20 + (1801088541, 7), // 21 + (2494357888, 7), // 22 + (3404825447, 7), // 23 + (191102976, 6), // 24 + (244140625, 6), // 25 + (308915776, 6), // 26 + (387420489, 6), // 27 + (481890304, 6), // 28 + (594823321, 6), // 29 + (729000000, 6), // 30 + (887503681, 6), // 31 + (0, 0), // 32 + (1291467969, 6), // 33 + (1544804416, 6), // 34 + (1838265625, 6), // 35 + (2176782336, 6) // 36 + ]; + + let (base, power) = BASES[radix as usize]; + (base as BigDigit, power) + } + 64 => { + const BASES: [(u64, usize); 37] = [(0, 0), (0, 0), + (9223372036854775808, 63), // 2 + (12157665459056928801, 40), // 3 + (4611686018427387904, 31), // 4 + (7450580596923828125, 27), // 5 + (4738381338321616896, 24), // 6 + (3909821048582988049, 22), // 7 + (9223372036854775808, 21), // 8 + (12157665459056928801, 20), // 9 + (10000000000000000000, 19), // 10 + (5559917313492231481, 18), // 11 + (2218611106740436992, 17), // 12 + (8650415919381337933, 17), // 13 + (2177953337809371136, 16), // 14 + (6568408355712890625, 16), // 15 + (1152921504606846976, 15), // 16 + (2862423051509815793, 15), // 17 + (6746640616477458432, 15), // 18 + (15181127029874798299, 15), // 19 + (1638400000000000000, 14), // 20 + (3243919932521508681, 14), // 21 + (6221821273427820544, 14), // 22 + (11592836324538749809, 14), // 23 + (876488338465357824, 13), // 24 + (1490116119384765625, 13), // 25 + (2481152873203736576, 13), // 26 + (4052555153018976267, 13), // 27 + (6502111422497947648, 13), // 28 + (10260628712958602189, 13), // 29 + (15943230000000000000, 13), // 30 + (787662783788549761, 12), // 31 + (1152921504606846976, 12), // 32 + (1667889514952984961, 12), // 33 + (2386420683693101056, 12), // 34 + (3379220508056640625, 12), // 35 + (4738381338321616896, 12), // 36 + ]; + + let (base, power) = BASES[radix as usize]; + (base as BigDigit, power) + } + _ => panic!("Invalid bigdigit size") + } } /// A Sign is a `BigInt`'s composing element. @@ -3459,8 +3550,8 @@ mod biguint_tests { fn test_convert_i64() { fn check(b1: BigUint, i: i64) { let b2: BigUint = FromPrimitive::from_i64(i).unwrap(); - assert!(b1 == b2); - assert!(b1.to_i64().unwrap() == i); + assert_eq!(b1, b2); + assert_eq!(b1.to_i64().unwrap(), i); } check(Zero::zero(), 0); @@ -3484,8 +3575,8 @@ mod biguint_tests { fn test_convert_u64() { fn check(b1: BigUint, u: u64) { let b2: BigUint = FromPrimitive::from_u64(u).unwrap(); - assert!(b1 == b2); - assert!(b1.to_u64().unwrap() == u); + assert_eq!(b1, b2); + assert_eq!(b1.to_u64().unwrap(), u); } check(Zero::zero(), 0); @@ -3976,6 +4067,7 @@ mod biguint_tests { format!("2{}1", repeat("0").take(bits / 2 - 1).collect::())), (10, match bits { + 64 => "36893488147419103233".to_string(), 32 => "8589934593".to_string(), 16 => "131073".to_string(), _ => panic!(), @@ -3993,12 +4085,14 @@ mod biguint_tests { repeat("0").take(bits / 2 - 1).collect::())), (8, match bits { + 64 => "14000000000000000000004000000000000000000001".to_string(), 32 => "6000000000100000000001".to_string(), 16 => "140000400001".to_string(), _ => panic!(), }), (10, match bits { + 64 => "1020847100762815390427017310442723737601".to_string(), 32 => "55340232229718589441".to_string(), 16 => "12885032961".to_string(), _ => panic!(),