malachite_base/num/arithmetic/
gcd.rs

1// Copyright © 2025 Mikhail Hogrefe
2//
3// Uses code adopted from the FLINT Library.
4//
5//      Copyright © 2009, 2016 William Hart
6//
7// This file is part of Malachite.
8//
9// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
10// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
11// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.
12
13use crate::num::arithmetic::traits::{Gcd, GcdAssign};
14use crate::num::basic::unsigneds::PrimitiveUnsigned;
15use core::cmp::min;
16
17#[cfg(feature = "test_build")]
18pub fn gcd_euclidean<T: PrimitiveUnsigned>(x: T, y: T) -> T {
19    if y == T::ZERO {
20        x
21    } else {
22        gcd_euclidean(y, x % y)
23    }
24}
25
26#[cfg(feature = "test_build")]
27pub fn gcd_binary<T: PrimitiveUnsigned>(x: T, y: T) -> T {
28    if x == y {
29        x
30    } else if x == T::ZERO {
31        y
32    } else if y == T::ZERO {
33        x
34    } else if x.even() {
35        if y.odd() {
36            gcd_binary(x >> 1, y)
37        } else {
38            gcd_binary(x >> 1, y >> 1) << 1
39        }
40    } else if y.even() {
41        gcd_binary(x, y >> 1)
42    } else if x > y {
43        gcd_binary((x - y) >> 1, y)
44    } else {
45        gcd_binary((y - x) >> 1, x)
46    }
47}
48
49pub_test! {gcd_fast_a<T: PrimitiveUnsigned>(mut x: T, mut y: T) -> T {
50    if x == T::ZERO {
51        return y;
52    }
53    if y == T::ZERO {
54        return x;
55    }
56    let x_zeros = x.trailing_zeros();
57    let y_zeros = y.trailing_zeros();
58    let f = min(x_zeros, y_zeros);
59    x >>= x_zeros;
60    y >>= y_zeros;
61    while x != y {
62        if x < y {
63            y -= x;
64            y >>= y.trailing_zeros();
65        } else {
66            x -= y;
67            x >>= x.trailing_zeros();
68        }
69    }
70    x << f
71}}
72
73#[cfg(feature = "test_build")]
74// This is a modified version of `n_xgcd` from `ulong_extras/xgcd.c`, FLINT 2.7.1.
75pub fn gcd_fast_b<T: PrimitiveUnsigned>(mut x: T, y: T) -> T {
76    let mut v;
77    if x >= y {
78        v = y;
79    } else {
80        v = x;
81        x = y;
82    }
83    let mut d;
84    // x and y both have their top bit set.
85    if (x & v).get_highest_bit() {
86        d = x - v;
87        x = v;
88        v = d;
89    }
90    // The second value has its second-highest set.
91    while (v << 1u32).get_highest_bit() {
92        d = x - v;
93        x = v;
94        if d < v {
95            v = d;
96        } else if d < (v << 1) {
97            v = d - x;
98        } else {
99            v = d - (x << 1);
100        }
101    }
102    while v != T::ZERO {
103        // Overflow is not possible due to top 2 bits of v not being set. Avoid divisions when
104        // quotient < 4.
105        if x < (v << 2) {
106            d = x - v;
107            x = v;
108            if d < v {
109                v = d;
110            } else if d < (v << 1) {
111                v = d - x;
112            } else {
113                v = d - (x << 1);
114            }
115        } else {
116            let rem = x % v;
117            x = v;
118            v = rem;
119        }
120    }
121    x
122}
123
124macro_rules! impl_gcd {
125    ($t:ident) => {
126        impl Gcd<$t> for $t {
127            type Output = $t;
128
129            /// Computes the GCD (greatest common divisor) of two numbers.
130            ///
131            /// The GCD of 0 and $n$, for any $n$, is 0. In particular, $\gcd(0, 0) = 0$, which
132            /// makes sense if we interpret "greatest" to mean "greatest by the divisibility order".
133            ///
134            /// $$
135            /// f(x, y) = \gcd(x, y).
136            /// $$
137            ///
138            /// # Worst-case complexity
139            /// $T(n) = O(n^2)$
140            ///
141            /// $M(n) = O(n)$
142            ///
143            /// where $T$ is time, $M$ is additional memory, and $n$ is
144            /// `max(self.significant_bits(), other.significant_bits())`.
145            ///
146            /// # Examples
147            /// See [here](super::gcd#gcd).
148            #[inline]
149            fn gcd(self, other: $t) -> $t {
150                gcd_fast_a(self, other)
151            }
152        }
153
154        impl GcdAssign<$t> for $t {
155            /// Replaces another with the GCD (greatest common divisor) of it and another number.
156            ///
157            /// The GCD of 0 and $n$, for any $n$, is 0. In particular, $\gcd(0, 0) = 0$, which
158            /// makes sense if we interpret "greatest" to mean "greatest by the divisibility order".
159            ///
160            /// $$
161            /// x \gets \gcd(x, y).
162            /// $$
163            ///
164            /// # Worst-case complexity
165            /// $T(n) = O(n^2)$
166            ///
167            /// $M(n) = O(n)$
168            ///
169            /// where $T$ is time, $M$ is additional memory, and $n$ is
170            /// `max(self.significant_bits(), other.significant_bits())`.
171            ///
172            /// # Examples
173            /// See [here](super::gcd#gcd_assign).
174            #[inline]
175            fn gcd_assign(&mut self, other: $t) {
176                *self = gcd_fast_a(*self, other);
177            }
178        }
179    };
180}
181apply_to_unsigneds!(impl_gcd);
malachite_base/num/arithmetic/gcd.rs

malachite_base/num/arithmetic/
gcd.rs
