From 1f2590656b80dbf33b38b64c81a61f1a423ab5b2 Mon Sep 17 00:00:00 2001 From: Manca Bizjak Date: Wed, 11 Jul 2018 12:14:46 +0200 Subject: [PATCH 1/4] Implement Roots for BigInt and BigUint This commit implements num-integer::Roots trait for BigInt and BigUint types, and also adds sqrt, cbrt, nth_root as inherent methods to allow access to them without importing Roots trait. For each type tests were added as submodules in the roots test module. Signed-off-by: Manca Bizjak --- Cargo.toml | 2 +- benches/bigint.rs | 32 +++++++++++++++++- src/bigint.rs | 30 ++++++++++++++++- src/biguint.rs | 69 ++++++++++++++++++++++++++++++++++++-- tests/roots.rs | 84 +++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 212 insertions(+), 5 deletions(-) create mode 100644 tests/roots.rs diff --git a/Cargo.toml b/Cargo.toml index df3a2233..3fe47f2d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,7 +31,7 @@ name = "shootout-pidigits" [dependencies] [dependencies.num-integer] -version = "0.1.38" +version = "0.1.39" default-features = false [dependencies.num-traits] diff --git a/benches/bigint.rs b/benches/bigint.rs index 7a763c2d..e731ce1a 100644 --- a/benches/bigint.rs +++ b/benches/bigint.rs @@ -4,13 +4,14 @@ extern crate test; extern crate num_bigint; extern crate num_traits; +extern crate num_integer; extern crate rand; use std::mem::replace; use test::Bencher; use num_bigint::{BigInt, BigUint, RandBigInt}; use num_traits::{Zero, One, FromPrimitive, Num}; -use rand::{SeedableRng, StdRng}; +use rand::{SeedableRng, StdRng, Rng}; fn get_rng() -> StdRng { let mut seed = [0; 32]; @@ -342,3 +343,32 @@ fn modpow_even(b: &mut Bencher) { b.iter(|| base.modpow(&e, &m)); } + +#[bench] +fn roots_sqrt(b: &mut Bencher) { + let mut rng = get_rng(); + let x = rng.gen_biguint(2048); + + b.iter(|| x.sqrt()); +} + +#[bench] +fn roots_cbrt(b: &mut Bencher) { + let mut rng = get_rng(); + let x = rng.gen_biguint(2048); + + b.iter(|| x.cbrt()); +} + +#[bench] +fn roots_nth(b: &mut Bencher) { + let mut rng = get_rng(); + let x = rng.gen_biguint(2048); + // Although n is u32, here we limit it to the set of u8 values since it + // hugely impacts the performance of nth_root due to exponentiation to + // the power of n-1. Using very large values for n is also not very realistic, + // and any n > x's bit size produces 1 as a result anyway. + let n: u8 = rng.gen(); + + b.iter(|| { x.nth_root(n as u32) }); +} diff --git a/src/bigint.rs b/src/bigint.rs index 3c8d2962..beeae2c0 100644 --- a/src/bigint.rs +++ b/src/bigint.rs @@ -16,7 +16,7 @@ use std::iter::{Product, Sum}; #[cfg(feature = "serde")] use serde; -use integer::Integer; +use integer::{Integer, Roots}; use traits::{ToPrimitive, FromPrimitive, Num, CheckedAdd, CheckedSub, CheckedMul, CheckedDiv, Signed, Zero, One}; @@ -1802,6 +1802,15 @@ impl Integer for BigInt { } } +impl Roots for BigInt { + fn nth_root(&self, n: u32) -> Self { + assert!(!(self.is_negative() && n.is_even()), + "n-th root is undefined for number (n={})", n); + + BigInt::from_biguint(self.sign, self.data.nth_root(n)) + } +} + impl ToPrimitive for BigInt { #[inline] fn to_i64(&self) -> Option { @@ -2538,6 +2547,25 @@ impl BigInt { }; BigInt::from_biguint(sign, mag) } + + /// Returns the truncated principal square root of `self` -- + /// see [Roots::sqrt](Roots::sqrt). + // struct.BigInt.html#trait.Roots + pub fn sqrt(&self) -> Self { + Roots::sqrt(self) + } + + /// Returns the truncated principal cube root of `self` -- + /// see [Roots::cbrt](Roots::cbrt). + pub fn cbrt(&self) -> Self { + Roots::cbrt(self) + } + + /// Returns the truncated principal `n`th root of `self` -- + /// See [Roots::nth_root](Roots::nth_root). + pub fn nth_root(&self, n: u32) -> Self { + Roots::nth_root(self, n) + } } impl_sum_iter_type!(BigInt); diff --git a/src/biguint.rs b/src/biguint.rs index 5d4aaf89..6e7c991c 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -17,9 +17,9 @@ use std::ascii::AsciiExt; #[cfg(feature = "serde")] use serde; -use integer::Integer; +use integer::{Integer, Roots}; use traits::{ToPrimitive, FromPrimitive, Float, Num, Unsigned, CheckedAdd, CheckedSub, CheckedMul, - CheckedDiv, Zero, One}; + CheckedDiv, Zero, One, pow}; use big_digit::{self, BigDigit, DoubleBigDigit}; @@ -1026,6 +1026,52 @@ impl Integer for BigUint { } } +impl Roots for BigUint { + fn nth_root(&self, n: u32) -> Self { + assert!(n > 0, "n must be at least 1"); + + let one = BigUint::one(); + + // Trivial cases + if self.is_zero() { + return BigUint::zero(); + } + + if self.is_one() { + return one; + } + + let n = n as usize; + let n_min_1 = (n as usize) - 1; + + // Newton's method to compute the nth root of an integer. + // + // Reference: + // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.14 + // + // Set initial guess to something definitely >= floor(nth_root of self) + // but as low as possible to speed up convergence. + let bit_len = self.len() * big_digit::BITS; + let guess = one << (bit_len/n + 1); + + let mut u = guess; + let mut s: BigUint; + + loop { + s = u; + let q = self / pow(s.clone(), n_min_1); + let t: BigUint = n_min_1 * &s + q; + + // Compute the candidate value for next iteration + u = t / n; + + if u >= s { break; } + } + + s + } +} + fn high_bits_to_u64(v: &BigUint) -> u64 { match v.data.len() { 0 => 0, @@ -1749,6 +1795,25 @@ impl BigUint { } acc } + + /// Returns the truncated principal square root of `self` -- + /// see [Roots::sqrt](Roots::sqrt). + // struct.BigInt.html#trait.Roots + pub fn sqrt(&self) -> Self { + Roots::sqrt(self) + } + + /// Returns the truncated principal cube root of `self` -- + /// see [Roots::cbrt](Roots::cbrt). + pub fn cbrt(&self) -> Self { + Roots::cbrt(self) + } + + /// Returns the truncated principal `n`th root of `self` -- + /// See [Roots::nth_root](Roots::nth_root). + pub fn nth_root(&self, n: u32) -> Self { + Roots::nth_root(self, n) + } } /// Returns the number of least-significant bits that are zero, diff --git a/tests/roots.rs b/tests/roots.rs new file mode 100644 index 00000000..4ceb5491 --- /dev/null +++ b/tests/roots.rs @@ -0,0 +1,84 @@ +extern crate num_bigint; +extern crate num_integer; +extern crate num_traits; + +mod biguint { + use num_bigint::BigUint; + use num_traits::FromPrimitive; + use std::str::FromStr; + + fn check(x: i32, n: u32, expected: i32) { + let big_x: BigUint = FromPrimitive::from_i32(x).unwrap(); + let big_expected: BigUint = FromPrimitive::from_i32(expected).unwrap(); + + assert_eq!(big_x.nth_root(n), big_expected); + } + + #[test] + fn test_sqrt() { + check(99, 2, 9); + check(100, 2, 10); + check(120, 2, 10); + } + + #[test] + fn test_cbrt() { + check(8, 3, 2); + check(26, 3, 2); + } + + #[test] + fn test_nth_root() { + check(0, 1, 0); + check(10, 1, 10); + check(100, 4, 3); + } + + #[test] + #[should_panic] + fn test_nth_root_n_is_zero() { + check(4, 0, 0); + } + + #[test] + fn test_nth_root_big() { + let x: BigUint = FromStr::from_str("123_456_789").unwrap(); + let expected : BigUint = FromPrimitive::from_i32(6).unwrap(); + + assert_eq!(x.nth_root(10), expected); + } +} + +mod bigint { + use num_bigint::BigInt; + use num_traits::FromPrimitive; + + fn check(x: i32, n: u32, expected: i32) { + let big_x: BigInt = FromPrimitive::from_i32(x).unwrap(); + let big_expected: BigInt = FromPrimitive::from_i32(expected).unwrap(); + + assert_eq!(big_x.nth_root(n), big_expected); + } + + #[test] + fn test_nth_root() { + check(-100, 3, -4); + } + + #[test] + #[should_panic] + fn test_nth_root_x_neg_n_even() { + check(-100, 4, 0); + } + + #[test] + #[should_panic] + fn test_sqrt_x_neg() { + check(-4, 2, -2); + } + + #[test] + fn test_cbrt() { + check(-8, 3, -2); + } +} From 2b473e940365e96833a1956344a29f9cea113e30 Mon Sep 17 00:00:00 2001 From: Manca Bizjak Date: Fri, 13 Jul 2018 12:52:40 +0200 Subject: [PATCH 2/4] Implement optimized sqrt, cbrt methods This commit overrides default implementations of Roots::sqrt and Roots::cbrt for BigInt and BigUint with optimized ones. It also improves tests and resolves minor inconsistencies. Signed-off-by: Manca Bizjak --- benches/bigint.rs | 11 ++---- src/bigint.rs | 12 ++++++- src/biguint.rs | 88 +++++++++++++++++++++++++++++++++++------------ tests/roots.rs | 74 ++++++++++++++++++++++++--------------- 4 files changed, 127 insertions(+), 58 deletions(-) diff --git a/benches/bigint.rs b/benches/bigint.rs index e731ce1a..6893f137 100644 --- a/benches/bigint.rs +++ b/benches/bigint.rs @@ -11,7 +11,7 @@ use std::mem::replace; use test::Bencher; use num_bigint::{BigInt, BigUint, RandBigInt}; use num_traits::{Zero, One, FromPrimitive, Num}; -use rand::{SeedableRng, StdRng, Rng}; +use rand::{SeedableRng, StdRng}; fn get_rng() -> StdRng { let mut seed = [0; 32]; @@ -361,14 +361,9 @@ fn roots_cbrt(b: &mut Bencher) { } #[bench] -fn roots_nth(b: &mut Bencher) { +fn roots_nth_100(b: &mut Bencher) { let mut rng = get_rng(); let x = rng.gen_biguint(2048); - // Although n is u32, here we limit it to the set of u8 values since it - // hugely impacts the performance of nth_root due to exponentiation to - // the power of n-1. Using very large values for n is also not very realistic, - // and any n > x's bit size produces 1 as a result anyway. - let n: u8 = rng.gen(); - b.iter(|| { x.nth_root(n as u32) }); + b.iter(|| x.nth_root(100)); } diff --git a/src/bigint.rs b/src/bigint.rs index beeae2c0..c5195fe7 100644 --- a/src/bigint.rs +++ b/src/bigint.rs @@ -1805,10 +1805,20 @@ impl Integer for BigInt { impl Roots for BigInt { fn nth_root(&self, n: u32) -> Self { assert!(!(self.is_negative() && n.is_even()), - "n-th root is undefined for number (n={})", n); + "root of degree {} is imaginary", n); BigInt::from_biguint(self.sign, self.data.nth_root(n)) } + + fn sqrt(&self) -> Self { + assert!(!self.is_negative(), "square root is imaginary"); + + BigInt::from_biguint(self.sign, self.data.sqrt()) + } + + fn cbrt(&self) -> Self { + BigInt::from_biguint(self.sign, self.data.cbrt()) + } } impl ToPrimitive for BigInt { diff --git a/src/biguint.rs b/src/biguint.rs index 6e7c991c..87ed7531 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -1027,32 +1027,30 @@ impl Integer for BigUint { } impl Roots for BigUint { - fn nth_root(&self, n: u32) -> Self { - assert!(n > 0, "n must be at least 1"); + // nth_root, sqrt and cbrt use Newton's method to compute + // principal root of a given degree for a given integer. - let one = BigUint::one(); + // Reference: + // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.14 + fn nth_root(&self, n: u32) -> Self { + assert!(n > 0, "root degree n must be at least 1"); - // Trivial cases - if self.is_zero() { - return BigUint::zero(); + if self.is_zero() || self.is_one() { + return self.clone() } - if self.is_one() { - return one; + match n { // Optimize for small n + 1 => return self.clone(), + 2 => return self.sqrt(), + 3 => return self.cbrt(), + _ => (), } let n = n as usize; - let n_min_1 = (n as usize) - 1; - - // Newton's method to compute the nth root of an integer. - // - // Reference: - // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.14 - // - // Set initial guess to something definitely >= floor(nth_root of self) - // but as low as possible to speed up convergence. + let n_min_1 = n - 1; + let bit_len = self.len() * big_digit::BITS; - let guess = one << (bit_len/n + 1); + let guess = BigUint::one() << (bit_len/n + 1); let mut u = guess; let mut s: BigUint; @@ -1062,7 +1060,6 @@ impl Roots for BigUint { let q = self / pow(s.clone(), n_min_1); let t: BigUint = n_min_1 * &s + q; - // Compute the candidate value for next iteration u = t / n; if u >= s { break; } @@ -1070,6 +1067,54 @@ impl Roots for BigUint { s } + + // Reference: + // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.13 + fn sqrt(&self) -> Self { + if self.is_zero() || self.is_one() { + return self.clone() + } + + let bit_len = self.len() * big_digit::BITS; + let guess = BigUint::one() << (bit_len/2 + 1); + + let mut u = guess; + let mut s: BigUint; + + loop { + s = u; + let q = self / &s; + let t: BigUint = &s + q; + u = t >> 1; + + if u >= s { break; } + } + + s + } + + fn cbrt(&self) -> Self { + if self.is_zero() || self.is_one() { + return self.clone() + } + + let bit_len = self.len() * big_digit::BITS; + let guess = BigUint::one() << (bit_len/3 + 1); + + let mut u = guess; + let mut s: BigUint; + + loop { + s = u; + let q = self / (&s * &s); + let t: BigUint = (&s << 1) + q; + u = t / 3u32; + + if u >= s { break; } + } + + s + } } fn high_bits_to_u64(v: &BigUint) -> u64 { @@ -1797,8 +1842,7 @@ impl BigUint { } /// Returns the truncated principal square root of `self` -- - /// see [Roots::sqrt](Roots::sqrt). - // struct.BigInt.html#trait.Roots + /// see [Roots::sqrt](Roots::sqrt) pub fn sqrt(&self) -> Self { Roots::sqrt(self) } @@ -1810,7 +1854,7 @@ impl BigUint { } /// Returns the truncated principal `n`th root of `self` -- - /// See [Roots::nth_root](Roots::nth_root). + /// see [Roots::nth_root](Roots::nth_root). pub fn nth_root(&self, n: u32) -> Self { Roots::nth_root(self, n) } diff --git a/tests/roots.rs b/tests/roots.rs index 4ceb5491..58838b76 100644 --- a/tests/roots.rs +++ b/tests/roots.rs @@ -4,46 +4,53 @@ extern crate num_traits; mod biguint { use num_bigint::BigUint; - use num_traits::FromPrimitive; + use num_traits::pow; use std::str::FromStr; - fn check(x: i32, n: u32, expected: i32) { - let big_x: BigUint = FromPrimitive::from_i32(x).unwrap(); - let big_expected: BigUint = FromPrimitive::from_i32(expected).unwrap(); + fn check(x: u64, n: u32) { + let big_x = BigUint::from(x); + let res = big_x.nth_root(n); - assert_eq!(big_x.nth_root(n), big_expected); + if n == 2 { + assert_eq!(&res, &big_x.sqrt()) + } else if n == 3 { + assert_eq!(&res, &big_x.cbrt()) + } + + assert!(pow(res.clone(), n as usize) <= big_x); + assert!(pow(res.clone() + 1u32, n as usize) > big_x); } #[test] fn test_sqrt() { - check(99, 2, 9); - check(100, 2, 10); - check(120, 2, 10); + check(99, 2); + check(100, 2); + check(120, 2); } #[test] fn test_cbrt() { - check(8, 3, 2); - check(26, 3, 2); + check(8, 3); + check(26, 3); } #[test] fn test_nth_root() { - check(0, 1, 0); - check(10, 1, 10); - check(100, 4, 3); + check(0, 1); + check(10, 1); + check(100, 4); } #[test] #[should_panic] fn test_nth_root_n_is_zero() { - check(4, 0, 0); + check(4, 0); } #[test] fn test_nth_root_big() { - let x: BigUint = FromStr::from_str("123_456_789").unwrap(); - let expected : BigUint = FromPrimitive::from_i32(6).unwrap(); + let x = BigUint::from_str("123_456_789").unwrap(); + let expected = BigUint::from(6u32); assert_eq!(x.nth_root(10), expected); } @@ -51,34 +58,47 @@ mod biguint { mod bigint { use num_bigint::BigInt; - use num_traits::FromPrimitive; - - fn check(x: i32, n: u32, expected: i32) { - let big_x: BigInt = FromPrimitive::from_i32(x).unwrap(); - let big_expected: BigInt = FromPrimitive::from_i32(expected).unwrap(); - - assert_eq!(big_x.nth_root(n), big_expected); + use num_traits::{Signed, pow}; + + fn check(x: i64, n: u32) { + let big_x = BigInt::from(x); + let res = big_x.nth_root(n); + + if n == 2 { + assert_eq!(&res, &big_x.sqrt()) + } else if n == 3 { + assert_eq!(&res, &big_x.cbrt()) + } + + if big_x.is_negative() { + assert!(pow(res.clone() - 1u32, n as usize) < big_x); + assert!(pow(res.clone(), n as usize) >= big_x); + } else { + assert!(pow(res.clone(), n as usize) <= big_x); + assert!(pow(res.clone() + 1u32, n as usize) > big_x); + } } #[test] fn test_nth_root() { - check(-100, 3, -4); + check(-100, 3); } #[test] #[should_panic] fn test_nth_root_x_neg_n_even() { - check(-100, 4, 0); + check(-100, 4); } #[test] #[should_panic] fn test_sqrt_x_neg() { - check(-4, 2, -2); + check(-4, 2); } #[test] fn test_cbrt() { - check(-8, 3, -2); + check(8, 3); + check(-8, 3); } } From bbe7b188196017f3de033f72e29ea1d3dbe12734 Mon Sep 17 00:00:00 2001 From: Manca Bizjak Date: Tue, 17 Jul 2018 15:06:55 +0200 Subject: [PATCH 3/4] Remove stale doc path, use self.bits() Signed-off-by: Manca Bizjak --- src/bigint.rs | 1 - src/biguint.rs | 9 +++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/bigint.rs b/src/bigint.rs index c5195fe7..9693e473 100644 --- a/src/bigint.rs +++ b/src/bigint.rs @@ -2560,7 +2560,6 @@ impl BigInt { /// Returns the truncated principal square root of `self` -- /// see [Roots::sqrt](Roots::sqrt). - // struct.BigInt.html#trait.Roots pub fn sqrt(&self) -> Self { Roots::sqrt(self) } diff --git a/src/biguint.rs b/src/biguint.rs index 87ed7531..5411f798 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -1049,8 +1049,7 @@ impl Roots for BigUint { let n = n as usize; let n_min_1 = n - 1; - let bit_len = self.len() * big_digit::BITS; - let guess = BigUint::one() << (bit_len/n + 1); + let guess = BigUint::one() << (self.bits()/n + 1); let mut u = guess; let mut s: BigUint; @@ -1075,8 +1074,7 @@ impl Roots for BigUint { return self.clone() } - let bit_len = self.len() * big_digit::BITS; - let guess = BigUint::one() << (bit_len/2 + 1); + let guess = BigUint::one() << (self.bits()/2 + 1); let mut u = guess; let mut s: BigUint; @@ -1098,8 +1096,7 @@ impl Roots for BigUint { return self.clone() } - let bit_len = self.len() * big_digit::BITS; - let guess = BigUint::one() << (bit_len/3 + 1); + let guess = BigUint::one() << (self.bits()/3 + 1); let mut u = guess; let mut s: BigUint; From 1d45ca9a6f9de9d2f517bb444de71a0e09a9f42c Mon Sep 17 00:00:00 2001 From: Manca Bizjak Date: Wed, 18 Jul 2018 08:22:51 +0200 Subject: [PATCH 4/4] Update links to Roots trait methods Signed-off-by: Manca Bizjak --- src/bigint.rs | 6 +++--- src/biguint.rs | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/bigint.rs b/src/bigint.rs index 9693e473..93bb6b26 100644 --- a/src/bigint.rs +++ b/src/bigint.rs @@ -2559,19 +2559,19 @@ impl BigInt { } /// Returns the truncated principal square root of `self` -- - /// see [Roots::sqrt](Roots::sqrt). + /// see [Roots::sqrt](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#method.sqrt). pub fn sqrt(&self) -> Self { Roots::sqrt(self) } /// Returns the truncated principal cube root of `self` -- - /// see [Roots::cbrt](Roots::cbrt). + /// see [Roots::cbrt](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#method.cbrt). pub fn cbrt(&self) -> Self { Roots::cbrt(self) } /// Returns the truncated principal `n`th root of `self` -- - /// See [Roots::nth_root](Roots::nth_root). + /// See [Roots::nth_root](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#tymethod.nth_root). pub fn nth_root(&self, n: u32) -> Self { Roots::nth_root(self, n) } diff --git a/src/biguint.rs b/src/biguint.rs index 5411f798..e7a3ce19 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -1839,19 +1839,19 @@ impl BigUint { } /// Returns the truncated principal square root of `self` -- - /// see [Roots::sqrt](Roots::sqrt) + /// see [Roots::sqrt](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#method.sqrt) pub fn sqrt(&self) -> Self { Roots::sqrt(self) } /// Returns the truncated principal cube root of `self` -- - /// see [Roots::cbrt](Roots::cbrt). + /// see [Roots::cbrt](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#method.cbrt). pub fn cbrt(&self) -> Self { Roots::cbrt(self) } /// Returns the truncated principal `n`th root of `self` -- - /// see [Roots::nth_root](Roots::nth_root). + /// see [Roots::nth_root](https://docs.rs/num-integer/0.1/num_integer/trait.Roots.html#tymethod.nth_root). pub fn nth_root(&self, n: u32) -> Self { Roots::nth_root(self, n) }