Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added BigUInt::factorial #229

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions src/biguint.rs
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ mod subtraction;

mod bits;
mod convert;
mod factorial;
mod iter;
mod monty;
mod power;
Expand Down Expand Up @@ -954,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 {
Expand Down
91 changes: 91 additions & 0 deletions src/biguint/factorial.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
use crate::BigUint;
use num_traits::One;

/// 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
];

pub fn factorial(x: usize) -> BigUint {
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
// 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]
}
18 changes: 18 additions & 0 deletions tests/biguint.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}