smart_big_rational/
denom_sparse.rs

1// Copyright 2026 The SmartBigRational Authors
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7//     https://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15use crate::primes::{
16    ODD_PRIME_DIVIDERS_U16, ODD_PRIME_DIVIDERS_U32, ODD_PRIME_DIVIDERS_U64, ODD_PRIMES,
17    known_odd_prime_factor_indices,
18};
19use crate::util::OddDivider;
20use crate::{Denom, DenomRef};
21use num_bigint::{BigInt, BigUint, Sign};
22use num_integer::Integer;
23use num_traits::{One, Zero};
24use smallvec::SmallVec;
25use std::cmp::Ordering;
26use std::iter::Peekable;
27use std::ops::{Div, DivAssign, Mul, MulAssign};
28
29/// Denominator representation that decomposes an integer as a product of the
30/// first `NUM_PRIMES` primes (up to 2^16), multiplied by a regular big integer
31/// when that's not sufficient.
32#[derive(Clone, Debug, PartialEq, Eq)]
33pub struct DenomSparseU16<const NUM_PRIMES: usize, const NUM_INLINE: usize> {
34    // Invariant: the representation is in canonical form, i.e. powers of small primes must
35    // saturate the primes array before overflowing into the remainder.
36    primes: SmallVec<[(u16, u16); NUM_INLINE]>,
37    remainder: Option<BigUint>,
38}
39
40/// Denominator representation that decomposes an integer as a product of the
41/// first 6542 primes (up to 0xfff1), multiplied by a regular big integer when
42/// that's not sufficient.
43pub type DenomSparse6542 = DenomSparseU16<6542, 8>;
44
45impl<const NUM_PRIMES: usize, const NUM_INLINE: usize>
46    DenomRef<DenomSparseU16<NUM_PRIMES, NUM_INLINE>> for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
47{
48}
49
50impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Denom
51    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
52{
53    const ONE: Self = Self {
54        primes: SmallVec::new_const(),
55        remainder: None,
56    };
57
58    fn into_biguint(self) -> BigUint {
59        self.into()
60    }
61
62    fn to_biguint(&self) -> BigUint {
63        self.into()
64    }
65
66    fn normalize(lnum: &mut BigInt, rnum: &mut BigInt, ldenom: &Self, rdenom: &Self) -> Self {
67        let mut primes = SmallVec::new();
68        let mut ltmp = 1_usize;
69        let mut rtmp = 1_usize;
70
71        for (p, (lcount, rcount)) in Zip(
72            ldenom.primes.iter().copied().peekable(),
73            rdenom.primes.iter().copied().peekable(),
74        ) {
75            match lcount.cmp(&rcount) {
76                Ordering::Equal => {
77                    primes.push((p, lcount));
78                }
79                Ordering::Less => {
80                    Self::accum_pow(lnum, &mut ltmp, p, rcount - lcount);
81                    primes.push((p, rcount));
82                }
83                Ordering::Greater => {
84                    Self::accum_pow(rnum, &mut rtmp, p, lcount - rcount);
85                    primes.push((p, lcount));
86                }
87            }
88        }
89
90        *lnum *= ltmp;
91        *rnum *= rtmp;
92        let remainder = match (&ldenom.remainder, &rdenom.remainder) {
93            (None, None) => None,
94            (None, Some(r)) => {
95                *lnum *= BigInt::from_biguint(Sign::Plus, r.clone());
96                Some(r.clone())
97            }
98            (Some(l), None) => {
99                *rnum *= BigInt::from_biguint(Sign::Plus, l.clone());
100                Some(l.clone())
101            }
102            (Some(l), Some(r)) => {
103                if l == r {
104                    Some(l.clone())
105                } else {
106                    *lnum *= BigInt::from_biguint(Sign::Plus, r.clone());
107                    *rnum *= BigInt::from_biguint(Sign::Plus, l.clone());
108                    Some(l * r)
109                }
110            }
111        };
112        Self { primes, remainder }
113    }
114
115    fn gcd_reduce(&mut self, num: &mut BigInt) {
116        if let Some(remainder) = &mut self.remainder {
117            let gcd = remainder.gcd(num.magnitude());
118            if !gcd.is_one() {
119                *remainder /= &gcd;
120                if remainder.is_one() {
121                    self.remainder = None;
122                }
123                let gcd: BigInt = gcd.into();
124                *num /= gcd;
125            }
126        }
127
128        'outer: for (p, i) in self.primes.iter_mut() {
129            let p = BigInt::from(*p);
130            while *i != 0 {
131                let (quo, rem) = num.div_rem(&p);
132                if !rem.is_zero() {
133                    break;
134                }
135                *num = quo;
136                *i -= 1;
137                if num.magnitude().is_one() {
138                    break 'outer;
139                }
140            }
141        }
142        self.primes.retain(|(_, i)| *i != 0);
143    }
144}
145
146impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
147    fn decompose(mut x: BigUint) -> Self {
148        const {
149            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
150        }
151
152        let bits = x.bits();
153        if bits <= 8 {
154            return Self::decompose_u8(x.try_into().unwrap());
155        } else if bits <= 16 {
156            return Self::decompose_u16(x.try_into().unwrap());
157        } else if bits <= 32 {
158            return Self::decompose_u32(x.try_into().unwrap());
159        } else if bits <= 64 {
160            return Self::decompose_u64(x.try_into().unwrap());
161        } else if bits <= 128 {
162            return Self::decompose_u128(x.try_into().unwrap());
163        }
164
165        let mut primes = SmallVec::new();
166
167        let mut count2 = x.trailing_zeros().unwrap();
168        if count2 != 0 {
169            x >>= count2;
170            if count2 <= u16::MAX as u64 {
171                primes.push((2, count2 as u16));
172                count2 = 0;
173            } else {
174                primes.push((2, u16::MAX));
175                count2 -= u16::MAX as u64;
176            }
177        }
178
179        'outer: for &p in ODD_PRIMES.iter().take(NUM_PRIMES - 1) {
180            let bigp = BigUint::from(p);
181            let mut count = 0;
182            while count != u16::MAX {
183                let (quo, rem) = x.div_rem(&bigp);
184                if !rem.is_zero() {
185                    break;
186                }
187                x = quo;
188
189                count += 1;
190                if x.is_one() {
191                    primes.push((p, count));
192                    break 'outer;
193                }
194            }
195            if count != 0 {
196                primes.push((p, count));
197            }
198        }
199
200        let remainder = if x.is_one() && count2 == 0 {
201            None
202        } else {
203            Some(x << count2)
204        };
205        Self { primes, remainder }
206    }
207
208    fn decompose_known_factors(
209        mut primes: SmallVec<[(u16, u16); NUM_INLINE]>,
210        p: u16,
211        mut pcount: u16,
212        factor_indices: impl Iterator<Item = (u16, u16)>,
213    ) -> Self {
214        let mut remainder = 1_usize;
215        for (index, mut qcount) in factor_indices {
216            let index = index as usize;
217            let q = ODD_PRIMES[index - 1];
218
219            if pcount != 0 {
220                if p == q {
221                    qcount += pcount;
222                } else {
223                    primes.push((p, pcount));
224                }
225                pcount = 0;
226            }
227
228            if index < NUM_PRIMES {
229                primes.push((q, qcount));
230            } else {
231                let q = q as usize;
232                for _ in 0..qcount {
233                    remainder *= q;
234                }
235            }
236        }
237
238        let remainder = if remainder == 1 {
239            None
240        } else {
241            Some(remainder.into())
242        };
243        Self { primes, remainder }
244    }
245
246    fn decompose_u8(mut x: u8) -> Self {
247        let mut primes = SmallVec::new();
248
249        let count2 = x.trailing_zeros();
250        if count2 != 0 {
251            x >>= count2;
252            primes.push((2, count2 as u16));
253        }
254
255        // 8-bit integers always fit in the look-up table.
256        Self::decompose_known_factors(
257            primes,
258            0,
259            0,
260            known_odd_prime_factor_indices(x.into()).unwrap(),
261        )
262    }
263
264    fn decompose_u16(mut x: u16) -> Self {
265        const {
266            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
267        }
268
269        let mut primes = SmallVec::new();
270
271        let count2 = x.trailing_zeros();
272        if count2 != 0 {
273            x >>= count2;
274            primes.push((2, count2 as u16));
275        }
276
277        'outer: for (i, &p) in ODD_PRIMES.iter().enumerate().take(NUM_PRIMES - 1) {
278            let divider = OddDivider {
279                divisor: p,
280                multiplier: ODD_PRIME_DIVIDERS_U16[i],
281                shift: p.ilog2(),
282            };
283            let mut count = 0;
284            while x != 1 {
285                // Use look-up table as soon as the remainder is small enough.
286                if let Some(iter) = known_odd_prime_factor_indices(x.into()) {
287                    return Self::decompose_known_factors(primes, p, count, iter);
288                }
289
290                let (quo, rem) = divider.div_rem(x);
291                if rem != 0 {
292                    break;
293                }
294                x = quo;
295
296                count += 1;
297                if x == 1 {
298                    primes.push((p, count));
299                    break 'outer;
300                }
301            }
302            if count != 0 {
303                primes.push((p, count));
304            }
305        }
306
307        let remainder = if x == 1 { None } else { Some(x.into()) };
308        Self { primes, remainder }
309    }
310
311    fn decompose_u32(mut x: u32) -> Self {
312        const {
313            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
314        }
315
316        let mut primes = SmallVec::new();
317
318        let count2 = x.trailing_zeros();
319        if count2 != 0 {
320            x >>= count2;
321            primes.push((2, count2 as u16));
322        }
323
324        'outer: for (i, &p) in ODD_PRIMES.iter().enumerate().take(NUM_PRIMES - 1) {
325            let bigp = p as u32;
326            let divider = OddDivider {
327                divisor: bigp,
328                multiplier: ODD_PRIME_DIVIDERS_U32[i],
329                shift: p.ilog2(),
330            };
331            let mut count = 0;
332            while x != 1 {
333                // Use look-up table as soon as the remainder is small enough.
334                if let Ok(xx) = x.try_into()
335                    && let Some(iter) = known_odd_prime_factor_indices(xx)
336                {
337                    return Self::decompose_known_factors(primes, p, count, iter);
338                }
339
340                let (quo, rem) = divider.div_rem(x);
341                if rem != 0 {
342                    break;
343                }
344                x = quo;
345
346                count += 1;
347                if x == 1 {
348                    primes.push((p, count));
349                    break 'outer;
350                }
351            }
352            if count != 0 {
353                primes.push((p, count));
354            }
355        }
356
357        let remainder = if x == 1 { None } else { Some(x.into()) };
358        Self { primes, remainder }
359    }
360
361    fn decompose_u64(mut x: u64) -> Self {
362        const {
363            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
364        }
365
366        let mut primes = SmallVec::new();
367
368        let count2 = x.trailing_zeros();
369        if count2 != 0 {
370            x >>= count2;
371            primes.push((2, count2 as u16));
372        }
373
374        'outer: for (i, &p) in ODD_PRIMES.iter().enumerate().take(NUM_PRIMES - 1) {
375            let bigp = p as u64;
376            let divider = OddDivider {
377                divisor: bigp,
378                multiplier: ODD_PRIME_DIVIDERS_U64[i],
379                shift: p.ilog2(),
380            };
381            let mut count = 0;
382            while x != 1 {
383                // Use look-up table as soon as the remainder is small enough.
384                if let Ok(xx) = x.try_into()
385                    && let Some(iter) = known_odd_prime_factor_indices(xx)
386                {
387                    return Self::decompose_known_factors(primes, p, count, iter);
388                }
389
390                let (quo, rem) = divider.div_rem(x);
391                if rem != 0 {
392                    break;
393                }
394                x = quo;
395
396                count += 1;
397                if x == 1 {
398                    primes.push((p, count));
399                    break 'outer;
400                }
401            }
402            if count != 0 {
403                primes.push((p, count));
404            }
405        }
406
407        let remainder = if x == 1 { None } else { Some(x.into()) };
408        Self { primes, remainder }
409    }
410
411    fn decompose_u128(mut x: u128) -> Self {
412        const {
413            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
414        }
415
416        let mut primes = SmallVec::new();
417
418        let count2 = x.trailing_zeros();
419        if count2 != 0 {
420            x >>= count2;
421            primes.push((2, count2 as u16));
422        }
423
424        'outer: for &p in ODD_PRIMES.iter().take(NUM_PRIMES - 1) {
425            let bigp = p as u128;
426            let mut count = 0;
427            while x != 1 {
428                // Use look-up table as soon as the remainder is small enough.
429                if let Ok(xx) = x.try_into()
430                    && let Some(iter) = known_odd_prime_factor_indices(xx)
431                {
432                    return Self::decompose_known_factors(primes, p, count, iter);
433                }
434
435                if !x.is_multiple_of(bigp) {
436                    break;
437                }
438                x /= bigp;
439
440                count += 1;
441                if x == 1 {
442                    primes.push((p, count));
443                    break 'outer;
444                }
445            }
446            if count != 0 {
447                primes.push((p, count));
448            }
449        }
450
451        let remainder = if x == 1 { None } else { Some(x.into()) };
452        Self { primes, remainder }
453    }
454
455    fn decompose_usize(mut x: usize) -> Self {
456        const {
457            assert!(NUM_PRIMES <= ODD_PRIMES.len() + 1);
458        }
459
460        let mut primes = SmallVec::new();
461
462        let count2 = x.trailing_zeros();
463        if count2 != 0 {
464            x >>= count2;
465            primes.push((2, count2 as u16));
466        }
467
468        'outer: for &p in ODD_PRIMES.iter().take(NUM_PRIMES - 1) {
469            let bigp = p as usize;
470            let mut count = 0;
471            while x != 1 {
472                // Use look-up table as soon as the remainder is small enough.
473                if let Some(iter) = known_odd_prime_factor_indices(x) {
474                    return Self::decompose_known_factors(primes, p, count, iter);
475                }
476
477                if !x.is_multiple_of(bigp) {
478                    break;
479                }
480                x /= bigp;
481
482                count += 1;
483                if x == 1 {
484                    primes.push((p, count));
485                    break 'outer;
486                }
487            }
488            if count != 0 {
489                primes.push((p, count));
490            }
491        }
492
493        let remainder = if x == 1 { None } else { Some(x.into()) };
494        Self { primes, remainder }
495    }
496
497    /// Decomposes the given remainder using only primes whose bit mask is set.
498    fn decompose_mask(
499        remainder: &mut Option<BigUint>,
500        primes: &mut SmallVec<[(u16, u16); NUM_INLINE]>,
501        mask: &[u16],
502    ) {
503        let x: &mut BigUint = match remainder {
504            None => return,
505            Some(x) => x,
506        };
507
508        let mut new_primes = SmallVec::new();
509        'outer: for (p, (mut count, flag)) in Zip(
510            primes.iter().copied().peekable(),
511            mask.iter().copied().map(|p| (p, true)).peekable(),
512        ) {
513            if flag && count != u16::MAX {
514                let bigp = BigUint::from(p);
515                while count != u16::MAX {
516                    let (quo, rem) = x.div_rem(&bigp);
517                    if !rem.is_zero() {
518                        break;
519                    }
520                    *x = quo;
521
522                    count += 1;
523                    if x.is_one() {
524                        new_primes.push((p, count));
525                        break 'outer;
526                    }
527                }
528            }
529            new_primes.push((p, count));
530        }
531        *primes = new_primes;
532
533        if x.is_one() {
534            *remainder = None;
535        }
536    }
537
538    /// Computes `prime.pow(exponent)` and multiplies it into the accumulated
539    /// `(numerator, tmp)`.
540    fn accum_pow(numerator: &mut BigInt, tmp: &mut usize, prime: u16, exponent: u16) {
541        let prime = prime as usize;
542        for _ in 0..exponent {
543            match tmp.checked_mul(prime) {
544                Some(prod) => *tmp = prod,
545                None => {
546                    *numerator *= *tmp;
547                    *tmp = prime;
548                }
549            }
550        }
551    }
552
553    fn accum_pow_option(
554        remainder: &mut Option<BigUint>,
555        tmp: &mut usize,
556        prime: u16,
557        exponent: u16,
558    ) {
559        let prime = prime as usize;
560        for _ in 0..exponent {
561            match tmp.checked_mul(prime) {
562                Some(prod) => *tmp = prod,
563                None => {
564                    match remainder {
565                        None => *remainder = Some(BigUint::from(*tmp)),
566                        Some(r) => *r *= *tmp,
567                    };
568                    *tmp = prime;
569                }
570            }
571        }
572    }
573
574    fn pow_primes(primes: &mut [(u16, u16)], exponent: u32, remainder: &mut Option<BigUint>) {
575        for (p, i) in primes.iter_mut() {
576            let product = (*i as u32).strict_mul(exponent);
577            if product <= u16::MAX as u32 {
578                *i = product as u16;
579            } else {
580                *i = u16::MAX;
581                let factor = BigUint::from(*p).pow(product - u16::MAX as u32);
582                match remainder {
583                    None => *remainder = Some(factor),
584                    Some(r) => *r *= factor,
585                };
586            }
587        }
588    }
589
590    fn mul_primes(
591        lhs: &[(u16, u16)],
592        rhs: &[(u16, u16)],
593        remainder: &mut Option<BigUint>,
594    ) -> SmallVec<[(u16, u16); NUM_INLINE]> {
595        let mut primes = SmallVec::new();
596
597        for (p, (lcount, rcount)) in Zip(
598            lhs.iter().copied().peekable(),
599            rhs.iter().copied().peekable(),
600        ) {
601            let sum = lcount as u32 + rcount as u32;
602            if sum <= u16::MAX as u32 {
603                primes.push((p, sum as u16));
604            } else {
605                primes.push((p, u16::MAX));
606                let factor = BigUint::from(p).pow(sum - u16::MAX as u32);
607                match remainder {
608                    None => *remainder = Some(factor),
609                    Some(r) => *r *= factor,
610                };
611            }
612        }
613        primes
614    }
615
616    fn div_primes(
617        num: &[(u16, u16)],
618        denom: &[(u16, u16)],
619        mut remainder: Option<BigUint>,
620    ) -> (Self, SmallVec<[u16; NUM_INLINE]>) {
621        let mut primes = SmallVec::new();
622        let mut mask = SmallVec::new();
623        let mut tmp = 1_usize;
624
625        for (p, (num_count, denom_count)) in Zip(
626            num.iter().copied().peekable(),
627            denom.iter().copied().peekable(),
628        ) {
629            if num_count == u16::MAX {
630                mask.push(p);
631            }
632            match num_count.cmp(&denom_count) {
633                Ordering::Greater => {
634                    primes.push((p, num_count - denom_count));
635                }
636                Ordering::Equal => (),
637                Ordering::Less => {
638                    Self::accum_pow_option(&mut remainder, &mut tmp, p, denom_count - num_count);
639                    let factor = BigUint::from(p).pow(denom_count as u32 - num_count as u32);
640                    remainder = match remainder {
641                        None => Some(factor),
642                        Some(r) => Some(r * factor),
643                    };
644                }
645            }
646        }
647
648        if tmp != 1 {
649            match remainder.as_mut() {
650                None => remainder = Some(BigUint::from(tmp)),
651                Some(r) => *r *= tmp,
652            };
653        }
654        (Self { primes, remainder }, mask)
655    }
656}
657
658impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<u8>
659    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
660{
661    fn from(value: u8) -> Self {
662        if value == 0 {
663            panic!("Attempted to create a denominator of zero");
664        }
665        Self::decompose_u8(value)
666    }
667}
668
669impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<u16>
670    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
671{
672    fn from(value: u16) -> Self {
673        if value == 0 {
674            panic!("Attempted to create a denominator of zero");
675        }
676        Self::decompose_u16(value)
677    }
678}
679
680impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<u32>
681    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
682{
683    fn from(value: u32) -> Self {
684        if value == 0 {
685            panic!("Attempted to create a denominator of zero");
686        }
687        Self::decompose_u32(value)
688    }
689}
690
691impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<u64>
692    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
693{
694    fn from(value: u64) -> Self {
695        if value == 0 {
696            panic!("Attempted to create a denominator of zero");
697        }
698        Self::decompose_u64(value)
699    }
700}
701
702impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<u128>
703    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
704{
705    fn from(value: u128) -> Self {
706        if value == 0 {
707            panic!("Attempted to create a denominator of zero");
708        }
709        Self::decompose_u128(value)
710    }
711}
712
713impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<usize>
714    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
715{
716    fn from(value: usize) -> Self {
717        if value == 0 {
718            panic!("Attempted to create a denominator of zero");
719        }
720        Self::decompose_usize(value)
721    }
722}
723
724impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<BigUint>
725    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
726{
727    fn from(value: BigUint) -> Self {
728        if value.is_zero() {
729            panic!("Attempted to create a denominator of zero");
730        }
731        Self::decompose(value)
732    }
733}
734
735impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<&BigUint>
736    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
737{
738    fn from(value: &BigUint) -> Self {
739        if value.is_zero() {
740            panic!("Attempted to create a denominator of zero");
741        }
742        Self::decompose(value.clone())
743    }
744}
745
746impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
747    for BigUint
748{
749    fn from(value: DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> BigUint {
750        let mut result = match value.remainder {
751            Some(x) => x,
752            None => BigUint::ONE,
753        };
754        let mut tmp = 1_usize;
755        for (p, count) in value.primes.into_iter() {
756            let p = p as usize;
757            for _ in 0..count {
758                match tmp.checked_mul(p) {
759                    Some(prod) => tmp = prod,
760                    None => {
761                        result *= tmp;
762                        tmp = p;
763                    }
764                }
765            }
766        }
767        result * tmp
768    }
769}
770
771impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> From<&DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
772    for BigUint
773{
774    fn from(value: &DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> BigUint {
775        let mut result = match &value.remainder {
776            Some(x) => x.clone(),
777            None => BigUint::ONE,
778        };
779        let mut tmp = 1_usize;
780        for &(p, count) in value.primes.iter() {
781            let p = p as usize;
782            for _ in 0..count {
783                match tmp.checked_mul(p) {
784                    Some(prod) => tmp = prod,
785                    None => {
786                        result *= tmp;
787                        tmp = p;
788                    }
789                }
790            }
791        }
792        result * tmp
793    }
794}
795
796impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> One
797    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
798{
799    fn one() -> Self {
800        Self::ONE
801    }
802}
803
804impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> num_traits::Pow<u32>
805    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
806{
807    type Output = Self;
808
809    fn pow(mut self, rhs: u32) -> Self {
810        if rhs == 0 {
811            return Self::ONE;
812        }
813        Self::pow_primes(&mut self.primes, rhs, &mut self.remainder);
814        self
815    }
816}
817
818impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> num_traits::Pow<u32>
819    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
820{
821    type Output = DenomSparseU16<NUM_PRIMES, NUM_INLINE>;
822
823    fn pow(self, rhs: u32) -> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
824        if rhs == 0 {
825            return DenomSparseU16::ONE;
826        }
827        let mut primes = self.primes.clone();
828        let mut remainder = self.remainder.clone();
829        DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::pow_primes(&mut primes, rhs, &mut remainder);
830        DenomSparseU16 { primes, remainder }
831    }
832}
833
834impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul
835    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
836{
837    type Output = Self;
838
839    fn mul(self, rhs: Self) -> Self {
840        let mut remainder = match (self.remainder, rhs.remainder) {
841            (None, None) => None,
842            (None, Some(r)) => Some(r),
843            (Some(l), None) => Some(l),
844            (Some(l), Some(r)) => Some(l * r),
845        };
846        let primes = Self::mul_primes(&self.primes, &rhs.primes, &mut remainder);
847        Self { primes, remainder }
848    }
849}
850
851impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<&Self>
852    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
853{
854    type Output = Self;
855
856    fn mul(self, rhs: &Self) -> Self {
857        let mut remainder = match (self.remainder, &rhs.remainder) {
858            (None, None) => None,
859            (None, Some(r)) => Some(r.clone()),
860            (Some(l), None) => Some(l),
861            (Some(l), Some(r)) => Some(l * r),
862        };
863        let primes = Self::mul_primes(&self.primes, &rhs.primes, &mut remainder);
864        Self { primes, remainder }
865    }
866}
867
868impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul
869    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
870{
871    type Output = DenomSparseU16<NUM_PRIMES, NUM_INLINE>;
872
873    fn mul(self, rhs: Self) -> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
874        let mut remainder = match (&self.remainder, &rhs.remainder) {
875            (None, None) => None,
876            (None, Some(r)) => Some(r.clone()),
877            (Some(l), None) => Some(l.clone()),
878            (Some(l), Some(r)) => Some(l * r),
879        };
880        let primes = DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::mul_primes(
881            &self.primes,
882            &rhs.primes,
883            &mut remainder,
884        );
885        DenomSparseU16 { primes, remainder }
886    }
887}
888
889impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
890    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
891{
892    type Output = DenomSparseU16<NUM_PRIMES, NUM_INLINE>;
893
894    fn mul(
895        self,
896        rhs: DenomSparseU16<NUM_PRIMES, NUM_INLINE>,
897    ) -> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
898        let mut remainder = match (&self.remainder, rhs.remainder) {
899            (None, None) => None,
900            (None, Some(r)) => Some(r),
901            (Some(l), None) => Some(l.clone()),
902            (Some(l), Some(r)) => Some(l * r),
903        };
904        let primes = DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::mul_primes(
905            &self.primes,
906            &rhs.primes,
907            &mut remainder,
908        );
909        DenomSparseU16 { primes, remainder }
910    }
911}
912
913impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> MulAssign
914    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
915{
916    fn mul_assign(&mut self, rhs: Self) {
917        match (&mut self.remainder, rhs.remainder) {
918            (_, None) => (),
919            (None, Some(r)) => self.remainder = Some(r),
920            (Some(l), Some(r)) => *l *= r,
921        };
922        self.primes = Self::mul_primes(&self.primes, &rhs.primes, &mut self.remainder);
923    }
924}
925
926impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> MulAssign<&Self>
927    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
928{
929    fn mul_assign(&mut self, rhs: &Self) {
930        match (&mut self.remainder, &rhs.remainder) {
931            (_, None) => (),
932            (None, Some(r)) => self.remainder = Some(r.clone()),
933            (Some(l), Some(r)) => *l *= r,
934        };
935        self.primes = Self::mul_primes(&self.primes, &rhs.primes, &mut self.remainder);
936    }
937}
938
939impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Div
940    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
941{
942    type Output = Self;
943
944    /// Divides this denominator by the other one.
945    ///
946    /// This function panics if the other one doesn't divide this one.
947    fn div(self, rhs: Self) -> Self {
948        let (
949            Self {
950                mut primes,
951                remainder: rhs_remainder,
952            },
953            mask,
954        ) = Self::div_primes(&self.primes, &rhs.primes, rhs.remainder);
955
956        let mut remainder = match (self.remainder, rhs_remainder) {
957            (None, None) => None,
958            (None, Some(r)) => {
959                if !r.is_one() {
960                    panic!("Attempted to divide a denominator by a non-divisor");
961                }
962                None
963            }
964            (Some(l), None) => Some(l),
965            (Some(l), Some(r)) => {
966                let (quo, rem) = l.div_rem(&r);
967                if !rem.is_zero() {
968                    panic!("Attempted to divide a denominator by a non-divisor");
969                }
970                if quo.is_one() { None } else { Some(quo) }
971            }
972        };
973
974        Self::decompose_mask(&mut remainder, &mut primes, &mask);
975        Self { primes, remainder }
976    }
977}
978
979impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Div<&Self>
980    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
981{
982    type Output = Self;
983
984    /// Divides this denominator by the other one.
985    ///
986    /// This function panics if the other one doesn't divide this one.
987    fn div(self, rhs: &Self) -> Self {
988        let (
989            Self {
990                mut primes,
991                remainder: rhs_remainder,
992            },
993            mask,
994        ) = Self::div_primes(&self.primes, &rhs.primes, rhs.remainder.clone());
995
996        let mut remainder = match (self.remainder, rhs_remainder) {
997            (None, None) => None,
998            (None, Some(r)) => {
999                if !r.is_one() {
1000                    panic!("Attempted to divide a denominator by a non-divisor");
1001                }
1002                None
1003            }
1004            (Some(l), None) => Some(l),
1005            (Some(l), Some(r)) => {
1006                let (quo, rem) = l.div_rem(&r);
1007                if !rem.is_zero() {
1008                    panic!("Attempted to divide a denominator by a non-divisor");
1009                }
1010                if quo.is_one() { None } else { Some(quo) }
1011            }
1012        };
1013
1014        Self::decompose_mask(&mut remainder, &mut primes, &mask);
1015        Self { primes, remainder }
1016    }
1017}
1018
1019impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Div
1020    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1021{
1022    type Output = DenomSparseU16<NUM_PRIMES, NUM_INLINE>;
1023
1024    /// Divides this denominator by the other one.
1025    ///
1026    /// This function panics if the other one doesn't divide this one.
1027    fn div(self, rhs: Self) -> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
1028        let (
1029            DenomSparseU16 {
1030                mut primes,
1031                remainder: rhs_remainder,
1032            },
1033            mask,
1034        ) = DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::div_primes(
1035            &self.primes,
1036            &rhs.primes,
1037            rhs.remainder.clone(),
1038        );
1039
1040        let mut remainder = match (&self.remainder, rhs_remainder) {
1041            (None, None) => None,
1042            (None, Some(r)) => {
1043                if !r.is_one() {
1044                    panic!("Attempted to divide a denominator by a non-divisor");
1045                }
1046                None
1047            }
1048            (Some(l), None) => Some(l.clone()),
1049            (Some(l), Some(r)) => {
1050                let (quo, rem) = l.div_rem(&r);
1051                if !rem.is_zero() {
1052                    panic!("Attempted to divide a denominator by a non-divisor");
1053                }
1054                if quo.is_one() { None } else { Some(quo) }
1055            }
1056        };
1057
1058        DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::decompose_mask(
1059            &mut remainder,
1060            &mut primes,
1061            &mask,
1062        );
1063        DenomSparseU16 { primes, remainder }
1064    }
1065}
1066
1067impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Div<DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
1068    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1069{
1070    type Output = DenomSparseU16<NUM_PRIMES, NUM_INLINE>;
1071
1072    /// Divides this denominator by the other one.
1073    ///
1074    /// This function panics if the other one doesn't divide this one.
1075    fn div(
1076        self,
1077        rhs: DenomSparseU16<NUM_PRIMES, NUM_INLINE>,
1078    ) -> DenomSparseU16<NUM_PRIMES, NUM_INLINE> {
1079        let (
1080            DenomSparseU16 {
1081                mut primes,
1082                remainder: rhs_remainder,
1083            },
1084            mask,
1085        ) = DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::div_primes(
1086            &self.primes,
1087            &rhs.primes,
1088            rhs.remainder,
1089        );
1090
1091        let mut remainder = match (&self.remainder, rhs_remainder) {
1092            (None, None) => None,
1093            (None, Some(r)) => {
1094                if !r.is_one() {
1095                    panic!("Attempted to divide a denominator by a non-divisor");
1096                }
1097                None
1098            }
1099            (Some(l), None) => Some(l.clone()),
1100            (Some(l), Some(r)) => {
1101                let (quo, rem) = l.div_rem(&r);
1102                if !rem.is_zero() {
1103                    panic!("Attempted to divide a denominator by a non-divisor");
1104                }
1105                if quo.is_one() { None } else { Some(quo) }
1106            }
1107        };
1108
1109        DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::decompose_mask(
1110            &mut remainder,
1111            &mut primes,
1112            &mask,
1113        );
1114        DenomSparseU16 { primes, remainder }
1115    }
1116}
1117
1118impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> DivAssign
1119    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1120{
1121    /// Divides this denominator by the other one.
1122    ///
1123    /// This function panics if the other one doesn't divide this one.
1124    fn div_assign(&mut self, rhs: Self) {
1125        let (
1126            Self {
1127                primes,
1128                remainder: rhs_remainder,
1129            },
1130            mask,
1131        ) = Self::div_primes(&self.primes, &rhs.primes, rhs.remainder);
1132        self.primes = primes;
1133
1134        match (&mut self.remainder, rhs_remainder) {
1135            (_, None) => (),
1136            (None, Some(r)) => {
1137                if !r.is_one() {
1138                    panic!("Attempted to divide a denominator by a non-divisor");
1139                }
1140            }
1141            (Some(l), Some(r)) => {
1142                let (quo, rem) = l.div_rem(&r);
1143                if !rem.is_zero() {
1144                    panic!("Attempted to divide a denominator by a non-divisor");
1145                }
1146                if quo.is_one() {
1147                    self.remainder = None
1148                } else {
1149                    *l = quo;
1150                }
1151            }
1152        };
1153
1154        Self::decompose_mask(&mut self.remainder, &mut self.primes, &mask);
1155    }
1156}
1157
1158impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> DivAssign<&Self>
1159    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1160{
1161    /// Divides this denominator by the other one.
1162    ///
1163    /// This function panics if the other one doesn't divide this one.
1164    fn div_assign(&mut self, rhs: &Self) {
1165        let (
1166            Self {
1167                primes,
1168                remainder: rhs_remainder,
1169            },
1170            mask,
1171        ) = Self::div_primes(&self.primes, &rhs.primes, rhs.remainder.clone());
1172        self.primes = primes;
1173
1174        match (&mut self.remainder, rhs_remainder) {
1175            (_, None) => (),
1176            (None, Some(r)) => {
1177                if !r.is_one() {
1178                    panic!("Attempted to divide a denominator by a non-divisor");
1179                }
1180            }
1181            (Some(l), Some(r)) => {
1182                let (quo, rem) = l.div_rem(&r);
1183                if !rem.is_zero() {
1184                    panic!("Attempted to divide a denominator by a non-divisor");
1185                }
1186                if quo.is_one() {
1187                    self.remainder = None
1188                } else {
1189                    *l = quo;
1190                }
1191            }
1192        };
1193
1194        Self::decompose_mask(&mut self.remainder, &mut self.primes, &mask);
1195    }
1196}
1197
1198impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
1199    for BigInt
1200{
1201    type Output = Self;
1202
1203    fn mul(mut self, rhs: DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> Self {
1204        let mut tmp = 1_usize;
1205        for (p, count) in rhs.primes.into_iter() {
1206            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut self, &mut tmp, p, count);
1207        }
1208
1209        self *= tmp;
1210        if let Some(remainder) = rhs.remainder {
1211            self *= BigInt::from_biguint(Sign::Plus, remainder);
1212        }
1213        self
1214    }
1215}
1216
1217impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<&DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
1218    for BigInt
1219{
1220    type Output = Self;
1221
1222    fn mul(mut self, rhs: &DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> Self {
1223        let mut tmp = 1_usize;
1224        for &(p, count) in rhs.primes.iter() {
1225            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut self, &mut tmp, p, count);
1226        }
1227
1228        self *= tmp;
1229        if let Some(remainder) = &rhs.remainder {
1230            self *= BigInt::from_biguint(Sign::Plus, remainder.clone());
1231        }
1232        self
1233    }
1234}
1235
1236impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
1237    for &BigInt
1238{
1239    type Output = BigInt;
1240
1241    fn mul(self, rhs: DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> BigInt {
1242        let mut this = self.clone();
1243        let mut tmp = 1_usize;
1244        for (p, count) in rhs.primes.into_iter() {
1245            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut this, &mut tmp, p, count);
1246        }
1247
1248        this *= tmp;
1249        if let Some(remainder) = rhs.remainder {
1250            this *= BigInt::from_biguint(Sign::Plus, remainder);
1251        }
1252        this
1253    }
1254}
1255
1256impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<&DenomSparseU16<NUM_PRIMES, NUM_INLINE>>
1257    for &BigInt
1258{
1259    type Output = BigInt;
1260
1261    fn mul(self, rhs: &DenomSparseU16<NUM_PRIMES, NUM_INLINE>) -> BigInt {
1262        let mut this = self.clone();
1263        let mut tmp = 1_usize;
1264        for &(p, count) in rhs.primes.iter() {
1265            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut this, &mut tmp, p, count);
1266        }
1267
1268        this *= tmp;
1269        if let Some(remainder) = &rhs.remainder {
1270            this *= BigInt::from_biguint(Sign::Plus, remainder.clone());
1271        }
1272        this
1273    }
1274}
1275
1276impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<BigInt>
1277    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1278{
1279    type Output = BigInt;
1280
1281    fn mul(self, mut rhs: BigInt) -> BigInt {
1282        let mut tmp = 1_usize;
1283        for (p, count) in self.primes.into_iter() {
1284            Self::accum_pow(&mut rhs, &mut tmp, p, count);
1285        }
1286
1287        rhs *= tmp;
1288        if let Some(remainder) = self.remainder {
1289            rhs *= BigInt::from_biguint(Sign::Plus, remainder);
1290        }
1291        rhs
1292    }
1293}
1294
1295impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<&BigInt>
1296    for DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1297{
1298    type Output = BigInt;
1299
1300    fn mul(self, rhs: &BigInt) -> BigInt {
1301        let mut rhs = rhs.clone();
1302        let mut tmp = 1_usize;
1303        for (p, count) in self.primes.into_iter() {
1304            Self::accum_pow(&mut rhs, &mut tmp, p, count);
1305        }
1306
1307        rhs *= tmp;
1308        if let Some(remainder) = self.remainder {
1309            rhs *= BigInt::from_biguint(Sign::Plus, remainder);
1310        }
1311        rhs
1312    }
1313}
1314
1315impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<BigInt>
1316    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1317{
1318    type Output = BigInt;
1319
1320    fn mul(self, mut rhs: BigInt) -> BigInt {
1321        let mut tmp = 1_usize;
1322        for &(p, count) in self.primes.iter() {
1323            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut rhs, &mut tmp, p, count);
1324        }
1325
1326        rhs *= tmp;
1327        if let Some(remainder) = &self.remainder {
1328            rhs *= BigInt::from_biguint(Sign::Plus, remainder.clone());
1329        }
1330        rhs
1331    }
1332}
1333
1334impl<const NUM_PRIMES: usize, const NUM_INLINE: usize> Mul<&BigInt>
1335    for &DenomSparseU16<NUM_PRIMES, NUM_INLINE>
1336{
1337    type Output = BigInt;
1338
1339    fn mul(self, rhs: &BigInt) -> BigInt {
1340        let mut rhs = rhs.clone();
1341        let mut tmp = 1_usize;
1342        for &(p, count) in self.primes.iter() {
1343            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(&mut rhs, &mut tmp, p, count);
1344        }
1345
1346        rhs *= tmp;
1347        if let Some(remainder) = &self.remainder {
1348            rhs *= BigInt::from_biguint(Sign::Plus, remainder.clone());
1349        }
1350        rhs
1351    }
1352}
1353
1354impl<const NUM_PRIMES: usize, const NUM_INLINE: usize>
1355    MulAssign<DenomSparseU16<NUM_PRIMES, NUM_INLINE>> for BigInt
1356{
1357    fn mul_assign(&mut self, rhs: DenomSparseU16<NUM_PRIMES, NUM_INLINE>) {
1358        let mut tmp = 1_usize;
1359        for (p, count) in rhs.primes.into_iter() {
1360            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(self, &mut tmp, p, count);
1361        }
1362
1363        *self *= tmp;
1364        if let Some(remainder) = rhs.remainder {
1365            *self *= BigInt::from_biguint(Sign::Plus, remainder);
1366        }
1367    }
1368}
1369
1370impl<const NUM_PRIMES: usize, const NUM_INLINE: usize>
1371    MulAssign<&DenomSparseU16<NUM_PRIMES, NUM_INLINE>> for BigInt
1372{
1373    fn mul_assign(&mut self, rhs: &DenomSparseU16<NUM_PRIMES, NUM_INLINE>) {
1374        let mut tmp = 1_usize;
1375        for &(p, count) in rhs.primes.iter() {
1376            DenomSparseU16::<NUM_PRIMES, NUM_INLINE>::accum_pow(self, &mut tmp, p, count);
1377        }
1378
1379        *self *= tmp;
1380        if let Some(remainder) = &rhs.remainder {
1381            *self *= BigInt::from_biguint(Sign::Plus, remainder.clone());
1382        }
1383    }
1384}
1385
1386struct Zip<A, B>(A, B);
1387
1388impl<U, V, A, B> Iterator for Zip<Peekable<A>, Peekable<B>>
1389where
1390    A: Iterator<Item = (u16, U)>,
1391    B: Iterator<Item = (u16, V)>,
1392    U: Default,
1393    V: Default,
1394{
1395    type Item = (u16, (U, V));
1396
1397    fn next(&mut self) -> Option<Self::Item> {
1398        match (self.0.peek(), self.1.peek()) {
1399            (None, None) => None,
1400            (None, Some(_)) => {
1401                let (p, v) = self.1.next().unwrap();
1402                Some((p, (Default::default(), v)))
1403            }
1404            (Some(_), None) => {
1405                let (p, u) = self.0.next().unwrap();
1406                Some((p, (u, Default::default())))
1407            }
1408            (Some((p, _)), Some((q, _))) => match p.cmp(q) {
1409                Ordering::Less => {
1410                    let (p, u) = self.0.next().unwrap();
1411                    Some((p, (u, Default::default())))
1412                }
1413                Ordering::Greater => {
1414                    let (q, v) = self.1.next().unwrap();
1415                    Some((q, (Default::default(), v)))
1416                }
1417                Ordering::Equal => {
1418                    let (p, u) = self.0.next().unwrap();
1419                    let (_, v) = self.1.next().unwrap();
1420                    Some((p, (u, v)))
1421                }
1422            },
1423        }
1424    }
1425}
1426
1427#[cfg(test)]
1428mod tests {
1429    use super::*;
1430    use smallvec::smallvec;
1431
1432    macro_rules! tests {
1433        (
1434            $mod:ident,
1435            $num_primes:expr,
1436            $( $case:ident ,)*
1437        ) => {
1438            mod $mod {
1439                $(
1440                    #[test]
1441                    fn $case() {
1442                        super::$case::<$num_primes>();
1443                    }
1444                )*
1445            }
1446        };
1447    }
1448
1449    macro_rules! all_tests {
1450        (
1451            $mod:ident,
1452            $num_primes:expr
1453        ) => {
1454            tests!(
1455                $mod,
1456                $num_primes,
1457                test_decompose_to_biguint,
1458                test_decompose_small_prime_power,
1459                test_mul_prime_powers,
1460                test_div_prime_powers,
1461                test_product,
1462                test_normalize,
1463                test_gcd_reduce,
1464            );
1465        };
1466    }
1467
1468    all_tests!(denom24, 24);
1469    all_tests!(denom6542, 6542);
1470
1471    fn test_decompose_to_biguint<const NUM_PRIMES: usize>() {
1472        for i in 1_usize..=(1 << 20) {
1473            let bigi = BigUint::from(i);
1474            let x = DenomSparseU16::<NUM_PRIMES, 8>::from(&bigi);
1475            assert_eq!(x.to_biguint(), bigi);
1476        }
1477    }
1478
1479    #[test]
1480    fn test_decompose_no_remainder() {
1481        for i in 1_usize..=65536 {
1482            let bigi = BigUint::from(i);
1483            let x = DenomSparse6542::from(&bigi);
1484            assert_eq!(x.remainder, None, "Remainder for {i}");
1485        }
1486    }
1487
1488    #[test]
1489    fn test_decompose_known_values() {
1490        assert_eq!(
1491            DenomSparse6542::from(0xfff1_u16),
1492            DenomSparse6542 {
1493                primes: smallvec![(0xfff1, 1)],
1494                remainder: None,
1495            }
1496        );
1497        assert_eq!(
1498            DenomSparse6542::from(0xffff_fffb_u32),
1499            DenomSparse6542 {
1500                primes: smallvec![],
1501                remainder: Some(0xffff_fffb_u32.into()),
1502            }
1503        );
1504        assert_eq!(
1505            DenomSparse6542::from(0xffff_ffff_ffff_ffc5_u64),
1506            DenomSparse6542 {
1507                primes: smallvec![],
1508                remainder: Some(0xffff_ffff_ffff_ffc5_u64.into()),
1509            }
1510        );
1511        assert_eq!(
1512            DenomSparse6542::from(0xffff_ffff_ffff_ffff_ffff_ffff_ffff_ff61_u128),
1513            DenomSparse6542 {
1514                primes: smallvec![],
1515                remainder: Some(0xffff_ffff_ffff_ffff_ffff_ffff_ffff_ff61_u128.into()),
1516            }
1517        );
1518        assert_eq!(
1519            DenomSparse6542::from(BigUint::from_slice(&[
1520                0xffff_ff43,
1521                0xffff_ffff,
1522                0xffff_ffff,
1523                0xffff_ffff,
1524                0xffff_ffff,
1525                0xffff_ffff,
1526                0xffff_ffff,
1527                0xffff_ffff,
1528            ])),
1529            DenomSparse6542 {
1530                primes: smallvec![],
1531                remainder: Some(BigUint::from_slice(&[
1532                    0xffff_ff43,
1533                    0xffff_ffff,
1534                    0xffff_ffff,
1535                    0xffff_ffff,
1536                    0xffff_ffff,
1537                    0xffff_ffff,
1538                    0xffff_ffff,
1539                    0xffff_ffff,
1540                ])),
1541            }
1542        );
1543
1544        assert_eq!(
1545            DenomSparse6542::from(2u8 * 3 * 5 * 7),
1546            DenomSparse6542 {
1547                primes: smallvec![(2, 1), (3, 1), (5, 1), (7, 1)],
1548                remainder: None,
1549            }
1550        );
1551        assert_eq!(
1552            DenomSparse6542::from(2u16 * 3 * 5 * 7 * 11 * 13),
1553            DenomSparse6542 {
1554                primes: smallvec![(2, 1), (3, 1), (5, 1), (7, 1), (11, 1), (13, 1)],
1555                remainder: None,
1556            }
1557        );
1558        assert_eq!(
1559            DenomSparse6542::from(2u32 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23),
1560            DenomSparse6542 {
1561                primes: smallvec![
1562                    (2, 1),
1563                    (3, 1),
1564                    (5, 1),
1565                    (7, 1),
1566                    (11, 1),
1567                    (13, 1),
1568                    (17, 1),
1569                    (19, 1),
1570                    (23, 1)
1571                ],
1572                remainder: None,
1573            }
1574        );
1575        assert_eq!(
1576            DenomSparse6542::from(
1577                2u64 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 31 * 37 * 41 * 43 * 47 * 53
1578            ),
1579            DenomSparse6542 {
1580                primes: smallvec![
1581                    (2, 1),
1582                    (3, 1),
1583                    (5, 1),
1584                    (7, 1),
1585                    (11, 1),
1586                    (13, 1),
1587                    (17, 1),
1588                    (19, 1),
1589                    (23, 1),
1590                    (31, 1),
1591                    (37, 1),
1592                    (41, 1),
1593                    (43, 1),
1594                    (47, 1),
1595                    (53, 1)
1596                ],
1597                remainder: None,
1598            }
1599        );
1600        assert_eq!(
1601            DenomSparse6542::from(
1602                2u128
1603                    * 3
1604                    * 5
1605                    * 7
1606                    * 11
1607                    * 13
1608                    * 17
1609                    * 19
1610                    * 23
1611                    * 31
1612                    * 37
1613                    * 41
1614                    * 43
1615                    * 47
1616                    * 53
1617                    * 59
1618                    * 61
1619                    * 67
1620                    * 71
1621                    * 73
1622                    * 79
1623                    * 83
1624                    * 89
1625                    * 97
1626                    * 101,
1627            ),
1628            DenomSparse6542 {
1629                primes: smallvec![
1630                    (2, 1),
1631                    (3, 1),
1632                    (5, 1),
1633                    (7, 1),
1634                    (11, 1),
1635                    (13, 1),
1636                    (17, 1),
1637                    (19, 1),
1638                    (23, 1),
1639                    (31, 1),
1640                    (37, 1),
1641                    (41, 1),
1642                    (43, 1),
1643                    (47, 1),
1644                    (53, 1),
1645                    (59, 1),
1646                    (61, 1),
1647                    (67, 1),
1648                    (71, 1),
1649                    (73, 1),
1650                    (79, 1),
1651                    (83, 1),
1652                    (89, 1),
1653                    (97, 1),
1654                    (101, 1),
1655                ],
1656                remainder: None,
1657            }
1658        );
1659        assert_eq!(
1660            DenomSparse6542::from(
1661                BigUint::try_from(
1662                    BigInt::from(2)
1663                        * 3
1664                        * 5
1665                        * 7
1666                        * 11
1667                        * 13
1668                        * 17
1669                        * 19
1670                        * 23
1671                        * 31
1672                        * 37
1673                        * 41
1674                        * 43
1675                        * 47
1676                        * 53
1677                        * 59
1678                        * 61
1679                        * 67
1680                        * 71
1681                        * 73
1682                        * 79
1683                        * 83
1684                        * 89
1685                        * 97
1686                        * 101
1687                        * 103
1688                        * 107
1689                        * 109
1690                        * 113
1691                        * 127
1692                        * 131
1693                        * 137
1694                        * 139
1695                        * 149
1696                        * 151
1697                        * 157
1698                        * 163
1699                        * 167
1700                        * 173
1701                        * 179
1702                        * 181
1703                        * 191
1704                        * 193
1705                        * 197
1706                        * 199
1707                        * 211
1708                        * 223
1709                        * 227
1710                        * 229
1711                        * 233
1712                )
1713                .unwrap(),
1714            ),
1715            DenomSparse6542 {
1716                primes: smallvec![
1717                    (2, 1),
1718                    (3, 1),
1719                    (5, 1),
1720                    (7, 1),
1721                    (11, 1),
1722                    (13, 1),
1723                    (17, 1),
1724                    (19, 1),
1725                    (23, 1),
1726                    (31, 1),
1727                    (37, 1),
1728                    (41, 1),
1729                    (43, 1),
1730                    (47, 1),
1731                    (53, 1),
1732                    (59, 1),
1733                    (61, 1),
1734                    (67, 1),
1735                    (71, 1),
1736                    (73, 1),
1737                    (79, 1),
1738                    (83, 1),
1739                    (89, 1),
1740                    (97, 1),
1741                    (101, 1),
1742                    (103, 1),
1743                    (107, 1),
1744                    (109, 1),
1745                    (113, 1),
1746                    (127, 1),
1747                    (131, 1),
1748                    (137, 1),
1749                    (139, 1),
1750                    (149, 1),
1751                    (151, 1),
1752                    (157, 1),
1753                    (163, 1),
1754                    (167, 1),
1755                    (173, 1),
1756                    (179, 1),
1757                    (181, 1),
1758                    (191, 1),
1759                    (193, 1),
1760                    (197, 1),
1761                    (199, 1),
1762                    (211, 1),
1763                    (223, 1),
1764                    (227, 1),
1765                    (229, 1),
1766                    (233, 1),
1767                ],
1768                remainder: None,
1769            }
1770        );
1771
1772        assert_eq!(
1773            DenomSparse6542::from(0xfb_u16 * 0xf1),
1774            DenomSparse6542 {
1775                primes: smallvec![(0xf1, 1), (0xfb, 1)],
1776                remainder: None,
1777            }
1778        );
1779        assert_eq!(
1780            DenomSparse6542::from(0xfff1_u32 * 0xffef),
1781            DenomSparse6542 {
1782                primes: smallvec![(0xffef, 1), (0xfff1, 1)],
1783                remainder: None,
1784            }
1785        );
1786        assert_eq!(
1787            DenomSparse6542::from(0xfff1_u64 * 0xffef * 0xffd9 * 0xffc7),
1788            DenomSparse6542 {
1789                primes: smallvec![(0xffc7, 1), (0xffd9, 1), (0xffef, 1), (0xfff1, 1)],
1790                remainder: None,
1791            }
1792        );
1793        assert_eq!(
1794            DenomSparse6542::from(
1795                0xfff1_u128 * 0xffef * 0xffd9 * 0xffc7 * 0xffa9 * 0xffa7 * 0xff9d * 0xff8f,
1796            ),
1797            DenomSparse6542 {
1798                primes: smallvec![
1799                    (0xff8f, 1),
1800                    (0xff9d, 1),
1801                    (0xffa7, 1),
1802                    (0xffa9, 1),
1803                    (0xffc7, 1),
1804                    (0xffd9, 1),
1805                    (0xffef, 1),
1806                    (0xfff1, 1)
1807                ],
1808                remainder: None,
1809            }
1810        );
1811        assert_eq!(
1812            DenomSparse6542::from(
1813                BigUint::from(0xfff1_u16)
1814                    * 0xffef_u16
1815                    * 0xffd9_u16
1816                    * 0xffc7_u16
1817                    * 0xffa9_u16
1818                    * 0xffa7_u16
1819                    * 0xff9d_u16
1820                    * 0xff8f_u16
1821                    * 0xff8b_u16
1822                    * 0xff85_u16
1823                    * 0xff7f_u16
1824                    * 0xff71_u16
1825                    * 0xff65_u16
1826                    * 0xff5b_u16
1827                    * 0xff4d_u16
1828                    * 0xff49_u16,
1829            ),
1830            DenomSparse6542 {
1831                primes: smallvec![
1832                    (0xff49, 1),
1833                    (0xff4d, 1),
1834                    (0xff5b, 1),
1835                    (0xff65, 1),
1836                    (0xff71, 1),
1837                    (0xff7f, 1),
1838                    (0xff85, 1),
1839                    (0xff8b, 1),
1840                    (0xff8f, 1),
1841                    (0xff9d, 1),
1842                    (0xffa7, 1),
1843                    (0xffa9, 1),
1844                    (0xffc7, 1),
1845                    (0xffd9, 1),
1846                    (0xffef, 1),
1847                    (0xfff1, 1)
1848                ],
1849                remainder: None,
1850            }
1851        );
1852
1853        assert_eq!(
1854            DenomSparse6542::from(BigUint::from(128_usize)),
1855            DenomSparse6542 {
1856                primes: smallvec![(2, 7)],
1857                remainder: None,
1858            }
1859        );
1860        assert_eq!(
1861            DenomSparse6542::from(BigUint::from(89_usize)),
1862            DenomSparse6542 {
1863                primes: smallvec![(89, 1)],
1864                remainder: None,
1865            }
1866        );
1867        assert_eq!(
1868            DenomSparse6542::from(BigUint::from(97_usize)),
1869            DenomSparse6542 {
1870                primes: smallvec![(97, 1)],
1871                remainder: None,
1872            }
1873        );
1874        assert_eq!(
1875            DenomSparse6542::from(BigUint::from(97000_usize)),
1876            DenomSparse6542 {
1877                primes: smallvec![(2, 3), (5, 3), (97, 1)],
1878                remainder: None,
1879            }
1880        );
1881    }
1882
1883    fn test_decompose_small_prime_power<const NUM_PRIMES: usize>() {
1884        for p in std::iter::once(2).chain(ODD_PRIMES).take(NUM_PRIMES) {
1885            let bigp = BigUint::from(p);
1886            assert_eq!(
1887                DenomSparseU16::<NUM_PRIMES, 8>::from(bigp.pow(100)),
1888                DenomSparseU16 {
1889                    primes: smallvec![(p, 100)],
1890                    remainder: None,
1891                }
1892            );
1893        }
1894    }
1895
1896    fn test_mul_prime_powers<const NUM_PRIMES: usize>() {
1897        let p = BigUint::from(2u32);
1898        for a in 1..=256 {
1899            let denom_a = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(a));
1900            for b in 1..=256 {
1901                let denom_b = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(b));
1902                let denom_ab = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(a + b));
1903                assert_eq!(&denom_a * denom_b, denom_ab);
1904            }
1905        }
1906    }
1907
1908    fn test_div_prime_powers<const NUM_PRIMES: usize>() {
1909        let p = BigUint::from(2u32);
1910        for a in 1..=256 {
1911            let denom_a = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(a));
1912            for b in 1..=256 {
1913                let denom_b = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(b));
1914                let denom_ab = DenomSparseU16::<NUM_PRIMES, 8>::from(p.pow(a + b));
1915                assert_eq!(denom_ab / denom_b, denom_a);
1916            }
1917        }
1918    }
1919
1920    fn test_product<const NUM_PRIMES: usize>() {
1921        let values = (100..200)
1922            .map(|i: usize| DenomSparseU16::<NUM_PRIMES, 8>::from(BigUint::from(i)))
1923            .collect::<Vec<_>>();
1924        for (i, x) in values.iter().enumerate().map(|(i, x)| (i + 100, x)) {
1925            for (j, y) in values.iter().enumerate().map(|(j, y)| (j + 100, y)) {
1926                let z = x * y;
1927                assert_eq!(
1928                    z,
1929                    DenomSparseU16::<NUM_PRIMES, 8>::from(BigUint::from(i * j))
1930                );
1931
1932                for (p, (zcount, (xcount, ycount))) in Zip(
1933                    z.primes.into_iter().peekable(),
1934                    Zip(
1935                        x.primes.iter().copied().peekable(),
1936                        y.primes.iter().copied().peekable(),
1937                    )
1938                    .peekable(),
1939                ) {
1940                    assert_eq!(
1941                        zcount,
1942                        xcount + ycount,
1943                        "{zcount} != {xcount} + {ycount} for {p}"
1944                    );
1945                }
1946            }
1947        }
1948    }
1949
1950    fn test_normalize<const NUM_PRIMES: usize>() {
1951        let values = (100..200)
1952            .map(|i: usize| DenomSparseU16::<NUM_PRIMES, 8>::from(BigUint::from(i)))
1953            .collect::<Vec<_>>();
1954        for x in &values {
1955            for y in &values {
1956                let mut xnum = BigInt::one();
1957                let mut ynum = BigInt::one();
1958                let lcm = DenomSparseU16::<NUM_PRIMES, 8>::normalize(&mut xnum, &mut ynum, x, y);
1959                let lcm_bigint = lcm.to_biguint();
1960                let xnum = xnum.to_biguint().unwrap();
1961                let ynum = ynum.to_biguint().unwrap();
1962
1963                assert_eq!(xnum * x.to_biguint(), lcm_bigint);
1964                assert_eq!(ynum * y.to_biguint(), lcm_bigint);
1965
1966                for (p, (lcm_count, (xcount, ycount))) in Zip(
1967                    lcm.primes.into_iter().peekable(),
1968                    Zip(
1969                        x.primes.iter().copied().peekable(),
1970                        y.primes.iter().copied().peekable(),
1971                    )
1972                    .peekable(),
1973                ) {
1974                    assert_eq!(
1975                        lcm_count,
1976                        std::cmp::max(xcount, ycount),
1977                        "{lcm_count} != max({xcount}, {ycount}) for {p}"
1978                    );
1979                }
1980            }
1981        }
1982    }
1983
1984    fn test_gcd_reduce<const NUM_PRIMES: usize>() {
1985        let mut num = BigInt::from(-3 * 97);
1986        let mut denom = DenomSparseU16::<NUM_PRIMES, 8>::from(BigUint::from(3u32 * 5 * 97));
1987        denom.gcd_reduce(&mut num);
1988        assert_eq!(num, BigInt::from(-1));
1989        assert_eq!(
1990            denom,
1991            DenomSparseU16::<NUM_PRIMES, 8>::from(BigUint::from(5u32))
1992        );
1993    }
1994}
smart_big_rational/denom_sparse.rs

smart_big_rational/
denom_sparse.rs
