relp_num/rational/big/
gcd.rs

1//! # The greatest common divisor of two arbitrary precision rationals
2use smallvec::SmallVec;
3
4use crate::integer::big::ops::normalize::{gcd, simplify_fraction_without_info};
5use crate::integer::big::{NonZeroUbig, Ubig};
6use crate::rational::big::Big;
7use crate::sign::Sign;
8use crate::traits::Gcd;
9
10/// Never fails: the numerator of the result divides both numerators and the denominator is a
11/// product of the two denominators, and neither can run out of room.
12impl<const S: usize> Gcd for Big<S> {
13    fn gcd(&self, other: &Self) -> Option<Self> {
14        Some(match (self.sign, other.sign) {
15            (Sign::Zero, Sign::Zero) => <Self as num_traits::Zero>::zero(),
16            // Every integer multiple of a value is a multiple of it, and zero is one of them, so
17            // the other value divides the pair. Only its magnitude does the dividing.
18            (Sign::Zero, _) => Self {
19                sign: Sign::Positive,
20                numerator: other.numerator.clone(),
21                denominator: other.denominator.clone(),
22            },
23            (_, Sign::Zero) => Self {
24                sign: Sign::Positive,
25                numerator: self.numerator.clone(),
26                denominator: self.denominator.clone(),
27            },
28            _ => Self {
29                sign: Sign::Positive,
30                numerator: numerator_gcd(&self.numerator, &other.numerator),
31                denominator: denominator_lcm(&self.denominator, &other.denominator),
32            },
33        })
34    }
35}
36
37/// The greatest common divisor of two values that are not zero.
38///
39/// The routine this calls refuses the cases that cost it an allocation to no purpose, so they are
40/// answered here.
41fn numerator_gcd<const S: usize>(left: &Ubig<S>, right: &Ubig<S>) -> Ubig<S> {
42    if num_traits::One::is_one(left) || num_traits::One::is_one(right) {
43        return Ubig::new(1);
44    }
45    if **left == **right {
46        return left.clone();
47    }
48
49    // SAFETY: Both are well formed and not empty, neither is one, and they differ.
50    unsafe { Ubig::from_inner_unchecked(gcd(left, right)) }
51}
52
53/// The least common multiple of two values that are not zero.
54///
55/// Written as `(left / gcd(left, right)) * right`, which is the order that keeps the intermediate
56/// no larger than the answer.
57fn denominator_lcm<const S: usize>(left: &NonZeroUbig<S>, right: &NonZeroUbig<S>) -> NonZeroUbig<S> {
58    let mut reduced = SmallVec::<[usize; S]>::from_slice(left);
59    let mut other = SmallVec::<[usize; S]>::from_slice(right);
60    // SAFETY: Both are well formed and not zero.
61    unsafe { simplify_fraction_without_info(&mut reduced, &mut other) };
62
63    // SAFETY: Dividing a value that is not zero by something that divides it leaves it not zero,
64    // and well formed.
65    let reduced = unsafe { NonZeroUbig::from_inner_unchecked(reduced) };
66    reduced * right.clone()
67}
68
69#[cfg(test)]
70mod test {
71    use std::str::FromStr;
72
73    use crate::{Gcd, RationalBig, RB};
74
75    fn big(value: &str) -> RationalBig {
76        RationalBig::from_str(value).unwrap()
77    }
78
79    #[test]
80    fn integers_get_the_ordinary_answer() {
81        assert_eq!(RB!(4).gcd(&RB!(6)), Some(RB!(2)));
82        assert_eq!(RB!(7).gcd(&RB!(13)), Some(RB!(1)));
83        assert_eq!(RB!(12).gcd(&RB!(12)), Some(RB!(12)));
84        assert_eq!(RB!(1).gcd(&RB!(9)), Some(RB!(1)));
85        assert_eq!(RB!(9).gcd(&RB!(1)), Some(RB!(1)));
86    }
87
88    #[test]
89    fn the_sign_of_the_operands_says_nothing_about_what_divides_them() {
90        assert_eq!(RB!(-4).gcd(&RB!(6)), Some(RB!(2)));
91        assert_eq!(RB!(-4).gcd(&RB!(-6)), Some(RB!(2)));
92    }
93
94    #[test]
95    fn zero_is_a_multiple_of_everything() {
96        assert_eq!(RB!(0).gcd(&RB!(0)), Some(RB!(0)));
97        assert_eq!(RB!(0).gcd(&RB!(-6)), Some(RB!(6)));
98        assert_eq!(RB!(6, 5).gcd(&RB!(0)), Some(RB!(6, 5)));
99    }
100
101    /// `gcd(p1 / q1, p2 / q2) = gcd(p1, p2) / lcm(q1, q2)`.
102    #[test]
103    fn fractions_take_the_numerator_gcd_over_the_denominator_lcm() {
104        assert_eq!(RB!(1, 2).gcd(&RB!(1, 3)), Some(RB!(1, 6)));
105        assert_eq!(RB!(2, 3).gcd(&RB!(4, 9)), Some(RB!(2, 9)));
106        assert_eq!(RB!(3, 4).gcd(&RB!(5, 6)), Some(RB!(1, 12)));
107        assert_eq!(RB!(1, 2).gcd(&RB!(1, 2)), Some(RB!(1, 2)));
108    }
109
110    /// Dividing by the answer has to leave integers, and integers with nothing left in common.
111    #[test]
112    fn dividing_by_it_leaves_coprime_integers() {
113        for (left, right) in [
114            (RB!(1, 2), RB!(1, 3)),
115            (RB!(3, 4), RB!(5, 6)),
116            (RB!(-7, 15), RB!(14, 25)),
117            (RB!(9), RB!(6)),
118        ] {
119            let divisor = left.gcd(&right).unwrap();
120            let one = &left / &divisor;
121            let other = &right / &divisor;
122
123            assert_eq!(&one * &divisor, left);
124            assert_eq!(&other * &divisor, right);
125            assert_eq!(one.gcd(&other), Some(RB!(1)));
126        }
127    }
128
129    /// Against the definition, over every small fraction pair.
130    ///
131    /// The answer has to divide both, and nothing larger may: a candidate twice as large fails to
132    /// divide one of them, or the answer was not the greatest.
133    #[test]
134    fn every_small_pair_against_the_definition() {
135        // A value in lowest terms is an integer exactly when its denominator is one.
136        let divides = |value: &RationalBig, divisor: &RationalBig| {
137            num_traits::One::is_one(&(value / divisor).denominator)
138        };
139
140        for p1 in 1..10_i64 {
141            for q1 in 1..10_i64 {
142                for p2 in 1..10_i64 {
143                    for q2 in 1..10_i64 {
144                        let left = RationalBig::from((p1, q1));
145                        let right = RationalBig::from((p2, q2));
146                        let divisor = left.gcd(&right).unwrap();
147
148                        assert!(divides(&left, &divisor), "{left} / {divisor}");
149                        assert!(divides(&right, &divisor), "{right} / {divisor}");
150
151                        // Anything strictly larger has to fail on one of the two, or the answer
152                        // was not the greatest common one.
153                        for multiple in [2_i64, 3, 5] {
154                            let larger = &divisor * &RationalBig::from(multiple);
155                            assert!(
156                                !divides(&left, &larger) || !divides(&right, &larger),
157                                "{larger} divides both {left} and {right}",
158                            );
159                        }
160                    }
161                }
162            }
163        }
164    }
165
166    #[test]
167    fn values_that_do_not_fit_in_a_word() {
168        let left = big("340282366920938463463374607431768211456"); // 2 ** 128
169        let right = big("170141183460469231731687303715884105728"); // 2 ** 127
170        assert_eq!(left.gcd(&right), Some(right.clone()));
171
172        let third = &RationalBig::from(1) / &left;
173        assert_eq!(third.gcd(&third), Some(third));
174    }
175}
relp_num/rational/big/gcd.rs

relp_num/rational/big/
gcd.rs
