feanor_math/algorithms/
resultant.rs

1
2use std::mem::swap;
3
4use crate::delegate::{UnwrapHom, WrapHom};
5use crate::reduce_lift::poly_eval::{EvalPolyLocallyRing, EvaluatePolyLocallyReductionMap};
6use crate::divisibility::{DivisibilityRingStore, Domain};
7use crate::pid::{EuclideanRing, PrincipalIdealRing, PrincipalIdealRingStore};
8use crate::algorithms::eea::signed_lcm;
9use crate::rings::field::{AsField, AsFieldBase};
10use crate::rings::fraction::FractionFieldStore;
11use crate::rings::poly::*;
12use crate::rings::finite::*;
13use crate::pid::EuclideanRingStore;
14use crate::specialization::FiniteRingOperation;
15use crate::algorithms;
16use crate::integer::*;
17use crate::rings::rational::*;
18use crate::homomorphism::*;
19use crate::ring::*;
20use crate::rings::poly::dense_poly::DensePolyRing;
21
22///
23/// Computes the resultant of `f` and `g` over the base ring, using
24/// a direct variant of Eulid's algorithm. 
25/// 
26/// This function is deprecated, you should instead use
27/// [`ComputeResultantRing::resultant()`] or [`resultant_finite_field()`].
28/// 
29/// The resultant is the determinant of the linear map
30/// ```text
31///   R[X]_deg(g)  x  R[X]_deg(f)  ->  R[X]_deg(fg),
32///        a       ,       b       ->    af + bg
33/// ```
34/// where `R[X]_d` refers to the vector space of polynomials in `R[X]` of degree
35/// less than `d`.
36/// 
37/// # Example
38/// 
39/// ```rust
40/// use feanor_math::ring::*;
41/// use feanor_math::primitive_int::*;
42/// use feanor_math::rings::poly::dense_poly::DensePolyRing;
43/// use feanor_math::rings::poly::*;
44/// use feanor_math::algorithms::resultant::*;
45/// let ZZ = StaticRing::<i64>::RING;
46/// let ZZX = DensePolyRing::new(ZZ, "X");
47/// let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1]);
48/// // the discrimiant is the resultant of f and f'
49/// let discriminant = resultant_global(&ZZX, ZZX.clone_el(&f), derive_poly(&ZZX, &f));
50/// assert_eq!(0, discriminant);
51/// ```
52/// 
53#[deprecated]
54pub fn resultant_global<P>(ring: P, mut f: El<P>, mut g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
55    where P: RingStore + Copy,
56        P::Type: PolyRing,
57        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Domain + PrincipalIdealRing
58{
59    let base_ring = ring.base_ring();
60    if ring.is_zero(&g) || ring.is_zero(&f) {
61        return base_ring.zero();
62    }
63    let mut scale_den = base_ring.one();
64    let mut scale_num = base_ring.one();
65
66    if ring.degree(&g).unwrap() < ring.degree(&f).unwrap() {
67        if ring.degree(&g).unwrap() % 2 != 0 && ring.degree(&f).unwrap() % 2 != 0 {
68            base_ring.negate_inplace(&mut scale_num);
69        }
70        std::mem::swap(&mut f, &mut g);
71    }
72
73    while ring.degree(&f).unwrap_or(0) >= 1 {
74
75        let balance_factor = ring.get_ring().balance_poly(&mut f);
76        if let Some(balance_factor) = balance_factor {
77            base_ring.mul_assign(&mut scale_num, base_ring.pow(balance_factor, ring.degree(&g).unwrap()));
78        }
79
80        // use here that `res(f, g) = a^(-deg(f)) lc(f)^(deg(g) - deg(ag - fh)) res(f, ag - fh)` if `deg(fh) <= deg(g)`
81        let deg_g = ring.degree(&g).unwrap();
82        let (_q, r, a) = algorithms::poly_div::poly_div_rem_domain(ring, g, &f);
83        let deg_r = ring.degree(&r).unwrap_or(0);
84
85        // adjust the scaling factor - we cancel out gcd's to prevent excessive number growth
86        base_ring.mul_assign(&mut scale_den, base_ring.pow(a, ring.degree(&f).unwrap()));
87        base_ring.mul_assign(&mut scale_num, base_ring.pow(base_ring.clone_el(ring.lc(&f).unwrap()), deg_g - deg_r));
88        let gcd = base_ring.ideal_gen(&scale_den, &scale_num);
89        scale_den = base_ring.checked_div(&scale_den, &gcd).unwrap();
90        scale_num = base_ring.checked_div(&scale_num, &gcd).unwrap();
91
92        g = f;
93        f = r;
94
95        if ring.degree(&g).unwrap() % 2 != 0 && ring.degree(&f).unwrap_or(0) % 2 != 0 {
96            base_ring.negate_inplace(&mut scale_num);
97        }
98    }
99
100    if ring.is_zero(&f) {
101        return base_ring.zero();
102    } else {
103        let mut result = base_ring.clone_el(&ring.coefficient_at(&f, 0));
104        result = base_ring.pow(result, ring.degree(&g).unwrap());
105        base_ring.mul_assign(&mut result, scale_num);
106        return base_ring.checked_div(&result, &scale_den).unwrap();
107    }
108}
109
110///
111/// Computes the resultant of two polynomials `f` and `g` over a finite field.
112/// 
113/// Usually you should use [`ComputeResultantRing::resultant()`], unless you
114/// are implementing said method for a custom ring.
115/// 
116#[stability::unstable(feature = "enable")]
117pub fn resultant_finite_field<P>(ring: P, mut f: El<P>, mut g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
118    where P: RingStore + Copy,
119        P::Type: PolyRing + EuclideanRing,
120        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Domain + FiniteRing
121{
122    let base_ring = ring.base_ring();
123    if ring.is_zero(&g) || ring.is_zero(&f) {
124        return base_ring.zero();
125    }
126    let mut scale = base_ring.one();
127    if ring.degree(&g).unwrap() < ring.degree(&f).unwrap() {
128        if ring.degree(&g).unwrap() % 2 != 0 && ring.degree(&f).unwrap() % 2 != 0 {
129            base_ring.negate_inplace(&mut scale);
130        }
131        swap(&mut f, &mut g);
132    }
133
134    while ring.degree(&f).unwrap_or(0) >= 1 {
135        assert!(ring.degree(&g).unwrap() >= ring.degree(&f).unwrap());
136        let deg_g = ring.degree(&g).unwrap();
137        let r = ring.euclidean_rem(g, &f);
138        g = r;
139        base_ring.mul_assign(&mut scale, base_ring.pow(base_ring.clone_el(ring.lc(&f).unwrap()), deg_g - ring.degree(&g).unwrap_or(0)));
140
141        swap(&mut f, &mut g);
142        if ring.degree(&g).unwrap() % 2 != 0 && ring.degree(&f).unwrap_or(0) % 2 != 0 {
143            base_ring.negate_inplace(&mut scale);
144        }
145    }
146
147    if ring.is_zero(&f) {
148        return base_ring.zero();
149    } else {
150        let mut result = base_ring.clone_el(&ring.coefficient_at(&f, 0));
151        result = base_ring.pow(result, ring.degree(&g).unwrap());
152        base_ring.mul_assign(&mut result, scale);
153        return result;
154    }
155}
156
157///
158/// Trait for rings that support computing resultants of polynomials
159/// over the ring.
160/// 
161#[stability::unstable(feature = "enable")]
162pub trait ComputeResultantRing: RingBase {
163
164    ///
165    /// Computes the resultant of `f` and `g` over the base ring.
166    /// 
167    /// The resultant is the determinant of the linear map
168    /// ```text
169    ///   R[X]_deg(g)  x  R[X]_deg(f)  ->  R[X]_deg(fg),
170    ///        a       ,       b       ->    af + bg
171    /// ```
172    /// where `R[X]_d` refers to the vector space of polynomials in `R[X]` of degree
173    /// less than `d`.
174    /// 
175    /// # Example
176    /// ```rust
177    /// use feanor_math::assert_el_eq;
178    /// use feanor_math::ring::*;
179    /// use feanor_math::integer::*;
180    /// use feanor_math::rings::poly::dense_poly::DensePolyRing;
181    /// use feanor_math::rings::poly::*;
182    /// use feanor_math::algorithms::resultant::*;
183    /// let ZZ = BigIntRing::RING;
184    /// let ZZX = DensePolyRing::new(ZZ, "X");
185    /// let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1]);
186    /// // the discrimiant is the resultant of f and f'
187    /// let discriminant = <_ as ComputeResultantRing>::resultant(&ZZX, ZZX.clone_el(&f), derive_poly(&ZZX, &f));
188    /// assert_el_eq!(ZZ, ZZ.zero(), discriminant);
189    /// ```
190    /// 
191    fn resultant<P>(poly_ring: P, f: El<P>, g: El<P>) -> Self::Element
192        where P: RingStore + Copy,
193            P::Type: PolyRing,
194            <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>;
195}
196
197impl<R: ?Sized + EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso> ComputeResultantRing for R {
198
199    default fn resultant<P>(ring: P, f: El<P>, g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
200        where P: RingStore + Copy,
201            P::Type: PolyRing,
202            <P::Type as RingExtension>::BaseRing: RingStore<Type = R>
203    {
204        struct ComputeResultant<P>
205            where P: RingStore + Copy,
206                P::Type: PolyRing,
207                <<P::Type as RingExtension>::BaseRing as RingStore>::Type: EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso
208        {
209            ring: P,
210            f: El<P>,
211            g: El<P>
212        }
213        impl<P> FiniteRingOperation<<<P::Type as RingExtension>::BaseRing as RingStore>::Type> for ComputeResultant<P>
214            where P: RingStore + Copy,
215                P::Type: PolyRing,
216                <<P::Type as RingExtension>::BaseRing as RingStore>::Type: EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso
217        {
218            type Output = El<<P::Type as RingExtension>::BaseRing>;
219
220            fn execute(self) -> Self::Output
221                where <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing
222            {
223                let new_poly_ring = DensePolyRing::new(AsField::from(AsFieldBase::promise_is_perfect_field(self.ring.base_ring())), "X");
224                let hom = new_poly_ring.lifted_hom(&self.ring, WrapHom::new(new_poly_ring.base_ring().get_ring()));
225                let result = resultant_finite_field(&new_poly_ring, hom.map(self.f), hom.map(self.g));
226                return UnwrapHom::new(new_poly_ring.base_ring().get_ring()).map(result);
227            }
228
229            fn fallback(self) -> Self::Output {
230                let ring = self.ring;
231                let f = self.f;
232                let g = self.g;
233                let base_ring = ring.base_ring();
234                if ring.is_zero(&f) || ring.is_zero(&g) {
235                    return base_ring.zero();
236                }
237                let n = ring.degree(&f).unwrap() + ring.degree(&g).unwrap();
238                let coeff_bound_ln = ring.terms(&f).chain(ring.terms(&g)).map(|(c, _)| base_ring.get_ring().ln_pseudo_norm(c)).max_by(f64::total_cmp).unwrap();
239                let ln_max_norm = coeff_bound_ln * n as f64 + base_ring.get_ring().ln_pseudo_norm(&base_ring.int_hom().map(n as i32)) * n as f64 / 2.;
240
241                let work_locally = base_ring.get_ring().local_computation(ln_max_norm);
242                let mut resultants = Vec::new();
243                for i in 0..base_ring.get_ring().local_ring_count(&work_locally) {
244                    let embedding = EvaluatePolyLocallyReductionMap::new(base_ring.get_ring(), &work_locally, i);
245                    let new_poly_ring = DensePolyRing::new(embedding.codomain(), "X");
246                    let poly_ring_embedding = new_poly_ring.lifted_hom(ring, &embedding);
247                    let local_f = poly_ring_embedding.map_ref(&f);
248                    let local_g = poly_ring_embedding.map_ref(&g);
249                    let local_resultant = <_ as ComputeResultantRing>::resultant(&new_poly_ring, local_f, local_g);
250                    resultants.push(local_resultant);
251                }
252                return base_ring.get_ring().lift_combine(&work_locally, &resultants);
253            }
254        }
255        R::specialize(ComputeResultant { ring, f, g })
256    }
257}
258
259impl<I> ComputeResultantRing for RationalFieldBase<I>
260    where I: RingStore,
261        I::Type: IntegerRing 
262{
263    fn resultant<P>(ring: P, f: El<P>, g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
264        where P: RingStore + Copy,
265            P::Type: PolyRing,
266            <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>
267    {
268        if ring.is_zero(&g) || ring.is_zero(&f) {
269            return ring.base_ring().zero();
270        }
271        let QQ = ring.base_ring();
272        let ZZ = QQ.base_ring();
273        let f_deg = ring.degree(&f).unwrap();
274        let g_deg = ring.degree(&g).unwrap();
275        let f_den_lcm = ring.terms(&f).map(|(c, _)| ZZ.clone_el(QQ.get_ring().den(c))).fold(ZZ.one(), |a, b| signed_lcm(a, b, ZZ));
276        let g_den_lcm = ring.terms(&g).map(|(c, _)| ZZ.clone_el(QQ.get_ring().den(c))).fold(ZZ.one(), |a, b| signed_lcm(a, b, ZZ));
277        let ZZX = DensePolyRing::new(ZZ, "X");
278        let f_int = ZZX.from_terms(ring.terms(&f).map(|(c, i)| { let (a, b) = (QQ.get_ring().num(c), QQ.get_ring().den(c)); (ZZ.checked_div(&ZZ.mul_ref(&f_den_lcm, a), b).unwrap(), i) }));
279        let g_int = ZZX.from_terms(ring.terms(&g).map(|(c, i)| { let (a, b) = (QQ.get_ring().num(c), QQ.get_ring().den(c)); (ZZ.checked_div(&ZZ.mul_ref(&f_den_lcm, a), b).unwrap(), i) }));
280        let result_unscaled = <_ as ComputeResultantRing>::resultant(&ZZX, f_int, g_int);
281        return QQ.from_fraction(result_unscaled, ZZ.mul(ZZ.pow(f_den_lcm, g_deg), ZZ.pow(g_den_lcm, f_deg)));
282    }
283}
284
285#[cfg(test)]
286use crate::primitive_int::StaticRing;
287#[cfg(test)]
288use crate::rings::multivariate::multivariate_impl::MultivariatePolyRingImpl;
289#[cfg(test)]
290use crate::rings::multivariate::*;
291#[cfg(test)]
292use crate::algorithms::buchberger::buchberger_simple;
293#[cfg(test)]
294use crate::integer::BigIntRing;
295#[cfg(test)]
296use crate::algorithms::poly_gcd::PolyTFracGCDRing;
297#[cfg(test)]
298use crate::rings::zn::ZnRingStore;
299#[cfg(test)]
300use crate::rings::zn::zn_64::Zn;
301
302#[test]
303#[allow(deprecated)]
304fn test_resultant_global() {
305    let ZZ = StaticRing::<i64>::RING;
306    let ZZX = DensePolyRing::new(ZZ, "X");
307
308    // a quadratic polynomial and its derivative - the resultant should be the discriminant
309    let f = ZZX.from_terms([(3, 0), (-5, 1), (1, 2)].into_iter());
310    let g = ZZX.from_terms([(-5, 0), (2, 1)].into_iter());
311
312    assert_el_eq!(ZZ, -13, resultant_global(&ZZX, ZZX.clone_el(&f), ZZX.clone_el(&g)));
313    assert_el_eq!(ZZ, -13, resultant_global(&ZZX, g, f));
314
315    // if f and g have common factors, this should be zero
316    let f = ZZX.from_terms([(1, 0), (-2, 1), (1, 2)].into_iter());
317    let g = ZZX.from_terms([(-1, 0), (1, 2)].into_iter());
318    assert_el_eq!(ZZ, 0, resultant_global(&ZZX, f, g));
319
320    // a slightly larger example
321    let f = ZZX.from_terms([(5, 0), (-1, 1), (3, 2), (1, 4)].into_iter());
322    let g = ZZX.from_terms([(-1, 0), (-1, 2), (1, 3), (4, 5)].into_iter());
323    assert_el_eq!(ZZ, 642632, resultant_global(&ZZX, f, g));
324
325    let Fp = Zn::new(512409557603043077).as_field().ok().unwrap();
326    let FpX = DensePolyRing::new(Fp, "X");
327    let f = FpX.from_terms([(Fp.one(), 4), (Fp.int_hom().map(-2), 1), (Fp.int_hom().map(2), 0)].into_iter());
328    let g = FpX.from_terms([(Fp.one(), 64), (Fp.one(), 0)].into_iter());
329    assert_el_eq!(Fp, Fp.coerce(&StaticRing::<i64>::RING, 148105674794572153), resultant_global(&FpX, f, g));
330}
331
332#[test]
333#[allow(deprecated)]
334fn test_resultant_finite_field() {
335    let Fp = Zn::new(17).as_field().ok().unwrap();
336    let FpX = DensePolyRing::new(Fp, "X");
337
338    let [f, g] = FpX.with_wrapped_indeterminate(|X| [X.pow_ref(5) + 13 * X.pow_ref(4) + 4 * X.pow_ref(3) + X.pow_ref(2) + 14 * X + 7, X.pow_ref(3) + 1]);
339    assert_el_eq!(Fp, Fp.int_hom().map(8), resultant_finite_field(&FpX, FpX.clone_el(&f), FpX.clone_el(&g)));
340    assert_el_eq!(Fp, Fp.int_hom().map(9), resultant_finite_field(&FpX, FpX.clone_el(&g), FpX.clone_el(&f)));
341    assert_el_eq!(Fp, Fp.int_hom().map(8), resultant_global(&FpX, FpX.clone_el(&f), FpX.clone_el(&g)));
342    assert_el_eq!(Fp, Fp.int_hom().map(9), resultant_global(&FpX, FpX.clone_el(&g), FpX.clone_el(&f)));
343
344    let [f, g] = FpX.with_wrapped_indeterminate(|X| [X.pow_ref(4) + 4 * X.pow_ref(3) + X.pow_ref(2) + 14 * X + 7, X.pow_ref(3) + 1]);
345    assert_el_eq!(Fp, Fp.int_hom().map(5), resultant_finite_field(&FpX, FpX.clone_el(&f), FpX.clone_el(&g)));
346    assert_el_eq!(Fp, Fp.int_hom().map(5), resultant_finite_field(&FpX, FpX.clone_el(&g), FpX.clone_el(&f)));
347    assert_el_eq!(Fp, Fp.int_hom().map(5), resultant_global(&FpX, FpX.clone_el(&f), FpX.clone_el(&g)));
348    assert_el_eq!(Fp, Fp.int_hom().map(5), resultant_global(&FpX, FpX.clone_el(&g), FpX.clone_el(&f)));
349}
350
351#[test]
352#[allow(deprecated)]
353fn test_resultant_local_polynomial() {
354    let ZZ = BigIntRing::RING;
355    let QQ = RationalField::new(ZZ);
356    static_assert_impls!(RationalFieldBase<BigIntRing>: PolyTFracGCDRing);
357    // we eliminate `Y`, so add it as the outer indeterminate
358    let QQX = DensePolyRing::new(&QQ, "X");
359    let QQXY = DensePolyRing::new(&QQX, "Y");
360    let ZZ_to_QQ = QQ.int_hom();
361
362    // 1 + X^2 + 2 Y + (1 + X) Y^2
363    let f= QQXY.from_terms([
364        (vec![(1, 0), (1, 2)], 0),
365        (vec![(2, 0)], 1),
366        (vec![(1, 0), (1, 1)], 2)
367    ].into_iter().map(|(v, i)| (QQX.from_terms(v.into_iter().map(|(c, j)| (ZZ_to_QQ.map(c), j))), i)));
368
369    // 3 + X + (2 + X) Y + (1 + X + X^2) Y^2
370    let g = QQXY.from_terms([
371        (vec![(3, 0), (1, 1)], 0),
372        (vec![(2, 0), (1, 1)], 1),
373        (vec![(1, 0), (1, 1), (1, 2)], 2)
374    ].into_iter().map(|(v, i)| (QQX.from_terms(v.into_iter().map(|(c, j)| (ZZ_to_QQ.map(c), j))), i)));
375
376    let actual = QQX.normalize(resultant_global(&QQXY, QQXY.clone_el(&f), QQXY.clone_el(&g)));
377    let actual_local = <_ as ComputeResultantRing>::resultant(&QQXY, QQXY.clone_el(&f), QQXY.clone_el(&g));
378    assert_el_eq!(&QQX, &actual, &actual_local);
379
380    let QQYX = MultivariatePolyRingImpl::new(&QQ, 2);
381    // reverse the order of indeterminates, so that we indeed eliminate `Y`
382    let [f, g] = QQYX.with_wrapped_indeterminates(|[Y, X]| [ 1 + X.pow_ref(2) + 2 * Y + (1 + X) * Y.pow_ref(2), 3 + X + (2 + X) * Y + (1 + X + X.pow_ref(2)) * Y.pow_ref(2) ]);
383
384    let gb = buchberger_simple::<_, _>(&QQYX, vec![f, g], Lex);
385    let expected = gb.into_iter().filter(|poly| QQYX.appearing_indeterminates(&poly).len() == 1).collect::<Vec<_>>();
386    assert!(expected.len() == 1);
387    let expected = QQX.normalize(QQX.from_terms(QQYX.terms(&expected[0]).map(|(c, m)| (QQ.clone_el(c), QQYX.exponent_at(m, 1)))));
388
389    assert_el_eq!(QQX, expected, actual);
390}
391
392#[test]
393fn test_resultant_local_integer() {
394    let ZZ = BigIntRing::RING;
395    let ZZX = DensePolyRing::new(ZZ, "X");
396
397    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [
398        X.pow_ref(32) + 1,
399        X.pow_ref(2) - X - 1
400    ]);
401    assert_el_eq!(ZZ, ZZ.int_hom().map(4870849), <_ as ComputeResultantRing>::resultant(&ZZX, f, g));
402
403    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [
404        X.pow_ref(4) - 2 * X + 2,
405        X.pow_ref(64) + 1
406    ]);
407    assert_el_eq!(ZZ, ZZ.parse("381380816531458621441", 10).unwrap(), <_ as ComputeResultantRing>::resultant(&ZZX, f, g));
408
409    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [
410        X.pow_ref(32) - 2 * X.pow_ref(8) + 2,
411        X.pow_ref(512) + 1
412    ]);
413    assert_el_eq!(ZZ, ZZ.pow(ZZ.parse("381380816531458621441", 10).unwrap(), 8), <_ as ComputeResultantRing>::resultant(&ZZX, f, g));
414}
415
416#[test]
417#[ignore]
418fn test_resultant_large() {
419    let ZZ = BigIntRing::RING;
420    let ZZX = DensePolyRing::new(ZZ, "X");
421
422    let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 2 * X + 2]);
423    let g = ZZX.from_terms([(ZZ.one(), 1 << 14), (ZZ.one(), 0)]);
424    println!("start");
425    let result = BigIntRingBase::resultant(&ZZX, f, g);
426    println!("{}", ZZ.format(&result))
427}
feanor_math/algorithms/resultant.rs

feanor_math/algorithms/
resultant.rs
