From 2f6eb3190347476abe6ffb7195bd9365891bdce3 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Sat, 11 Mar 2017 16:41:56 +0100 Subject: [PATCH 1/3] Implement binomial The implementation avoids overflows and fractions. --- integer/src/lib.rs | 75 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) diff --git a/integer/src/lib.rs b/integer/src/lib.rs index 00cc3ba..923c5b5 100644 --- a/integer/src/lib.rs +++ b/integer/src/lib.rs @@ -665,6 +665,36 @@ impl_integer_for_usize!(u32, test_integer_u32); impl_integer_for_usize!(u64, test_integer_u64); impl_integer_for_usize!(usize, test_integer_usize); +/// Calculate r * a / b, avoiding overflows and fractions. +fn multiply_and_divide(r: T, a: T, b: T) -> T { + // See http://blog.plover.com/math/choose-2.html for the idea. + let g = gcd(r.clone(), b.clone()); + (r/g.clone()) * (a / (b/g)) +} + +/// Calculate the binomial coefficient. +pub fn binomial(mut n: T, k: T) -> T { + // See http://blog.plover.com/math/choose.html for the idea. + if k > n { + return T::zero(); + } + let two = T::one() + T::one(); + if k > n.clone()/two { + return binomial(n.clone(), n - k); + } + let mut r = T::one(); + let mut d = T::one(); + loop { + if d > k { + break; + } + r = multiply_and_divide(r, n.clone(), d.clone()); + n = n - T::one(); + d = d + T::one(); + } + r +} + #[test] fn test_lcm_overflow() { macro_rules! check { @@ -692,3 +722,48 @@ fn test_lcm_overflow() { check!(i64, 0x4000_0000_0000_0000, 0x04, 0x4000_0000_0000_0000); check!(u64, 0x8000_0000_0000_0000, 0x02, 0x8000_0000_0000_0000); } + +#[test] +fn test_binomial() { + macro_rules! check { + ($t:ty, $x:expr, $y:expr, $r:expr) => { { + let x: $t = $x; + let y: $t = $y; + let expected: $t = $r; + assert_eq!(binomial(x, y), expected); + } } + } + check!(u8, 9, 4, 126); + check!(u8, 0, 0, 1); + + check!(i8, 9, 4, 126); + check!(i8, 0, 0, 1); + + check!(u16, 100, 2, 4950); + check!(u16, 14, 4, 1001); + check!(u16, 0, 0, 1); + + check!(i16, 100, 2, 4950); + check!(i16, 14, 4, 1001); + check!(i16, 0, 0, 1); + + check!(u32, 100, 2, 4950); + check!(u32, 35, 11, 417225900); + check!(u32, 14, 4, 1001); + check!(u32, 0, 0, 1); + + check!(i32, 100, 2, 4950); + check!(i32, 35, 11, 417225900); + check!(i32, 14, 4, 1001); + check!(i32, 0, 0, 1); + + check!(u64, 100, 2, 4950); + check!(u64, 35, 11, 417225900); + check!(u64, 14, 4, 1001); + check!(u64, 0, 0, 1); + + check!(i64, 100, 2, 4950); + check!(i64, 35, 11, 417225900); + check!(i64, 14, 4, 1001); + check!(i64, 0, 0, 1); +} From 07df43b0341ec296e3539e3effc5271ceabd784d Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Sun, 12 Mar 2017 10:08:26 +0100 Subject: [PATCH 2/3] Address binomial review feedback * Use substraction instead of division in comparison. * More tests. * Add comment. --- integer/src/lib.rs | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/integer/src/lib.rs b/integer/src/lib.rs index 923c5b5..a9c5895 100644 --- a/integer/src/lib.rs +++ b/integer/src/lib.rs @@ -668,6 +668,9 @@ impl_integer_for_usize!(usize, test_integer_usize); /// Calculate r * a / b, avoiding overflows and fractions. fn multiply_and_divide(r: T, a: T, b: T) -> T { // See http://blog.plover.com/math/choose-2.html for the idea. + // + // This depends on the fact that `b` must evenly divide `r*a`, as that's + // what lets you know that `b/gcd(r, b)` divides `a` evenly. let g = gcd(r.clone(), b.clone()); (r/g.clone()) * (a / (b/g)) } @@ -678,8 +681,7 @@ pub fn binomial(mut n: T, k: T) -> T { if k > n { return T::zero(); } - let two = T::one() + T::one(); - if k > n.clone()/two { + if k > n.clone() - k.clone() { return binomial(n.clone(), n - k); } let mut r = T::one(); @@ -731,39 +733,50 @@ fn test_binomial() { let y: $t = $y; let expected: $t = $r; assert_eq!(binomial(x, y), expected); + if y <= x { + assert_eq!(binomial(x, x - y), expected); + } } } } check!(u8, 9, 4, 126); check!(u8, 0, 0, 1); + check!(u8, 2, 3, 0); check!(i8, 9, 4, 126); check!(i8, 0, 0, 1); + check!(i8, 2, 3, 0); check!(u16, 100, 2, 4950); check!(u16, 14, 4, 1001); check!(u16, 0, 0, 1); + check!(u16, 2, 3, 0); check!(i16, 100, 2, 4950); check!(i16, 14, 4, 1001); check!(i16, 0, 0, 1); + check!(i16, 2, 3, 0); check!(u32, 100, 2, 4950); check!(u32, 35, 11, 417225900); check!(u32, 14, 4, 1001); check!(u32, 0, 0, 1); + check!(u32, 2, 3, 0); check!(i32, 100, 2, 4950); check!(i32, 35, 11, 417225900); check!(i32, 14, 4, 1001); check!(i32, 0, 0, 1); + check!(i32, 2, 3, 0); check!(u64, 100, 2, 4950); check!(u64, 35, 11, 417225900); check!(u64, 14, 4, 1001); check!(u64, 0, 0, 1); + check!(u64, 2, 3, 0); check!(i64, 100, 2, 4950); check!(i64, 35, 11, 417225900); check!(i64, 14, 4, 1001); check!(i64, 0, 0, 1); + check!(i64, 2, 3, 0); } From 40cbe1ff200ab9b19ea39ae31a63d1cb018cf2ea Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Sun, 12 Mar 2017 10:23:27 +0100 Subject: [PATCH 3/3] Use correct formula --- integer/src/lib.rs | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/integer/src/lib.rs b/integer/src/lib.rs index a9c5895..b84b7af 100644 --- a/integer/src/lib.rs +++ b/integer/src/lib.rs @@ -668,11 +668,8 @@ impl_integer_for_usize!(usize, test_integer_usize); /// Calculate r * a / b, avoiding overflows and fractions. fn multiply_and_divide(r: T, a: T, b: T) -> T { // See http://blog.plover.com/math/choose-2.html for the idea. - // - // This depends on the fact that `b` must evenly divide `r*a`, as that's - // what lets you know that `b/gcd(r, b)` divides `a` evenly. let g = gcd(r.clone(), b.clone()); - (r/g.clone()) * (a / (b/g)) + (r/g.clone() * a) / (b/g) } /// Calculate the binomial coefficient.