randperm_crt/
lib.rs

1#![warn(clippy::must_use_candidate)]
2#![deny(clippy::use_self)]
3#![deny(clippy::double_must_use)]
4#![deny(clippy::if_not_else)]
5#![deny(clippy::inconsistent_struct_constructor)]
6#![deny(clippy::iter_not_returning_iterator)]
7#![deny(clippy::map_unwrap_or)]
8#![deny(clippy::mod_module_files)]
9#![deny(clippy::semicolon_if_nothing_returned)]
10
11mod crt;
12
13use rand::Rng;
14
15#[derive(Debug, Clone, PartialEq, Eq)]
16struct FactoredInteger {
17    factors: Vec<(u8, u8)>,
18}
19
20impl FactoredInteger {
21    fn new(mut n: u64) -> Option<Self> {
22        let mut factors = Vec::new();
23
24        let pow2 = n.trailing_zeros() as u8;
25        if pow2 != 0 {
26            n >>= pow2;
27            factors.push((2, pow2));
28        }
29
30        for p in (3..u8::MAX).step_by(2) {
31            let q = p as u64;
32
33            let mut counter = 0;
34            while n % q == 0 {
35                counter += 1;
36                n /= q;
37            }
38
39            if counter > 0 {
40                factors.push((p, counter));
41            }
42
43            if n == 1 {
44                break;
45            }
46        }
47
48        if n == 1 {
49            Some(Self { factors })
50        } else {
51            None
52        }
53    }
54}
55
56pub trait Permutation: Sized {
57    fn num_points(&self) -> u64;
58    fn nth(&self, n: u64) -> Option<u64>;
59
60    fn iter(&self) -> PermutationIter<'_, Self> {
61        PermutationIter { perm: self, idx: 0 }
62    }
63}
64
65#[derive(Debug, Clone, PartialEq, Eq)]
66pub struct RandomPermutation {
67    num_points: u64,
68    sub_perms: Vec<Vec<u64>>,
69}
70
71impl RandomPermutation {
72    #[cfg(feature = "thread_rng")]
73    #[must_use]
74    pub fn new(n: u64) -> Option<Self> {
75        Self::with_rng(n, &mut rand::rng())
76    }
77
78    pub fn with_rng<R: Rng>(n: u64, rng: &mut R) -> Option<Self> {
79        let factored_n = FactoredInteger::new(n)?;
80        let num_prime_powers = factored_n.factors.len();
81
82        let mut order = (0..num_prime_powers).collect::<Vec<_>>();
83        for a in 0..num_prime_powers {
84            let b = rng.random_range(a..num_prime_powers);
85            order.swap(a, b);
86        }
87
88        let sub_perms = (0..num_prime_powers)
89            .map(|i| {
90                let (p, k) = factored_n.factors[order[i]];
91                let pk = (p as u64).pow(k as u32);
92                let mut vec = (0..pk).collect::<Vec<_>>();
93
94                let pk = pk as usize;
95                for a in 0..pk {
96                    let b = rng.random_range(a..pk);
97                    vec.swap(a, b);
98                }
99
100                vec
101            })
102            .collect();
103
104        Some(Self {
105            num_points: n,
106            sub_perms,
107        })
108    }
109
110    #[must_use]
111    pub fn inverse(&self) -> Inverse<'_> {
112        Inverse { perm: self }
113    }
114}
115
116impl Permutation for RandomPermutation {
117    fn num_points(&self) -> u64 {
118        self.num_points
119    }
120
121    fn nth(&self, mut n: u64) -> Option<u64> {
122        if n >= self.num_points {
123            return None;
124        }
125
126        let remainders = self.sub_perms.iter().fold(Vec::new(), |mut rem, perm| {
127            let pk = perm.len() as u64;
128            rem.push(perm[(n % pk) as usize]);
129            n /= pk;
130            rem
131        });
132
133        let moduli = self
134            .sub_perms
135            .iter()
136            .map(|perm| perm.len() as u64)
137            .collect::<Vec<_>>();
138
139        Some(crt::chinese_remainder(&remainders, &moduli).unwrap())
140    }
141}
142
143pub struct Inverse<'a> {
144    perm: &'a RandomPermutation,
145}
146
147impl Permutation for Inverse<'_> {
148    fn num_points(&self) -> u64 {
149        self.perm.num_points()
150    }
151
152    fn nth(&self, n: u64) -> Option<u64> {
153        if n >= self.num_points() {
154            None
155        } else {
156            Some(self.perm.sub_perms.iter().rev().fold(0, |idx, perm| {
157                let pk = perm.len() as u64;
158                let pos = perm.iter().position(|&a| a == n % pk).unwrap() as u64;
159                idx * pk + pos
160            }))
161        }
162    }
163}
164
165pub struct PermutationIter<'a, P: Permutation> {
166    perm: &'a P,
167    idx: u64,
168}
169
170impl<P: Permutation> Iterator for PermutationIter<'_, P> {
171    type Item = u64;
172
173    fn next(&mut self) -> Option<Self::Item> {
174        let a = self.perm.nth(self.idx);
175        self.idx += 1;
176        a
177    }
178
179    fn nth(&mut self, n: usize) -> Option<Self::Item> {
180        self.idx += n as u64;
181        self.perm.nth(self.idx)
182    }
183}
184
185pub struct Composition<'a> {
186    perms: &'a [RandomPermutation],
187}
188
189impl<'a> Composition<'a> {
190    #[must_use]
191    pub fn new(perms: &'a [RandomPermutation]) -> Option<Self> {
192        if perms.is_empty() {
193            return None;
194        }
195
196        if !perms
197            .iter()
198            .all(|p| p.num_points() == perms[0].num_points())
199        {
200            return None;
201        }
202
203        Some(Self { perms })
204    }
205}
206
207impl Permutation for Composition<'_> {
208    fn num_points(&self) -> u64 {
209        self.perms[0].num_points()
210    }
211
212    fn nth(&self, n: u64) -> Option<u64> {
213        self.perms.iter().try_fold(n, |n, perm| perm.nth(n))
214    }
215}
216
217#[cfg(test)]
218mod tests {
219    use rand::SeedableRng;
220    use rand_xoshiro::Xoshiro256StarStar;
221
222    use crate::*;
223
224    const SEED: [u8; 32] = [
225        144, 115, 104, 224, 226, 59, 231, 208, 100, 18, 137, 138, 234, 236, 129, 82, 184, 196, 19,
226        43, 145, 94, 60, 77, 184, 198, 244, 164, 174, 224, 59, 152,
227    ];
228
229    mod factored_integer {
230        use super::*;
231
232        #[test]
233        fn test_new_1() {
234            let n = FactoredInteger::new(14237396402848819200);
235            assert_eq!(
236                n,
237                Some(FactoredInteger {
238                    factors: vec![
239                        (2, 11),
240                        (3, 3),
241                        (5, 2),
242                        (7, 3),
243                        (11, 4),
244                        (13, 1),
245                        (19, 3),
246                        (23, 1)
247                    ]
248                })
249            );
250        }
251
252        #[test]
253        fn test_new_2() {
254            let n = FactoredInteger::new(8929777156897433877);
255            assert_eq!(
256                n,
257                Some(FactoredInteger {
258                    factors: vec![(3, 25), (199, 1), (211, 1), (251, 1)]
259                })
260            );
261        }
262
263        #[test]
264        fn test_new_3() {
265            let n = FactoredInteger::new(2u64.pow(63));
266            assert_eq!(
267                n,
268                Some(FactoredInteger {
269                    factors: vec![(2, 63)]
270                })
271            );
272        }
273
274        #[test]
275        fn test_new_4() {
276            let n = FactoredInteger::new(257);
277            assert_eq!(n, None);
278        }
279
280        #[test]
281        fn test_new_5() {
282            let n = FactoredInteger::new(1297068779 * 3196491187);
283            assert_eq!(n, None);
284        }
285    }
286
287    mod random_permutation {
288        use super::*;
289
290        #[test]
291        fn test_random_permutation() {
292            let mut rng = Xoshiro256StarStar::from_seed(SEED);
293
294            let p = RandomPermutation::with_rng(362880, &mut rng).unwrap();
295            let mut vec = p.iter().collect::<Vec<_>>();
296            vec.sort();
297
298            assert!(vec.iter().copied().eq(0..362880));
299        }
300
301        #[test]
302        fn test_nth() {
303            let mut rng = Xoshiro256StarStar::from_seed(SEED);
304            let p = RandomPermutation::with_rng((1..=20).product(), &mut rng).unwrap();
305
306            assert_eq!(p.nth(0), Some(1319942158626894223));
307            assert_eq!(p.nth(1), Some(2412265509236814223));
308            assert_eq!(p.nth(2), Some(2263312325062734223));
309            assert_eq!(p.nth(1000000000000000000), Some(320224176323398978));
310            assert_eq!(p.nth(1146506666437230775), Some(0));
311            assert_eq!(p.nth(2432902008176639999), Some(725036512928904481));
312            assert_eq!(p.nth(2432902008176640000), None);
313            assert_eq!(p.nth(u64::MAX), None);
314        }
315
316        #[test]
317        fn test_nth_2() {
318            let mut rng = Xoshiro256StarStar::from_seed(SEED);
319            let p = RandomPermutation::with_rng(300, &mut rng).unwrap();
320            let v = p.iter().collect::<Vec<_>>();
321
322            assert_eq!(
323                v,
324                &[
325                    176, 276, 76, 101, 201, 1, 26, 126, 226, 251, 51, 151, 200, 0, 100, 125, 225,
326                    25, 50, 150, 250, 275, 75, 175, 152, 252, 52, 77, 177, 277, 2, 102, 202, 227,
327                    27, 127, 224, 24, 124, 149, 249, 49, 74, 174, 274, 299, 99, 199, 8, 108, 208,
328                    233, 33, 133, 158, 258, 58, 83, 183, 283, 284, 84, 184, 209, 9, 109, 134, 234,
329                    34, 59, 159, 259, 236, 36, 136, 161, 261, 61, 86, 186, 286, 11, 111, 211, 104,
330                    204, 4, 29, 129, 229, 254, 54, 154, 179, 279, 79, 212, 12, 112, 137, 237, 37,
331                    62, 162, 262, 287, 87, 187, 272, 72, 172, 197, 297, 97, 122, 222, 22, 47, 147,
332                    247, 68, 168, 268, 293, 93, 193, 218, 18, 118, 143, 243, 43, 56, 156, 256, 281,
333                    81, 181, 206, 6, 106, 131, 231, 31, 80, 180, 280, 5, 105, 205, 230, 30, 130,
334                    155, 255, 55, 44, 144, 244, 269, 69, 169, 194, 294, 94, 119, 219, 19, 128, 228,
335                    28, 53, 153, 253, 278, 78, 178, 203, 3, 103, 20, 120, 220, 245, 45, 145, 170,
336                    270, 70, 95, 195, 295, 164, 264, 64, 89, 189, 289, 14, 114, 214, 239, 39, 139,
337                    296, 96, 196, 221, 21, 121, 146, 246, 46, 71, 171, 271, 260, 60, 160, 185, 285,
338                    85, 110, 210, 10, 35, 135, 235, 116, 216, 16, 41, 141, 241, 266, 66, 166, 191,
339                    291, 91, 92, 192, 292, 17, 117, 217, 242, 42, 142, 167, 267, 67, 248, 48, 148,
340                    173, 273, 73, 98, 198, 298, 23, 123, 223, 32, 132, 232, 257, 57, 157, 182, 282,
341                    82, 107, 207, 7, 140, 240, 40, 65, 165, 265, 290, 90, 190, 215, 15, 115, 188,
342                    288, 88, 113, 213, 13, 38, 138, 238, 263, 63, 163
343                ]
344            );
345        }
346    }
347
348    mod inverse {
349        use super::*;
350
351        #[test]
352        fn test_nth_1() {
353            let mut rng = Xoshiro256StarStar::from_seed(SEED);
354            let p = RandomPermutation::with_rng((1..=20).product(), &mut rng).unwrap();
355            let inv = p.inverse();
356
357            assert_eq!(inv.nth(1319942158626894223), Some(0));
358            assert_eq!(inv.nth(2412265509236814223), Some(1));
359            assert_eq!(inv.nth(2263312325062734223), Some(2));
360            assert_eq!(inv.nth(320224176323398978), Some(1000000000000000000));
361            assert_eq!(inv.nth(0), Some(1146506666437230775));
362            assert_eq!(inv.nth(725036512928904481), Some(2432902008176639999));
363            assert_eq!(inv.nth(2432902008176640000), None);
364            assert_eq!(inv.nth(u64::MAX), None);
365        }
366
367        #[test]
368        fn test_nth_2() {
369            let mut rng = Xoshiro256StarStar::from_seed(SEED);
370            let p = RandomPermutation::with_rng(300, &mut rng).unwrap();
371            let inv = p.inverse();
372            let v = inv.iter().collect::<Vec<_>>();
373
374            assert_eq!(
375                v,
376                &[
377                    13, 5, 30, 178, 86, 147, 139, 275, 48, 64, 224, 81, 97, 293, 198, 286, 230,
378                    243, 127, 167, 180, 208, 116, 261, 37, 17, 6, 34, 170, 87, 151, 143, 264, 52,
379                    68, 225, 73, 101, 294, 202, 278, 231, 247, 131, 156, 184, 212, 117, 253, 41,
380                    18, 10, 26, 171, 91, 155, 132, 268, 56, 69, 217, 77, 102, 298, 194, 279, 235,
381                    251, 120, 160, 188, 213, 109, 257, 42, 22, 2, 27, 175, 95, 144, 136, 272, 57,
382                    61, 221, 78, 106, 290, 195, 283, 239, 240, 124, 164, 189, 205, 113, 258, 46,
383                    14, 3, 31, 179, 84, 148, 140, 273, 49, 65, 222, 82, 98, 291, 199, 287, 228,
384                    244, 128, 165, 181, 209, 114, 262, 38, 15, 7, 35, 168, 88, 152, 141, 265, 53,
385                    66, 226, 74, 99, 295, 203, 276, 232, 248, 129, 157, 185, 210, 118, 254, 39, 19,
386                    11, 24, 172, 92, 153, 133, 269, 54, 70, 218, 75, 103, 299, 192, 280, 236, 249,
387                    121, 161, 186, 214, 110, 255, 43, 23, 0, 28, 176, 93, 145, 137, 270, 58, 62,
388                    219, 79, 107, 288, 196, 284, 237, 241, 125, 162, 190, 206, 111, 259, 47, 12, 4,
389                    32, 177, 85, 149, 138, 274, 50, 63, 223, 83, 96, 292, 200, 285, 229, 245, 126,
390                    166, 182, 207, 115, 263, 36, 16, 8, 33, 169, 89, 150, 142, 266, 51, 67, 227,
391                    72, 100, 296, 201, 277, 233, 246, 130, 158, 183, 211, 119, 252, 40, 20, 9, 25,
392                    173, 90, 154, 134, 267, 55, 71, 216, 76, 104, 297, 193, 281, 234, 250, 122,
393                    159, 187, 215, 108, 256, 44, 21, 1, 29, 174, 94, 146, 135, 271, 59, 60, 220,
394                    80, 105, 289, 197, 282, 238, 242, 123, 163, 191, 204, 112, 260, 45
395                ]
396            );
397        }
398
399        #[test]
400        fn test_nth_3() {
401            let mut rng = Xoshiro256StarStar::from_seed(SEED);
402            let p = RandomPermutation::with_rng((1..=20).product(), &mut rng).unwrap();
403            let inv = p.inverse();
404
405            for i in 0..1000 {
406                assert_eq!(p.nth(inv.nth(i).unwrap()), Some(i));
407            }
408        }
409    }
410
411    mod iterator {
412        use super::*;
413
414        #[test]
415        fn test_next() {
416            let mut rng = Xoshiro256StarStar::from_seed(SEED);
417            let p = RandomPermutation::with_rng(3113510400, &mut rng).unwrap();
418            let mut iter = p.iter();
419
420            for i in 0..1000 {
421                assert_eq!(iter.next(), p.nth(i));
422            }
423        }
424
425        #[test]
426        fn test_nth() {
427            let mut rng = Xoshiro256StarStar::from_seed(SEED);
428            let p = RandomPermutation::with_rng(3113510400, &mut rng).unwrap();
429            let mut iter = p.iter();
430
431            for i in 0..1000 {
432                assert_eq!(iter.nth(1000000), p.nth((i + 1) * 1000000));
433            }
434        }
435    }
436
437    mod composition {
438        use super::*;
439
440        #[test]
441        fn test_new_1() {
442            let comp = Composition::new(&[]);
443            assert!(comp.is_none());
444        }
445
446        #[test]
447        fn test_new_2() {
448            let mut rng = Xoshiro256StarStar::from_seed(SEED);
449            let p1 = RandomPermutation::with_rng(300, &mut rng).unwrap();
450            let p2 = RandomPermutation::with_rng(400, &mut rng).unwrap();
451
452            let v = vec![p1, p2];
453            let comp = Composition::new(&v);
454
455            assert!(comp.is_none());
456        }
457
458        #[test]
459        fn test_nth() {
460            let mut rng = Xoshiro256StarStar::from_seed(SEED);
461            let p1 = RandomPermutation::with_rng(300, &mut rng).unwrap();
462            let p2 = RandomPermutation::with_rng(300, &mut rng).unwrap();
463
464            let v = vec![p1.clone(), p2.clone()];
465            let comp = Composition::new(&v).unwrap();
466
467            for i in 0..300 {
468                assert_eq!(comp.nth(i), p2.nth(p1.nth(i).unwrap()));
469            }
470        }
471    }
472}