feanor_math/algorithms/eea/
mod.rs

1use std::cmp::Ordering;
2use std::mem::swap;
3
4use crate::integer::*;
5use crate::ordered::{OrderedRing, OrderedRingStore};
6use crate::pid::*;
7use crate::ring::*;
8
9/// For `a, b` computes `s, t, d` such that `s*a + t*b == d` is a greatest
10/// common divisor of `a` and `b`.
11///
12/// The gcd `d` is only unique up to units, and `s, t` are not unique at all.
13/// No guarantees are given on which of these solutions is returned. For integers,
14/// see [`signed_eea()`] which gives more guarantees.
15///
16/// Note that this function always uses the euclidean algorithm to compute these values.
17/// In most cases, it is instead recommended to use [`PrincipalIdealRing::extended_ideal_gen()`],
18/// which uses a ring-specific algorithm to compute the Bezout identity (which will of
19/// course be [`eea()`] in some cases).
20pub fn eea<R>(a: El<R>, b: El<R>, ring: R) -> (El<R>, El<R>, El<R>)
21where
22    R: RingStore,
23    R::Type: EuclideanRing,
24{
25    let (mut a, mut b) = (a, b);
26
27    let (mut sa, mut ta) = (ring.one(), ring.zero());
28    let (mut sb, mut tb) = (ring.zero(), ring.one());
29
30    while !ring.is_zero(&b) {
31        // loop invariants; unfortunately, this evaluation might cause an integer overflow
32        // debug_assert!(ring.eq_el(&a, &ring.add(ring.mul_ref(&sa, &fst), ring.mul_ref(&ta,
33        // &snd)))); debug_assert!(ring.eq_el(&b, &ring.add(ring.mul_ref(&sb, &fst),
34        // ring.mul_ref(&tb, &snd))));
35
36        let (quo, rem) = ring.euclidean_div_rem(a, &b);
37        ta = ring.sub(ta, ring.mul_ref(&quo, &tb));
38        sa = ring.sub(sa, ring.mul_ref_snd(quo, &sb));
39        a = rem;
40
41        swap(&mut a, &mut b);
42        swap(&mut sa, &mut sb);
43        swap(&mut ta, &mut tb);
44    }
45    return (sa, ta, a);
46}
47
48/// Computes the gcd `d` of `a` and `b`, together with "half a Bezout identity", i.e.
49/// some `s` such that `s * a = d mod b`.
50///
51/// For details, see [`eea()`].
52#[stability::unstable(feature = "enable")]
53pub fn half_eea<R>(a: El<R>, b: El<R>, ring: R) -> (El<R>, El<R>)
54where
55    R: RingStore,
56    R::Type: EuclideanRing,
57{
58    let (mut a, mut b) = (a, b);
59    let (mut s, mut t) = (ring.one(), ring.zero());
60
61    // invariant: `s * a == a mod b` and `t * a == b mod b`
62    while !ring.is_zero(&b) {
63        let (q, r) = ring.euclidean_div_rem(a, &b);
64        a = r;
65        ring.sub_assign(&mut s, ring.mul_ref_snd(q, &t));
66        swap(&mut a, &mut b);
67        swap(&mut s, &mut t);
68    }
69    return (s, a);
70}
71
72/// For integers a, b finds the smallest integers s, t so that
73/// `s*a + t*b == gcd(a, b)` is the greatest common divisor of a, b.
74///
75/// Details: s and t are not unique, this function will return
76/// the smallest tuple (s, t) (ordered by the total ordering
77/// given by `(s, t) ≤ (u, v) :<=> |s| ≤ |u| and |t| ≤ |v|`).
78/// In the case |a| = |b|, there are two minimal elements, in this case, it is
79/// unspecified whether this function returns (±1, 0, a) or (0, ±1, a).
80/// We define the greatest common divisor gcd(a, b) as the minimal
81/// element of the set of integers dividing a and b (ordered by divisibility),
82/// whose sign matches the sign of a.
83///
84/// In particular, have
85/// ```rust
86/// # use feanor_math::algorithms::eea::signed_gcd;
87/// # use feanor_math::primitive_int::*;
88/// assert_eq!(2, signed_gcd(6, 8, &StaticRing::<i64>::RING));
89/// assert_eq!(0, signed_gcd(0, 0, &StaticRing::<i64>::RING));
90/// assert_eq!(5, signed_gcd(0, -5, &StaticRing::<i64>::RING));
91/// assert_eq!(-5, signed_gcd(-5, 0, &StaticRing::<i64>::RING));
92/// assert_eq!(-1, signed_gcd(-1, 1, &StaticRing::<i64>::RING));
93/// assert_eq!(1, signed_gcd(1, -1, &StaticRing::<i64>::RING));
94/// ```
95/// and therefore `signed_eea(6, 8) == (-1, 1, 2)`,
96/// `signed_eea(-6, 8) == (-1, -1, -2)`,
97/// `signed_eea(8, -6) == (1, 1, 2)`,
98/// `signed_eea(0, 0) == (0, 0, 0)`
99pub fn signed_eea<R>(fst: El<R>, snd: El<R>, ring: R) -> (El<R>, El<R>, El<R>)
100where
101    R: RingStore,
102    R::Type: EuclideanRing + OrderedRing,
103{
104    if ring.is_zero(&fst) {
105        return match ring.cmp(&snd, &ring.zero()) {
106            Ordering::Equal => (ring.zero(), ring.zero(), ring.zero()),
107            Ordering::Less => (ring.zero(), ring.negate(ring.one()), ring.negate(snd)),
108            Ordering::Greater => (ring.zero(), ring.one(), snd),
109        };
110    }
111    let fst_negative = ring.cmp(&fst, &ring.zero());
112
113    let (s, t, d) = eea(fst, snd, &ring);
114
115    // the sign is not consistent (potentially toggled each iteration),
116    // so normalize here
117    if ring.cmp(&d, &ring.zero()) == fst_negative {
118        return (s, t, d);
119    } else {
120        return (ring.negate(s), ring.negate(t), ring.negate(d));
121    }
122}
123
124/// Finds a greatest common divisor of a and b.
125///  
126/// The gcd of two elements `a, b` in a euclidean ring is the (w.r.t divisibility) greatest
127/// element that divides both elements, i.e. the greatest element (w.r.t. divisibility) `g` such
128/// that `g | a, b`.
129///
130/// Note that this function always uses the euclidean algorithm to compute the gcd. In most
131/// cases, it is instead recommended to use [`PrincipalIdealRing::ideal_gen()`], which uses
132/// a ring-specific algorithm to compute the gcd (which will of course be [`gcd()`] in some cases).
133///
134/// In general, the gcd is only unique up to multiplication by units. For integers, the function
135/// [`signed_gcd()`] gives more guarantees.
136pub fn gcd<R>(a: El<R>, b: El<R>, ring: R) -> El<R>
137where
138    R: RingStore,
139    R::Type: EuclideanRing,
140{
141    let (mut a, mut b) = (a, b);
142
143    // invariant: `gcd(a, b) = gcd(original_a, original_b)`
144    while !ring.is_zero(&b) {
145        let (_, r) = ring.euclidean_div_rem(a, &b);
146        a = b;
147        b = r;
148    }
149    return a;
150}
151
152/// Finds the greatest common divisor of a and b.
153///
154/// The gcd is only unique up to multiplication by units, so in this case up to sign.
155/// However, this function guarantees the following behavior w.r.t different signs:
156///
157/// ```text
158///   a < 0 => gcd(a, b) < 0
159///   a > 0 => gcd(a, b) > 0
160///   sign of b is irrelevant
161///   gcd(0, 0) = 0
162/// ```
163pub fn signed_gcd<R>(a: El<R>, b: El<R>, ring: R) -> El<R>
164where
165    R: RingStore,
166    R::Type: EuclideanRing + OrderedRing,
167{
168    let (_, _, d) = signed_eea(a, b, ring);
169    return d;
170}
171
172/// Finds the least common multiple of two elements in an ordered euclidean ring,
173/// e.g. of two integers.
174///
175/// The general lcm is only unique up to multiplication by units. For `signed_lcm`,
176/// the following behavior is guaranteed:
177/// ```text
178///   b > 0 => lcm(a, b) >= 0
179///   b < 0 => lcm(a, b) <= 0
180///   lcm(0, b) = lcm(a, 0) = lcm(0, 0) = 0
181/// ```
182pub fn signed_lcm<R>(fst: El<R>, snd: El<R>, ring: R) -> El<R>
183where
184    R: RingStore,
185    R::Type: EuclideanRing + OrderedRing,
186{
187    if ring.is_zero(&fst) || ring.is_zero(&snd) {
188        ring.zero()
189    } else {
190        ring.mul(
191            ring.euclidean_div(ring.clone_el(&fst), &signed_gcd(fst, ring.clone_el(&snd), &ring)),
192            snd,
193        )
194    }
195}
196
197/// Finds the least common multiple of two elements `a, b` in a euclidean ring, i.e. the smallest
198/// (w.r.t. divisibility) element `y` with `a, b | y`.
199///
200/// In general, the lcm is only unique up to multiplication by units. For integers, the function
201/// [`signed_lcm()`] gives more guarantees.
202pub fn lcm<R>(fst: El<R>, snd: El<R>, ring: R) -> El<R>
203where
204    R: RingStore,
205    R::Type: EuclideanRing,
206{
207    if ring.is_zero(&fst) || ring.is_zero(&snd) {
208        ring.zero()
209    } else {
210        ring.euclidean_div(ring.mul_ref(&fst, &snd), &gcd(fst, snd, &ring))
211    }
212}
213
214/// Computes x such that `x = a mod p` and `x = b mod q`. Requires that p and q are coprime.
215pub fn inv_crt<I>(a: El<I>, b: El<I>, p: &El<I>, q: &El<I>, ZZ: I) -> El<I>
216where
217    I: RingStore,
218    I::Type: IntegerRing,
219{
220    let (s, t, d) = signed_eea(ZZ.clone_el(p), ZZ.clone_el(q), &ZZ);
221    assert!(ZZ.is_one(&d) || ZZ.is_neg_one(&d));
222    let mut result = ZZ.add(ZZ.prod([a, t, ZZ.clone_el(q)]), ZZ.prod([b, s, ZZ.clone_el(p)]));
223
224    let n = ZZ.mul_ref(p, q);
225    result = ZZ.euclidean_rem(result, &n);
226    if ZZ.is_neg(&result) {
227        ZZ.add_assign(&mut result, n);
228    }
229    return result;
230}
231
232#[cfg(test)]
233use crate::primitive_int::*;
234
235#[test]
236fn test_gcd() {
237    assert_eq!(3, gcd(15, 6, &StaticRing::<i64>::RING).abs());
238    assert_eq!(3, gcd(6, 15, &StaticRing::<i64>::RING).abs());
239
240    assert_eq!(7, gcd(0, 7, &StaticRing::<i64>::RING).abs());
241    assert_eq!(7, gcd(7, 0, &StaticRing::<i64>::RING).abs());
242    assert_eq!(0, gcd(0, 0, &StaticRing::<i64>::RING).abs());
243
244    assert_eq!(1, gcd(9, 1, &StaticRing::<i64>::RING).abs());
245    assert_eq!(1, gcd(1, 9, &StaticRing::<i64>::RING).abs());
246
247    assert_eq!(1, gcd(13, 300, &StaticRing::<i64>::RING).abs());
248    assert_eq!(1, gcd(300, 13, &StaticRing::<i64>::RING).abs());
249
250    assert_eq!(3, gcd(-15, 6, &StaticRing::<i64>::RING).abs());
251    assert_eq!(3, gcd(6, -15, &StaticRing::<i64>::RING).abs());
252    assert_eq!(3, gcd(-6, -15, &StaticRing::<i64>::RING).abs());
253}
254
255#[test]
256fn test_signed_gcd() {
257    assert_eq!(3, signed_gcd(15, 6, &StaticRing::<i64>::RING));
258    assert_eq!(3, signed_gcd(6, 15, &StaticRing::<i64>::RING));
259
260    assert_eq!(7, signed_gcd(0, 7, &StaticRing::<i64>::RING));
261    assert_eq!(7, signed_gcd(7, 0, &StaticRing::<i64>::RING));
262    assert_eq!(0, signed_gcd(0, 0, &StaticRing::<i64>::RING));
263
264    assert_eq!(1, signed_gcd(9, 1, &StaticRing::<i64>::RING));
265    assert_eq!(1, signed_gcd(1, 9, &StaticRing::<i64>::RING));
266
267    assert_eq!(1, signed_gcd(13, 300, &StaticRing::<i64>::RING));
268    assert_eq!(1, signed_gcd(300, 13, &StaticRing::<i64>::RING));
269
270    assert_eq!(-3, signed_gcd(-15, 6, &StaticRing::<i64>::RING));
271    assert_eq!(3, signed_gcd(6, -15, &StaticRing::<i64>::RING));
272    assert_eq!(-3, signed_gcd(-6, -15, &StaticRing::<i64>::RING));
273}
274
275#[test]
276fn test_eea_sign() {
277    assert_eq!((2, -1, 1), signed_eea(3, 5, &StaticRing::<i64>::RING));
278    assert_eq!((-1, 2, 1), signed_eea(5, 3, &StaticRing::<i64>::RING));
279    assert_eq!((2, 1, -1), signed_eea(-3, 5, &StaticRing::<i64>::RING));
280    assert_eq!((-1, -2, 1), signed_eea(5, -3, &StaticRing::<i64>::RING));
281    assert_eq!((2, 1, 1), signed_eea(3, -5, &StaticRing::<i64>::RING));
282    assert_eq!((-1, -2, -1), signed_eea(-5, 3, &StaticRing::<i64>::RING));
283    assert_eq!((2, -1, -1), signed_eea(-3, -5, &StaticRing::<i64>::RING));
284    assert_eq!((-1, 2, -1), signed_eea(-5, -3, &StaticRing::<i64>::RING));
285    assert_eq!((0, 0, 0), signed_eea(0, 0, &StaticRing::<i64>::RING));
286    assert_eq!((1, 0, 4), signed_eea(4, 0, &StaticRing::<i64>::RING));
287    assert_eq!((0, 1, 4), signed_eea(0, 4, &StaticRing::<i64>::RING));
288    assert_eq!((1, 0, -4), signed_eea(-4, 0, &StaticRing::<i64>::RING));
289    assert_eq!((0, -1, 4), signed_eea(0, -4, &StaticRing::<i64>::RING));
290}
291
292#[test]
293fn test_signed_eea() {
294    assert_eq!((-1, 1, 2), signed_eea(6, 8, &StaticRing::<i64>::RING));
295    assert_eq!((2, -1, 5), signed_eea(15, 25, &StaticRing::<i64>::RING));
296    assert_eq!((4, -7, 2), signed_eea(32, 18, &StaticRing::<i64>::RING));
297}
298
299#[test]
300fn test_signed_lcm() {
301    assert_eq!(24, signed_lcm(6, 8, &StaticRing::<i64>::RING));
302    assert_eq!(24, signed_lcm(-6, 8, &StaticRing::<i64>::RING));
303    assert_eq!(-24, signed_lcm(6, -8, &StaticRing::<i64>::RING));
304    assert_eq!(-24, signed_lcm(-6, -8, &StaticRing::<i64>::RING));
305    assert_eq!(0, signed_lcm(0, 0, &StaticRing::<i64>::RING));
306    assert_eq!(0, signed_lcm(6, 0, &StaticRing::<i64>::RING));
307    assert_eq!(0, signed_lcm(0, 8, &StaticRing::<i64>::RING));
308}
feanor_math/algorithms/eea/mod.rs

feanor_math/algorithms/eea/
mod.rs
