feanor_math/algorithms/eea/
mod.rs

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

feanor_math/algorithms/eea/
mod.rs
