From 7310f51996250dc0ef4251bd961fef7f58bfe24a Mon Sep 17 00:00:00 2001 From: Patrick Norton Date: Sun, 7 Nov 2021 10:22:55 -0800 Subject: [PATCH 1/3] Added initial factorial implementation --- src/biguint.rs | 1 + src/biguint/factorial.rs | 101 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+) create mode 100644 src/biguint/factorial.rs diff --git a/src/biguint.rs b/src/biguint.rs index 623823c8..9d898ca4 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -20,6 +20,7 @@ mod subtraction; mod bits; mod convert; +mod factorial; mod iter; mod monty; mod power; diff --git a/src/biguint/factorial.rs b/src/biguint/factorial.rs new file mode 100644 index 00000000..e813114c --- /dev/null +++ b/src/biguint/factorial.rs @@ -0,0 +1,101 @@ +use crate::BigUint; +use num_traits::{One, ToPrimitive}; + +/// All factorials that fit into a u64 +const SMALL_FACTORIALS: [u64; 21] = [ + 1, // 0! = 1 + 1, // 1! = 1 + 2, // 2! = 2 + 6, // 3! = 6 + 24, // 4! = 24 + 120, // 5! = 120 + 720, // 6! = 720 + 5040, // 7! = 5040 + 40320, // 8! = 40320 + 362880, // 9! = 362880 + 3628800, // 10! = 3628800 + 39916800, // 11! = 39916800 + 479001600, // 12! = 479001600 + 6227020800, // 13! = 6227020800 + 87178291200, // 14! = 87178291200 + 1307674368000, // 15! = 1307674368000 + 20922789888000, // 16! = 20922789888000 + 355687428096000, // 17! = 355687428096000 + 6402373705728000, // 18! = 6402373705728000 + 121645100408832000, // 19! = 121645100408832000 + 2432902008176640000, // 20! = 2432902008176640000 +]; + +/// The odd factors of each factorial: Each number in this table is odd, and +/// self[i] << k == i! for some k +const ODD_FACTORIALS: [u64; 25] = [ + 1, // 0! = 1 << 0 + 1, // 1! = 1 << 0 + 1, // 2! = 1 << 1 + 3, // 3! = 3 << 1 + 3, // 4! = 3 << 3 + 15, // 5! = 15 << 3 + 45, // 6! = 45 << 4 + 315, // 7! = 315 << 4 + 315, // 8! = 315 << 7 + 2835, // 9! = 2835 << 7 + 14175, // 10! = 14175 << 8 + 155925, // 11! = 155925 << 8 + 467775, // 12! = 467775 << 10 + 6081075, // 13! = 6081075 << 10 + 42567525, // 14! = 42567525 << 11 + 638512875, // 15! = 638512875 << 11 + 638512875, // 16! = 638512875 << 15 + 10854718875, // 17! = 10854718875 << 15 + 97692469875, // 18! = 97692469875 << 16 + 1856156927625, // 19! = 1856156927625 << 16 + 9280784638125, // 20! = 9280784638125 << 18 + 194896477400625, // 21! = 194896477400625 << 18 + 2143861251406875, // 22! = 2143861251406875 << 19 + 49308808782358125, // 23! = 49308808782358125 << 19 + 147926426347074375, // 24! = 147926426347074375 << 22 +]; + +impl BigUint { + pub fn factorial(&self) -> BigUint { + // usize::MAX! > 2**usize::MAX, so just error w/o crashing the system + match self.to_usize() { + Option::Some(x) => factorial(x), + Option::None => panic!("{} is too large to calculate the factorial", self), + } + } +} + +fn factorial(x: usize) -> BigUint { + if let Option::Some(&x) = SMALL_FACTORIALS.get(x) { + x.into() + } else { + // This uses a neat trick: The number of 2s that are a factor of x! is + // equal to x - x.count_ones() + let result = odd_factorial(x); + let power_two_count = x - x.count_ones() as usize; + result << power_two_count + } +} + +#[inline] +fn odd_factorial(mut x: usize) -> BigUint { + if let Option::Some(&x) = ODD_FACTORIALS.get(x) { + return x.into(); + } + // Algorithm: `odd_factorial(n) = odds(1..=n) * odd_factorial(n / 2) + // Prime sieving could be used to improve this in the future, see + // https://gmplib.org/manual/Factorial-Algorithm + // This uses a loop instead of recursion, to improve speed & reduce stack + // space + let mut result = BigUint::one(); + while x as usize > ODD_FACTORIALS.len() { + let mut i = 3; + while i <= x { + result *= i; + i += 2; + } + x /= 2; + } + result * ODD_FACTORIALS[x as usize] +} From 65528a365141cc949b9612c97514a1636b8cfb63 Mon Sep 17 00:00:00 2001 From: Patrick Norton Date: Sun, 7 Nov 2021 11:22:02 -0800 Subject: [PATCH 2/3] Added factorial test --- src/biguint.rs | 9 +++++++++ src/biguint/factorial.rs | 16 +++------------- tests/biguint.rs | 18 ++++++++++++++++++ 3 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/biguint.rs b/src/biguint.rs index 9d898ca4..3b2bdcb2 100644 --- a/src/biguint.rs +++ b/src/biguint.rs @@ -955,6 +955,15 @@ impl BigUint { self.normalize(); } } + + /// Computes the factorial of the `BigUInt`. + pub fn factorial(&self) -> BigUint { + // usize::MAX! > 2**usize::MAX, so just error w/o running the system out of memory + match self.to_usize() { + Option::Some(x) => factorial::factorial(x), + Option::None => panic!("{} is too large to calculate the factorial", self), + } + } } pub(crate) trait IntDigits { diff --git a/src/biguint/factorial.rs b/src/biguint/factorial.rs index e813114c..55ed3ed3 100644 --- a/src/biguint/factorial.rs +++ b/src/biguint/factorial.rs @@ -1,5 +1,5 @@ use crate::BigUint; -use num_traits::{One, ToPrimitive}; +use num_traits::One; /// All factorials that fit into a u64 const SMALL_FACTORIALS: [u64; 21] = [ @@ -56,17 +56,7 @@ const ODD_FACTORIALS: [u64; 25] = [ 147926426347074375, // 24! = 147926426347074375 << 22 ]; -impl BigUint { - pub fn factorial(&self) -> BigUint { - // usize::MAX! > 2**usize::MAX, so just error w/o crashing the system - match self.to_usize() { - Option::Some(x) => factorial(x), - Option::None => panic!("{} is too large to calculate the factorial", self), - } - } -} - -fn factorial(x: usize) -> BigUint { +pub fn factorial(x: usize) -> BigUint { if let Option::Some(&x) = SMALL_FACTORIALS.get(x) { x.into() } else { @@ -89,7 +79,7 @@ fn odd_factorial(mut x: usize) -> BigUint { // This uses a loop instead of recursion, to improve speed & reduce stack // space let mut result = BigUint::one(); - while x as usize > ODD_FACTORIALS.len() { + while x as usize >= ODD_FACTORIALS.len() { let mut i = 3; while i <= x { result *= i; diff --git a/tests/biguint.rs b/tests/biguint.rs index 821b754a..73dc9a3f 100644 --- a/tests/biguint.rs +++ b/tests/biguint.rs @@ -1862,3 +1862,21 @@ fn test_set_bit() { x.set_bit(1, false); assert_eq!(x, BigUint::zero()); } + +const FACTORIAL_100_BYTES: &[u32] = &[ + 0x00000000, 0x00000000, 0x00000000, 0x2735c61a, 0xee8b02ea, 0xb3b72ed2, 0x9420c6ec, 0x45570cca, + 0xdf103917, 0x943a321c, 0xeb21b5b2, 0x66ef9a70, 0xa40d16e9, 0x28d54bbd, 0xdc240695, 0x964ec395, + 0x00001b30, +]; + +#[test] +fn test_factorial() { + let mut x = BigUint::zero(); + assert_eq!(x.factorial(), BigUint::one()); + x = BigUint::one(); + assert_eq!(x.factorial(), BigUint::one()); + x = BigUint::from(10u8); + assert_eq!(x.factorial(), BigUint::from(3628800u32)); + x = BigUint::from(100u8); + assert_eq!(x.factorial(), BigUint::from_slice(FACTORIAL_100_BYTES)); +} From 7d7db9cef25485e00b65449b44f14bfd73931a36 Mon Sep 17 00:00:00 2001 From: PatrickNorton <36009535+PatrickNorton@users.noreply.github.com> Date: Fri, 15 Apr 2022 21:48:25 -0700 Subject: [PATCH 3/3] Remove unnecessary `Option::` Co-authored-by: pierwill <19642016+pierwill@users.noreply.github.com> --- src/biguint/factorial.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biguint/factorial.rs b/src/biguint/factorial.rs index 55ed3ed3..516b6620 100644 --- a/src/biguint/factorial.rs +++ b/src/biguint/factorial.rs @@ -57,7 +57,7 @@ const ODD_FACTORIALS: [u64; 25] = [ ]; pub fn factorial(x: usize) -> BigUint { - if let Option::Some(&x) = SMALL_FACTORIALS.get(x) { + if let Some(&x) = SMALL_FACTORIALS.get(x) { x.into() } else { // This uses a neat trick: The number of 2s that are a factor of x! is