Skip to main content

feanor_math/algorithms/
resultant.rs

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