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#[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 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 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#[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#[stability::unstable(feature = "enable")]
165pub trait ComputeResultantRing: RingBase {
166 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 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 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 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 let QQX = DensePolyRing::new(&QQ, "X");
441 let QQXY = DensePolyRing::new(&QQX, "Y");
442 let ZZ_to_QQ = QQ.int_hom();
443
444 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 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 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}