|
| 1 | +use once_cell::sync::Lazy; |
| 2 | +use std::collections::HashMap; |
| 3 | +use std::convert::TryFrom; |
| 4 | + |
| 5 | +const MARKER_TABLE_SIZE: usize = 200_000; |
| 6 | + |
| 7 | +// TODO: replace with const fn when it is possible |
| 8 | +// (for and if are not allowed in const fn on current stable) |
| 9 | +// https://github.com/rust-lang/rust/issues/87575 |
| 10 | +static MARKER_TABLES: Lazy<HashMap<u8, Vec<u128>>> = Lazy::new(|| { |
| 11 | + let mut m = HashMap::new(); |
| 12 | + for k in 1..10u8 { |
| 13 | + let mut table = vec![0u128; MARKER_TABLE_SIZE]; |
| 14 | + let table_size = if k == 1 { |
| 15 | + 128 |
| 16 | + } else if k == 2 { |
| 17 | + 8128 |
| 18 | + } else { |
| 19 | + table.len() |
| 20 | + }; |
| 21 | + |
| 22 | + table[0] = ((1 << k) - 1) as u128; |
| 23 | + for i in 1..table_size { |
| 24 | + table[i] = next_rank(table[i - 1]); |
| 25 | + } |
| 26 | + m.insert(k, table); |
| 27 | + } |
| 28 | + m |
| 29 | +}); |
| 30 | + |
| 31 | +/// https://en.wikipedia.org/wiki/Combinatorial_number_system |
| 32 | +pub fn rank(value: usize, k: u8) -> u128 { |
| 33 | + assert!(k > 0 && k < 10, "kappa needs to be less than 10"); |
| 34 | + // it's possible this may overflow if value > (128 choose k) or return |
| 35 | + // a bad value (0) if value > (128 choose k) and k == 1 or 2 |
| 36 | + if value >= MARKER_TABLE_SIZE { |
| 37 | + let mut marker = MARKER_TABLES[&k][MARKER_TABLE_SIZE - 1]; |
| 38 | + for _ in 0..(value - MARKER_TABLE_SIZE) { |
| 39 | + // next_rank would overflow if we pass 0, we return it instead |
| 40 | + if marker == 0 { |
| 41 | + return marker; |
| 42 | + } |
| 43 | + marker = next_rank(marker); |
| 44 | + } |
| 45 | + marker |
| 46 | + } else { |
| 47 | + MARKER_TABLES[&k][value] |
| 48 | + } |
| 49 | +} |
| 50 | + |
| 51 | +/// https://en.wikipedia.org/wiki/Combinatorial_number_system |
| 52 | +pub fn unrank(marker: u128) -> usize { |
| 53 | + // val = choose(rank(0), 1) + choose(rank(1), 2) + choose(rank(2), 3) + ... |
| 54 | + let mut working_marker = marker; |
| 55 | + let mut value = 0u64; |
| 56 | + let mut idx = 0; |
| 57 | + while working_marker != 0 { |
| 58 | + let rank = u64::from(working_marker.trailing_zeros()); |
| 59 | + working_marker -= 1 << rank; |
| 60 | + idx += 1; |
| 61 | + value += choose(rank, idx); |
| 62 | + } |
| 63 | + value as usize |
| 64 | +} |
| 65 | + |
| 66 | +/// (Hopefully) fast implementation of a binomial |
| 67 | +/// |
| 68 | +/// This uses a preset group of equations for k < 8 and then falls back to a |
| 69 | +/// multiplicative implementation that tries to prevent overflows while |
| 70 | +/// maintaining all results as exact integers. |
| 71 | +#[inline] |
| 72 | +pub fn choose(n: u64, k: u8) -> u64 { |
| 73 | + // (extra border condition for speed-up?) |
| 74 | + // if n == u64::from(k) { |
| 75 | + // return 1; |
| 76 | + // } |
| 77 | + match k { |
| 78 | + 0 => 1, |
| 79 | + 1 => n, |
| 80 | + 2 => n * (n - 1) / 2, |
| 81 | + 3 => n * (n - 1) * (n - 2) / 6, |
| 82 | + 4 => n * (n - 1) * (n - 2) * (n - 3) / 24, |
| 83 | + 5 => n * (n - 1) * (n - 2) * (n - 3) * (n - 4) / 120, |
| 84 | + 6 => n * (n - 1) * (n - 2) * (n - 3) * (n - 4) * (n - 5) / 720, |
| 85 | + 7 => n * (n - 1) * (n - 2) * (n - 3) * (n - 4) * (n - 5) * (n - 6) / 5040, |
| 86 | + _ => { |
| 87 | + let mut num: u128 = 1; |
| 88 | + let mut denom: u128 = 1; |
| 89 | + for i in 1..=u128::from(k) { |
| 90 | + num *= u128::from(n) + 1 - i; |
| 91 | + if num % i == 0 { |
| 92 | + num /= i; |
| 93 | + continue; |
| 94 | + } |
| 95 | + denom *= i; |
| 96 | + if num % denom == 0 { |
| 97 | + num /= denom; |
| 98 | + denom = 1; |
| 99 | + } |
| 100 | + } |
| 101 | + TryFrom::try_from(num / denom) |
| 102 | + .unwrap_or_else(|_| panic!("{} choose {} is greater than 2**64", n, k)) |
| 103 | + // (or recursively) choose(n - 1, k - 1) + choose(n-1, k) |
| 104 | + // for floats, this should work since they handle fractions: |
| 105 | + // (1..u64::from(k)).map(|i| (n + 1 - i) / i).product(), |
| 106 | + } |
| 107 | + } |
| 108 | +} |
| 109 | + |
| 110 | +#[inline] |
| 111 | +fn next_rank(marker: u128) -> u128 { |
| 112 | + if marker == 0 { |
| 113 | + unreachable!("Got next_rank called with marker == 0"); |
| 114 | + } |
| 115 | + let t = marker | (marker - 1); |
| 116 | + (t + 1) | (((!t & (t + 1)) - 1) >> (marker.trailing_zeros() + 1)) |
| 117 | +} |
| 118 | + |
| 119 | +#[cfg(test)] |
| 120 | +mod tests { |
| 121 | + use super::*; |
| 122 | + |
| 123 | + #[test] |
| 124 | + fn test_rank() { |
| 125 | + assert_eq!(rank(0, 3), 7); |
| 126 | + assert_eq!(rank(2, 3), 13); |
| 127 | + assert_eq!(rank(0, 3).count_ones(), 3); |
| 128 | + assert_eq!(rank(2, 3).count_ones(), 3); |
| 129 | + assert_eq!(rank(35001, 4).count_ones(), 4); |
| 130 | + |
| 131 | + // Maximum value of 64 choose 3 |
| 132 | + assert_eq!(rank(41663, 3).count_ones(), 3); |
| 133 | + } |
| 134 | + |
| 135 | + #[test] |
| 136 | + fn test_unrank() { |
| 137 | + // 3 bit markers |
| 138 | + assert_eq!(unrank(7), 0); |
| 139 | + assert_eq!(unrank(13), 2); |
| 140 | + } |
| 141 | + |
| 142 | + #[test] |
| 143 | + fn test_rank_and_unrank() { |
| 144 | + for k in 1..4u8 { |
| 145 | + for value in [1 as usize, 23, 45].iter() { |
| 146 | + assert_eq!(unrank(rank(*value, k)), *value); |
| 147 | + } |
| 148 | + } |
| 149 | + } |
| 150 | + |
| 151 | + #[test] |
| 152 | + fn test_choose() { |
| 153 | + assert_eq!(choose(1, 1), 1); |
| 154 | + assert_eq!(choose(10, 1), 10); |
| 155 | + |
| 156 | + assert_eq!(choose(5, 2), 10); |
| 157 | + |
| 158 | + assert_eq!(choose(5, 3), 10); |
| 159 | + |
| 160 | + assert_eq!(choose(5, 4), 5); |
| 161 | + |
| 162 | + assert_eq!(choose(5, 5), 1); |
| 163 | + assert_eq!(choose(20, 5), 15504); |
| 164 | + |
| 165 | + assert_eq!(choose(20, 6), 38760); |
| 166 | + |
| 167 | + assert_eq!(choose(20, 7), 77520); |
| 168 | + assert_eq!(choose(23, 7), 245157); |
| 169 | + |
| 170 | + // test the last branch |
| 171 | + assert_eq!(choose(8, 8), 1); |
| 172 | + assert_eq!(choose(9, 8), 9); |
| 173 | + |
| 174 | + // every value of 64 choose n should work |
| 175 | + assert_eq!(choose(64, 0), 1); |
| 176 | + assert_eq!(choose(64, 1), 64); |
| 177 | + assert_eq!(choose(64, 16), 488526937079580); |
| 178 | + assert_eq!(choose(64, 32), 1832624140942590534); |
| 179 | + assert_eq!(choose(64, 48), 488526937079580); |
| 180 | + assert_eq!(choose(64, 63), 64); |
| 181 | + assert_eq!(choose(64, 64), 1); |
| 182 | + |
| 183 | + // super high values can overflow; these are approaching the limit |
| 184 | + assert_eq!(choose(128, 11), 2433440563030400); |
| 185 | + assert_eq!(choose(128, 13), 211709328983644800); |
| 186 | + assert_eq!(choose(256, 9), 11288510714272000); |
| 187 | + } |
| 188 | + |
| 189 | + #[test] |
| 190 | + #[should_panic(expected = "256 choose 20 is greater than 2**64")] |
| 191 | + fn test_choose_overflow() { |
| 192 | + assert_eq!(choose(256, 20), 11288510714272000); |
| 193 | + } |
| 194 | + |
| 195 | + #[test] |
| 196 | + fn test_next_rank() { |
| 197 | + assert_eq!(next_rank(0b1), 0b10); |
| 198 | + assert_eq!(next_rank(0b100), 0b1000); |
| 199 | + |
| 200 | + assert_eq!(next_rank(0b111), 0b1011); |
| 201 | + assert_eq!(next_rank(0b1000101), 0b1000110); |
| 202 | + } |
| 203 | +} |
0 commit comments