stv_rs/arithmetic/
approx_rational.rs

1// Copyright 2023 Google LLC
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 super::{Rational, RationalRef};
16use num::traits::{One, Zero};
17use num::{BigInt, BigRational};
18use std::cmp::Ordering;
19use std::fmt::{Debug, Display};
20use std::iter::{Product, Sum};
21use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign};
22
23#[derive(Clone, Debug, PartialEq)]
24struct Denom {
25    primes: [u8; Self::NUM_PRIMES],
26    remainder: Option<BigInt>,
27}
28
29impl Denom {
30    const NUM_PRIMES: usize = 24;
31    const PRIMES: [u64; Self::NUM_PRIMES] = [
32        2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
33    ];
34
35    fn one() -> Self {
36        Denom {
37            primes: [0; Self::NUM_PRIMES],
38            remainder: None,
39        }
40    }
41
42    fn precision() -> Self {
43        let mut primes = [0; Denom::NUM_PRIMES];
44        // 1_000_000 = (2*5)^6
45        primes[0] = 6; // prime = 2
46        primes[2] = 6; // prime = 5
47        Denom {
48            primes,
49            remainder: None,
50        }
51    }
52
53    fn to_bigint(&self) -> BigInt {
54        let mut result = match &self.remainder {
55            Some(x) => x.clone(),
56            None => BigInt::one(),
57        };
58        let mut tmp = 1u64;
59        for (i, &count) in self.primes.iter().enumerate() {
60            let p = Self::PRIMES[i];
61            for _ in 0..count {
62                match tmp.checked_mul(p) {
63                    Some(prod) => tmp = prod,
64                    None => {
65                        result *= tmp;
66                        tmp = p;
67                    }
68                }
69            }
70        }
71        result * tmp
72    }
73
74    const fn decompose_small(mut x: u64) -> Self {
75        let mut primes = [0; Self::NUM_PRIMES];
76        // TODO: use a for loop once supported in `const fn` context.
77        let mut i = 0;
78        while x > 1 && i < Self::NUM_PRIMES {
79            let p = Self::PRIMES[i];
80            while x % p == 0 {
81                x /= p;
82                primes[i] += 1;
83            }
84            i += 1;
85        }
86
87        if x != 1 {
88            panic!("Failed to decompose small integer into small prime factors.");
89        }
90        Denom {
91            primes,
92            remainder: None,
93        }
94    }
95
96    // TODO: Use std::array::from_fn when available in const contexts.
97    const DECOMPOSED: [Self; 90] = [
98        Self::decompose_small(1),
99        Self::decompose_small(2),
100        Self::decompose_small(3),
101        Self::decompose_small(4),
102        Self::decompose_small(5),
103        Self::decompose_small(6),
104        Self::decompose_small(7),
105        Self::decompose_small(8),
106        Self::decompose_small(9),
107        Self::decompose_small(10),
108        Self::decompose_small(11),
109        Self::decompose_small(12),
110        Self::decompose_small(13),
111        Self::decompose_small(14),
112        Self::decompose_small(15),
113        Self::decompose_small(16),
114        Self::decompose_small(17),
115        Self::decompose_small(18),
116        Self::decompose_small(19),
117        Self::decompose_small(20),
118        Self::decompose_small(21),
119        Self::decompose_small(22),
120        Self::decompose_small(23),
121        Self::decompose_small(24),
122        Self::decompose_small(25),
123        Self::decompose_small(26),
124        Self::decompose_small(27),
125        Self::decompose_small(28),
126        Self::decompose_small(29),
127        Self::decompose_small(30),
128        Self::decompose_small(31),
129        Self::decompose_small(32),
130        Self::decompose_small(33),
131        Self::decompose_small(34),
132        Self::decompose_small(35),
133        Self::decompose_small(36),
134        Self::decompose_small(37),
135        Self::decompose_small(38),
136        Self::decompose_small(39),
137        Self::decompose_small(40),
138        Self::decompose_small(41),
139        Self::decompose_small(42),
140        Self::decompose_small(43),
141        Self::decompose_small(44),
142        Self::decompose_small(45),
143        Self::decompose_small(46),
144        Self::decompose_small(47),
145        Self::decompose_small(48),
146        Self::decompose_small(49),
147        Self::decompose_small(50),
148        Self::decompose_small(51),
149        Self::decompose_small(52),
150        Self::decompose_small(53),
151        Self::decompose_small(54),
152        Self::decompose_small(55),
153        Self::decompose_small(56),
154        Self::decompose_small(57),
155        Self::decompose_small(58),
156        Self::decompose_small(59),
157        Self::decompose_small(60),
158        Self::decompose_small(61),
159        Self::decompose_small(62),
160        Self::decompose_small(63),
161        Self::decompose_small(64),
162        Self::decompose_small(65),
163        Self::decompose_small(66),
164        Self::decompose_small(67),
165        Self::decompose_small(68),
166        Self::decompose_small(69),
167        Self::decompose_small(70),
168        Self::decompose_small(71),
169        Self::decompose_small(72),
170        Self::decompose_small(73),
171        Self::decompose_small(74),
172        Self::decompose_small(75),
173        Self::decompose_small(76),
174        Self::decompose_small(77),
175        Self::decompose_small(78),
176        Self::decompose_small(79),
177        Self::decompose_small(80),
178        Self::decompose_small(81),
179        Self::decompose_small(82),
180        Self::decompose_small(83),
181        Self::decompose_small(84),
182        Self::decompose_small(85),
183        Self::decompose_small(86),
184        Self::decompose_small(87),
185        Self::decompose_small(88),
186        Self::decompose_small(89),
187        Self::decompose_small(90),
188    ];
189
190    fn decompose(x: BigInt) -> Self {
191        if x > BigInt::zero() && x <= BigInt::from(90) {
192            return Self::DECOMPOSED[TryInto::<usize>::try_into(x).unwrap() - 1].clone();
193        }
194        Self::decompose_now(x)
195    }
196
197    fn decompose_now(mut x: BigInt) -> Self {
198        let mut primes = [0; Self::NUM_PRIMES];
199        'outer: for (i, &p) in Self::PRIMES.iter().enumerate() {
200            while (&x % p).is_zero() {
201                x /= p;
202                primes[i] += 1;
203                if x.is_one() {
204                    break 'outer;
205                }
206            }
207        }
208
209        if x.is_one() {
210            Denom {
211                primes,
212                remainder: None,
213            }
214        } else {
215            Denom {
216                primes,
217                remainder: Some(x),
218            }
219        }
220    }
221
222    /// Returns the least common multiple of two denominators, adjusting the
223    /// numerators accordingly.
224    fn normalize(lnum: &mut BigInt, rnum: &mut BigInt, ldenom: &Self, rdenom: &Self) -> Self {
225        let mut primes = [0; Self::NUM_PRIMES];
226        let mut ltmp = 1u64;
227        let mut rtmp = 1u64;
228        for (i, &p) in Self::PRIMES.iter().enumerate() {
229            let lcount = ldenom.primes[i];
230            let rcount = rdenom.primes[i];
231            match lcount.cmp(&rcount) {
232                Ordering::Equal => {
233                    primes[i] = lcount;
234                }
235                Ordering::Less => {
236                    Self::accum_pow(lnum, &mut ltmp, p, rcount - lcount);
237                    primes[i] = rcount;
238                }
239                Ordering::Greater => {
240                    Self::accum_pow(rnum, &mut rtmp, p, lcount - rcount);
241                    primes[i] = lcount;
242                }
243            }
244        }
245
246        *lnum *= ltmp;
247        *rnum *= rtmp;
248        let remainder = match (&ldenom.remainder, &rdenom.remainder) {
249            (None, None) => None,
250            (None, Some(r)) => {
251                *lnum *= r;
252                Some(r.clone())
253            }
254            (Some(l), None) => {
255                *rnum *= l;
256                Some(l.clone())
257            }
258            (Some(l), Some(r)) => {
259                if l == r {
260                    Some(l.clone())
261                } else {
262                    *lnum *= r;
263                    *rnum *= l;
264                    Some(l * r)
265                }
266            }
267        };
268        Denom { primes, remainder }
269    }
270
271    /// Computes `prime.pow(exponent)` and multiplies it into the accumulated
272    /// `(numerator, tmp)`.
273    fn accum_pow(numerator: &mut BigInt, tmp: &mut u64, prime: u64, exponent: u8) {
274        for _ in 0..exponent {
275            match tmp.checked_mul(prime) {
276                Some(prod) => *tmp = prod,
277                None => {
278                    *numerator *= *tmp;
279                    *tmp = prime;
280                }
281            }
282        }
283    }
284}
285
286impl Mul<&'_ Denom> for &'_ Denom {
287    type Output = Denom;
288    #[allow(clippy::needless_range_loop)]
289    fn mul(self, rhs: &'_ Denom) -> Denom {
290        let mut primes = [0; Denom::NUM_PRIMES];
291        for i in 0..Denom::NUM_PRIMES {
292            primes[i] = self.primes[i].checked_add(rhs.primes[i]).unwrap();
293        }
294        let remainder = match (&self.remainder, &rhs.remainder) {
295            (None, None) => None,
296            (None, Some(r)) => Some(r.clone()),
297            (Some(l), None) => Some(l.clone()),
298            (Some(l), Some(r)) => Some(l * r),
299        };
300        Denom { primes, remainder }
301    }
302}
303
304/// A [`BigRational`] that approximates to some precision in
305/// [`Rational::div_up_as_keep_factor()`]. The other operations behave exactly
306/// as [`BigRational`].
307#[derive(Clone, Debug)]
308pub struct ApproxRational {
309    num: BigInt,
310    denom: Denom,
311}
312
313impl PartialEq for ApproxRational {
314    fn eq(&self, rhs: &Self) -> bool {
315        &self.num * rhs.denom.to_bigint() == &rhs.num * self.denom.to_bigint()
316    }
317}
318
319impl Eq for ApproxRational {}
320
321impl PartialOrd for ApproxRational {
322    fn partial_cmp(&self, rhs: &Self) -> Option<Ordering> {
323        Some(self.cmp(rhs))
324    }
325}
326
327impl Ord for ApproxRational {
328    fn cmp(&self, rhs: &Self) -> Ordering {
329        (&self.num * rhs.denom.to_bigint()).cmp(&(&rhs.num * self.denom.to_bigint()))
330    }
331}
332
333impl ApproxRational {
334    fn reduced(&self) -> BigRational {
335        BigRational::new(self.num.clone(), self.denom.to_bigint())
336    }
337}
338
339impl Display for ApproxRational {
340    fn fmt(&self, f: &mut std::fmt::Formatter) -> Result<(), std::fmt::Error> {
341        Display::fmt(&self.reduced(), f)
342    }
343}
344
345impl Zero for ApproxRational {
346    fn zero() -> Self {
347        ApproxRational {
348            num: BigInt::zero(),
349            denom: Denom::one(),
350        }
351    }
352    fn is_zero(&self) -> bool {
353        self.num.is_zero()
354    }
355}
356impl One for ApproxRational {
357    fn one() -> Self {
358        ApproxRational {
359            num: BigInt::one(),
360            denom: Denom::one(),
361        }
362    }
363}
364
365impl Add for ApproxRational {
366    type Output = Self;
367    fn add(mut self, mut rhs: Self) -> Self {
368        let denom = Denom::normalize(&mut self.num, &mut rhs.num, &self.denom, &rhs.denom);
369        ApproxRational {
370            num: self.num + rhs.num,
371            denom,
372        }
373    }
374}
375impl Sub for ApproxRational {
376    type Output = Self;
377    fn sub(mut self, mut rhs: Self) -> Self {
378        let denom = Denom::normalize(&mut self.num, &mut rhs.num, &self.denom, &rhs.denom);
379        ApproxRational {
380            num: self.num - rhs.num,
381            denom,
382        }
383    }
384}
385impl Mul for ApproxRational {
386    type Output = Self;
387    fn mul(self, rhs: Self) -> Self {
388        ApproxRational {
389            num: self.num * rhs.num,
390            denom: &self.denom * &rhs.denom,
391        }
392    }
393}
394impl Mul<BigInt> for ApproxRational {
395    type Output = Self;
396    fn mul(self, rhs: BigInt) -> Self {
397        ApproxRational {
398            num: self.num * rhs,
399            denom: self.denom,
400        }
401    }
402}
403impl Div<BigInt> for ApproxRational {
404    type Output = Self;
405    #[allow(clippy::suspicious_arithmetic_impl)]
406    fn div(self, rhs: BigInt) -> Self {
407        ApproxRational {
408            num: self.num,
409            denom: &self.denom * &Denom::decompose(rhs),
410        }
411    }
412}
413
414impl Add<&'_ Self> for ApproxRational {
415    type Output = Self;
416    fn add(mut self, rhs: &'_ Self) -> Self {
417        let mut rnum = rhs.num.clone();
418        let denom = Denom::normalize(&mut self.num, &mut rnum, &self.denom, &rhs.denom);
419        ApproxRational {
420            num: self.num + rnum,
421            denom,
422        }
423    }
424}
425impl Sub<&'_ Self> for ApproxRational {
426    type Output = Self;
427    fn sub(mut self, rhs: &'_ Self) -> Self {
428        let mut rnum = rhs.num.clone();
429        let denom = Denom::normalize(&mut self.num, &mut rnum, &self.denom, &rhs.denom);
430        ApproxRational {
431            num: self.num - rnum,
432            denom,
433        }
434    }
435}
436impl Mul<&'_ Self> for ApproxRational {
437    type Output = Self;
438    fn mul(self, rhs: &'_ Self) -> Self {
439        ApproxRational {
440            num: self.num * &rhs.num,
441            denom: &self.denom * &rhs.denom,
442        }
443    }
444}
445
446impl Add<&'_ ApproxRational> for &'_ ApproxRational {
447    type Output = ApproxRational;
448    fn add(self, rhs: &'_ ApproxRational) -> ApproxRational {
449        let mut lnum = self.num.clone();
450        let mut rnum = rhs.num.clone();
451        let denom = Denom::normalize(&mut lnum, &mut rnum, &self.denom, &rhs.denom);
452        ApproxRational {
453            num: lnum + rnum,
454            denom,
455        }
456    }
457}
458impl Sub<&'_ ApproxRational> for &'_ ApproxRational {
459    type Output = ApproxRational;
460    fn sub(self, rhs: &'_ ApproxRational) -> ApproxRational {
461        let mut lnum = self.num.clone();
462        let mut rnum = rhs.num.clone();
463        let denom = Denom::normalize(&mut lnum, &mut rnum, &self.denom, &rhs.denom);
464        ApproxRational {
465            num: lnum - rnum,
466            denom,
467        }
468    }
469}
470impl Mul<&'_ ApproxRational> for &'_ ApproxRational {
471    type Output = ApproxRational;
472    fn mul(self, rhs: &'_ ApproxRational) -> ApproxRational {
473        ApproxRational {
474            num: &self.num * &rhs.num,
475            denom: &self.denom * &rhs.denom,
476        }
477    }
478}
479impl Mul<&'_ BigInt> for &'_ ApproxRational {
480    type Output = ApproxRational;
481    fn mul(self, rhs: &'_ BigInt) -> ApproxRational {
482        ApproxRational {
483            num: &self.num * rhs,
484            denom: self.denom.clone(),
485        }
486    }
487}
488impl Div<&'_ BigInt> for &'_ ApproxRational {
489    type Output = ApproxRational;
490    #[allow(clippy::suspicious_arithmetic_impl)]
491    fn div(self, rhs: &'_ BigInt) -> ApproxRational {
492        ApproxRational {
493            num: self.num.clone(),
494            denom: &self.denom * &Denom::decompose(rhs.clone()),
495        }
496    }
497}
498
499impl AddAssign for ApproxRational {
500    fn add_assign(&mut self, mut rhs: Self) {
501        self.denom = Denom::normalize(&mut self.num, &mut rhs.num, &self.denom, &rhs.denom);
502        self.num += rhs.num;
503    }
504}
505impl SubAssign for ApproxRational {
506    fn sub_assign(&mut self, mut rhs: Self) {
507        self.denom = Denom::normalize(&mut self.num, &mut rhs.num, &self.denom, &rhs.denom);
508        self.num -= rhs.num;
509    }
510}
511impl MulAssign for ApproxRational {
512    fn mul_assign(&mut self, rhs: Self) {
513        self.denom = &self.denom * &rhs.denom;
514        self.num *= rhs.num;
515    }
516}
517
518impl AddAssign<&'_ Self> for ApproxRational {
519    fn add_assign(&mut self, rhs: &'_ Self) {
520        let mut rnum = rhs.num.clone();
521        self.denom = Denom::normalize(&mut self.num, &mut rnum, &self.denom, &rhs.denom);
522        self.num += &rnum;
523    }
524}
525impl SubAssign<&'_ Self> for ApproxRational {
526    fn sub_assign(&mut self, rhs: &'_ Self) {
527        let mut rnum = rhs.num.clone();
528        self.denom = Denom::normalize(&mut self.num, &mut rnum, &self.denom, &rhs.denom);
529        self.num -= &rnum;
530    }
531}
532impl MulAssign<&'_ Self> for ApproxRational {
533    fn mul_assign(&mut self, rhs: &'_ Self) {
534        self.denom = &self.denom * &rhs.denom;
535        self.num *= &rhs.num;
536    }
537}
538impl DivAssign<&'_ BigInt> for ApproxRational {
539    #[allow(clippy::suspicious_op_assign_impl)]
540    fn div_assign(&mut self, rhs: &'_ BigInt) {
541        self.denom = &self.denom * &Denom::decompose(rhs.clone());
542    }
543}
544
545impl Sum for ApproxRational {
546    fn sum<I>(iter: I) -> Self
547    where
548        I: Iterator<Item = Self>,
549    {
550        iter.fold(Self::zero(), |acc, x| acc + x)
551    }
552}
553
554impl Product for ApproxRational {
555    fn product<I>(iter: I) -> Self
556    where
557        I: Iterator<Item = Self>,
558    {
559        iter.fold(Self::one(), |acc, x| acc * x)
560    }
561}
562
563impl<'a> Sum<&'a Self> for ApproxRational {
564    fn sum<I>(iter: I) -> Self
565    where
566        I: Iterator<Item = &'a Self>,
567    {
568        iter.fold(Self::zero(), |acc, x| acc + x)
569    }
570}
571
572impl RationalRef<&BigInt, ApproxRational> for &ApproxRational {}
573
574impl Rational<BigInt> for ApproxRational {
575    fn from_int(i: BigInt) -> Self {
576        ApproxRational {
577            num: i,
578            denom: Denom::one(),
579        }
580    }
581
582    fn ratio_i(num: BigInt, denom: BigInt) -> Self {
583        ApproxRational {
584            num,
585            denom: Denom::decompose(denom),
586        }
587    }
588
589    fn to_f64(&self) -> f64 {
590        Rational::to_f64(&self.reduced())
591    }
592
593    fn epsilon() -> Self {
594        Self::zero()
595    }
596
597    fn is_exact() -> bool {
598        true
599    }
600
601    fn description() -> &'static str {
602        "exact rational arithmetic with rounding of keep factors (6 decimal places)"
603    }
604
605    fn div_up_as_keep_factor(&self, rhs: &Self) -> Self {
606        let precision = BigInt::from(1_000_000);
607
608        let ldenom = self.denom.to_bigint();
609        let rdenom = rhs.denom.to_bigint();
610
611        let num = &self.num * precision * rdenom;
612        let denom = ldenom * &rhs.num;
613
614        Self {
615            num: (num + &denom - 1) / denom,
616            denom: Denom::precision(),
617        }
618    }
619
620    #[cfg(test)]
621    fn get_positive_test_values() -> Vec<Self> {
622        let mut result = Vec::new();
623        for i in 0..=30 {
624            result.push(Self::ratio(1 << i, 1));
625        }
626        for i in 0..=30 {
627            result.push(Self::ratio(1, 1 << i));
628        }
629        for i in 0..=30 {
630            result.push(Self::ratio(0x7FFF_FFFF - (1 << i), 1));
631        }
632        for i in 0..=30 {
633            result.push(Self::ratio(1, 0x7FFF_FFFF - (1 << i)));
634        }
635        result
636    }
637}
638
639#[cfg(test)]
640mod test {
641    use super::*;
642    use crate::{big_numeric_tests, numeric_benchmarks, numeric_tests};
643
644    fn make_ratio(num: i64, denom: i64) -> ApproxRational {
645        ApproxRational::ratio_i(BigInt::from(num), BigInt::from(denom))
646    }
647
648    numeric_tests!(
649        BigInt,
650        ApproxRational,
651        test_values_are_positive,
652        test_is_exact,
653        test_ratio,
654        test_ratio_invert,
655        test_is_zero,
656        test_zero_is_add_neutral,
657        test_add_is_commutative,
658        test_opposite,
659        test_sub_self,
660        test_add_sub,
661        test_sub_add,
662        test_one_is_mul_neutral,
663        test_mul_is_commutative,
664        test_mul_up_is_commutative,
665        test_mul_up_integers,
666        test_mul_up_wrt_mul,
667        test_one_is_div_up_neutral => fail(r"assertion `left == right` failed: div_up(a, 1) != a for 1/128
668  left: ApproxRational { num: 7813, denom: Denom { primes: [6, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], remainder: None } }
669 right: ApproxRational { num: 1, denom: Denom { primes: [7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], remainder: None } }"),
670        test_div_up_self,
671        test_mul_div_up => fail(r"assertion `left == right` failed: div_up(a * b, b) != a for 1/128, 1
672  left: ApproxRational { num: 7813, denom: Denom { primes: [6, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], remainder: None } }
673 right: ApproxRational { num: 1, denom: Denom { primes: [7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], remainder: None } }"),
674        test_mul_by_int,
675        test_mul_div_by_int,
676        test_references,
677        test_assign,
678    );
679
680    big_numeric_tests!(
681        BigInt,
682        ApproxRational,
683        Some(100_000),
684        test_add_is_associative,
685        test_mul_is_associative,
686        test_mul_is_distributive,
687        test_mul_by_int_is_associative,
688        test_mul_by_int_is_distributive,
689        test_div_by_int_is_associative,
690        test_div_by_int_is_distributive,
691        test_sum,
692        test_product,
693    );
694
695    numeric_benchmarks!(
696        BigInt,
697        ApproxRational,
698        bench_add,
699        bench_sub,
700        bench_mul,
701        bench_div_up,
702    );
703
704    #[test]
705    fn test_description() {
706        assert_eq!(
707            ApproxRational::description(),
708            "exact rational arithmetic with rounding of keep factors (6 decimal places)"
709        );
710    }
711
712    #[test]
713    fn test_display() {
714        assert_eq!(format!("{}", ApproxRational::zero()), "0");
715        assert_eq!(format!("{}", ApproxRational::one()), "1");
716        assert_eq!(format!("{}", make_ratio(0, 1)), "0");
717        assert_eq!(format!("{}", make_ratio(1, 1)), "1");
718        assert_eq!(format!("{}", make_ratio(1, 2)), "1/2");
719        assert_eq!(format!("{}", make_ratio(1, 23456789)), "1/23456789");
720        assert_eq!(format!("{}", make_ratio(123456789, 1)), "123456789");
721        assert_eq!(format!("{}", make_ratio(-1, 1)), "-1");
722        assert_eq!(format!("{}", make_ratio(-1, 2)), "-1/2");
723        assert_eq!(format!("{}", make_ratio(2, 2)), "1");
724        assert_eq!(format!("{}", make_ratio(60, 14)), "30/7");
725    }
726
727    #[test]
728    fn test_display_test_values() {
729        #[rustfmt::skip]
730        let expected_displays = [
731            "1", "2", "4", "8", "16", "32", "64", "128", "256", "512", "1024", "2048", "4096",
732            "8192", "16384", "32768", "65536", "131072", "262144", "524288", "1048576", "2097152",
733            "4194304", "8388608", "16777216", "33554432", "67108864", "134217728", "268435456",
734            "536870912", "1073741824",
735            "1", "1/2", "1/4", "1/8", "1/16", "1/32", "1/64", "1/128", "1/256", "1/512", "1/1024",
736            "1/2048", "1/4096", "1/8192", "1/16384", "1/32768", "1/65536", "1/131072", "1/262144",
737            "1/524288", "1/1048576", "1/2097152", "1/4194304", "1/8388608", "1/16777216",
738            "1/33554432", "1/67108864", "1/134217728", "1/268435456", "1/536870912", "1/1073741824",
739            "2147483646", "2147483645", "2147483643", "2147483639", "2147483631", "2147483615",
740            "2147483583", "2147483519", "2147483391", "2147483135", "2147482623", "2147481599",
741            "2147479551", "2147475455", "2147467263", "2147450879", "2147418111", "2147352575",
742            "2147221503", "2146959359", "2146435071", "2145386495", "2143289343", "2139095039",
743            "2130706431", "2113929215", "2080374783", "2013265919", "1879048191", "1610612735",
744            "1073741823",
745            "1/2147483646", "1/2147483645", "1/2147483643", "1/2147483639", "1/2147483631",
746            "1/2147483615", "1/2147483583", "1/2147483519", "1/2147483391", "1/2147483135",
747            "1/2147482623", "1/2147481599", "1/2147479551", "1/2147475455", "1/2147467263",
748            "1/2147450879", "1/2147418111", "1/2147352575", "1/2147221503", "1/2146959359",
749            "1/2146435071", "1/2145386495", "1/2143289343", "1/2139095039", "1/2130706431",
750            "1/2113929215", "1/2080374783", "1/2013265919", "1/1879048191", "1/1610612735",
751            "1/1073741823",
752        ];
753        let actual_displays: Vec<String> = ApproxRational::get_positive_test_values()
754            .iter()
755            .map(|x| format!("{x}"))
756            .collect();
757        assert_eq!(actual_displays, expected_displays);
758    }
759
760    mod denom {
761        use super::*;
762
763        #[test]
764        fn test_decompose_small() {
765            for (i, x) in Denom::DECOMPOSED.iter().enumerate() {
766                assert_eq!(x, &Denom::decompose_small(i as u64 + 1));
767                assert_eq!(x, &Denom::decompose(BigInt::from(i + 1)));
768                assert_eq!(x, &Denom::decompose_now(BigInt::from(i + 1)));
769            }
770        }
771
772        #[test]
773        #[should_panic(expected = "Failed to decompose small integer into small prime factors.")]
774        fn test_decompose_small_out_of_range() {
775            assert_eq!(Denom::decompose_small(97).to_bigint(), BigInt::from(97));
776        }
777
778        #[test]
779        fn test_decompose_is_correct() {
780            for i in 1..=1000 {
781                let bigi = BigInt::from(i);
782                let x = Denom::decompose(bigi.clone());
783                let mut recomposed = x.remainder.unwrap_or_else(BigInt::one);
784                for (i, &prime) in Denom::PRIMES.iter().enumerate() {
785                    for _ in 0..x.primes[i] {
786                        recomposed *= prime;
787                    }
788                }
789                assert_eq!(recomposed, bigi);
790            }
791        }
792
793        #[test]
794        fn test_decompose_known_values() {
795            assert_eq!(
796                Denom::decompose(BigInt::from(128)),
797                Denom {
798                    primes: [
799                        7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
800                    ],
801                    remainder: None,
802                }
803            );
804            assert_eq!(
805                Denom::decompose(BigInt::from(89)),
806                Denom {
807                    primes: [
808                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1
809                    ],
810                    remainder: None,
811                }
812            );
813            assert_eq!(
814                Denom::decompose(BigInt::from(97)),
815                Denom {
816                    primes: [
817                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
818                    ],
819                    remainder: Some(BigInt::from(97)),
820                }
821            );
822            assert_eq!(
823                Denom::decompose(BigInt::from(97000)),
824                Denom {
825                    primes: [
826                        3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
827                    ],
828                    remainder: Some(BigInt::from(97)),
829                }
830            );
831        }
832
833        #[test]
834        fn test_to_bigint() {
835            for i in 1..=1000 {
836                let bigi = BigInt::from(i);
837                let x = Denom::decompose(bigi.clone());
838                assert_eq!(x.to_bigint(), bigi);
839            }
840        }
841
842        #[test]
843        fn test_product() {
844            let values = (100..200)
845                .map(|i| Denom::decompose(BigInt::from(i)))
846                .collect::<Vec<_>>();
847            for (i, x) in values.iter().enumerate().map(|(i, x)| (i + 100, x)) {
848                for (j, y) in values.iter().enumerate().map(|(j, y)| (j + 100, y)) {
849                    let z = x * y;
850                    assert_eq!(z, Denom::decompose(BigInt::from(i * j)));
851                    for k in 0..Denom::NUM_PRIMES {
852                        assert_eq!(z.primes[k], x.primes[k] + y.primes[k]);
853                    }
854                }
855            }
856        }
857
858        #[test]
859        fn test_normalize() {
860            let values = (100..200)
861                .map(|i| Denom::decompose(BigInt::from(i)))
862                .collect::<Vec<_>>();
863            for x in &values {
864                for y in &values {
865                    let mut xnum = BigInt::one();
866                    let mut ynum = BigInt::one();
867                    let lcm = Denom::normalize(&mut xnum, &mut ynum, x, y);
868                    let lcm_bigint = lcm.to_bigint();
869
870                    assert_eq!(xnum * x.to_bigint(), lcm_bigint);
871                    assert_eq!(ynum * y.to_bigint(), lcm_bigint);
872                    for k in 0..Denom::NUM_PRIMES {
873                        assert_eq!(lcm.primes[k], std::cmp::max(x.primes[k], y.primes[k]));
874                    }
875                }
876            }
877        }
878    }
879}
stv_rs/arithmetic/approx_rational.rs

stv_rs/arithmetic/
approx_rational.rs
