feanor-math 3.5.18

use std::mem::swap;

use crate::algorithms;
use crate::algorithms::eea::signed_lcm;
use crate::delegate::{UnwrapHom, WrapHom};
use crate::divisibility::{DivisibilityRingStore, Domain};
use crate::homomorphism::*;
use crate::integer::*;
use crate::pid::{EuclideanRing, EuclideanRingStore, PrincipalIdealRing, PrincipalIdealRingStore};
use crate::reduce_lift::poly_eval::{EvalPolyLocallyRing, EvaluatePolyLocallyReductionMap};
use crate::ring::*;
use crate::rings::field::{AsField, AsFieldBase};
use crate::rings::finite::*;
use crate::rings::fraction::FractionFieldStore;
use crate::rings::poly::dense_poly::DensePolyRing;
use crate::rings::poly::*;
use crate::rings::rational::*;
use crate::specialization::FiniteRingOperation;

/// Computes the resultant of `f` and `g` over the base ring, using
/// a direct variant of Eulid's algorithm.
///
/// This function is deprecated, you should instead use
/// [`ComputeResultantRing::resultant()`] or [`resultant_finite_field()`].
///
/// The resultant is the determinant of the linear map
/// ```text
///   R[X]_deg(g)  x  R[X]_deg(f)  ->  R[X]_deg(fg),
///        a       ,       b       ->    af + bg
/// ```
/// where `R[X]_d` refers to the vector space of polynomials in `R[X]` of degree
/// less than `d`.
///
/// # Example
///
/// ```rust
/// use feanor_math::algorithms::resultant::*;
/// use feanor_math::primitive_int::*;
/// use feanor_math::ring::*;
/// use feanor_math::rings::poly::dense_poly::DensePolyRing;
/// use feanor_math::rings::poly::*;
/// let ZZ = StaticRing::<i64>::RING;
/// let ZZX = DensePolyRing::new(ZZ, "X");
/// let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1]);
/// // the discrimiant is the resultant of f and f'
/// let discriminant = resultant_global(&ZZX, ZZX.clone_el(&f), derive_poly(&ZZX, &f));
/// assert_eq!(0, discriminant);
/// ```
#[deprecated]
pub fn resultant_global<P>(ring: P, mut f: El<P>, mut g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
where
    P: RingStore + Copy,
    P::Type: PolyRing,
    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Domain + PrincipalIdealRing,
{
    let base_ring = ring.base_ring();
    if ring.is_zero(&g) || ring.is_zero(&f) {
        return base_ring.zero();
    }
    let mut scale_den = base_ring.one();
    let mut scale_num = base_ring.one();

    if ring.degree(&g).unwrap() < ring.degree(&f).unwrap() {
        if !ring.degree(&g).unwrap().is_multiple_of(2) && !ring.degree(&f).unwrap().is_multiple_of(2) {
            base_ring.negate_inplace(&mut scale_num);
        }
        std::mem::swap(&mut f, &mut g);
    }

    while ring.degree(&f).unwrap_or(0) >= 1 {
        let balance_factor = ring.get_ring().balance_poly(&mut f);
        if let Some(balance_factor) = balance_factor {
            base_ring.mul_assign(&mut scale_num, base_ring.pow(balance_factor, ring.degree(&g).unwrap()));
        }

        // 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)`
        let deg_g = ring.degree(&g).unwrap();
        let (_q, r, a) = algorithms::poly_div::poly_div_rem_domain(ring, g, &f);
        let deg_r = ring.degree(&r).unwrap_or(0);

        // adjust the scaling factor - we cancel out gcd's to prevent excessive number growth
        base_ring.mul_assign(&mut scale_den, base_ring.pow(a, ring.degree(&f).unwrap()));
        base_ring.mul_assign(
            &mut scale_num,
            base_ring.pow(base_ring.clone_el(ring.lc(&f).unwrap()), deg_g - deg_r),
        );
        let gcd = base_ring.ideal_gen(&scale_den, &scale_num);
        scale_den = base_ring.checked_div(&scale_den, &gcd).unwrap();
        scale_num = base_ring.checked_div(&scale_num, &gcd).unwrap();

        g = f;
        f = r;

        if !ring.degree(&g).unwrap().is_multiple_of(2) && !ring.degree(&f).unwrap_or(0).is_multiple_of(2) {
            base_ring.negate_inplace(&mut scale_num);
        }
    }

    if ring.is_zero(&f) {
        return base_ring.zero();
    } else {
        let mut result = base_ring.clone_el(ring.coefficient_at(&f, 0));
        result = base_ring.pow(result, ring.degree(&g).unwrap());
        base_ring.mul_assign(&mut result, scale_num);
        return base_ring.checked_div(&result, &scale_den).unwrap();
    }
}

/// Computes the resultant of two polynomials `f` and `g` over a finite field.
///
/// Usually you should use [`ComputeResultantRing::resultant()`], unless you
/// are implementing said method for a custom ring.
#[stability::unstable(feature = "enable")]
pub fn resultant_finite_field<P>(ring: P, mut f: El<P>, mut g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
where
    P: RingStore + Copy,
    P::Type: PolyRing + EuclideanRing,
    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Domain + FiniteRing,
{
    let base_ring = ring.base_ring();
    if ring.is_zero(&g) || ring.is_zero(&f) {
        return base_ring.zero();
    }
    let mut scale = base_ring.one();
    if ring.degree(&g).unwrap() < ring.degree(&f).unwrap() {
        if !ring.degree(&g).unwrap().is_multiple_of(2) && !ring.degree(&f).unwrap().is_multiple_of(2) {
            base_ring.negate_inplace(&mut scale);
        }
        swap(&mut f, &mut g);
    }

    while ring.degree(&f).unwrap_or(0) >= 1 {
        assert!(ring.degree(&g).unwrap() >= ring.degree(&f).unwrap());
        let deg_g = ring.degree(&g).unwrap();
        let r = ring.euclidean_rem(g, &f);
        g = r;
        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),
            ),
        );

        swap(&mut f, &mut g);
        if !ring.degree(&g).unwrap().is_multiple_of(2) && !ring.degree(&f).unwrap_or(0).is_multiple_of(2) {
            base_ring.negate_inplace(&mut scale);
        }
    }

    if ring.is_zero(&f) {
        return base_ring.zero();
    } else {
        let mut result = base_ring.clone_el(ring.coefficient_at(&f, 0));
        result = base_ring.pow(result, ring.degree(&g).unwrap());
        base_ring.mul_assign(&mut result, scale);
        return result;
    }
}

/// Trait for rings that support computing resultants of polynomials
/// over the ring.
#[stability::unstable(feature = "enable")]
pub trait ComputeResultantRing: RingBase {
    /// Computes the resultant of `f` and `g` over the base ring.
    ///
    /// The resultant is the determinant of the linear map
    /// ```text
    ///   R[X]_deg(g)  x  R[X]_deg(f)  ->  R[X]_deg(fg),
    ///        a       ,       b       ->    af + bg
    /// ```
    /// where `R[X]_d` refers to the vector space of polynomials in `R[X]` of degree
    /// less than `d`.
    ///
    /// # Example
    /// ```rust
    /// use feanor_math::algorithms::resultant::*;
    /// use feanor_math::assert_el_eq;
    /// use feanor_math::integer::*;
    /// use feanor_math::ring::*;
    /// use feanor_math::rings::poly::dense_poly::DensePolyRing;
    /// use feanor_math::rings::poly::*;
    /// let ZZ = BigIntRing::RING;
    /// let ZZX = DensePolyRing::new(ZZ, "X");
    /// let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1]);
    /// // the discrimiant is the resultant of f and f'
    /// let discriminant =
    ///     <_ as ComputeResultantRing>::resultant(&ZZX, ZZX.clone_el(&f), derive_poly(&ZZX, &f));
    /// assert_el_eq!(ZZ, ZZ.zero(), discriminant);
    /// ```
    fn resultant<P>(poly_ring: P, f: El<P>, g: El<P>) -> Self::Element
    where
        P: RingStore + Copy,
        P::Type: PolyRing,
        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>;
}

impl<R: ?Sized + EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso> ComputeResultantRing for R {
    default fn resultant<P>(ring: P, f: El<P>, g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
    where
        P: RingStore + Copy,
        P::Type: PolyRing,
        <P::Type as RingExtension>::BaseRing: RingStore<Type = R>,
    {
        struct ComputeResultant<P>
        where
            P: RingStore + Copy,
            P::Type: PolyRing,
            <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
                EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso,
        {
            ring: P,
            f: El<P>,
            g: El<P>,
        }
        impl<P> FiniteRingOperation<<<P::Type as RingExtension>::BaseRing as RingStore>::Type> for ComputeResultant<P>
        where
            P: RingStore + Copy,
            P::Type: PolyRing,
            <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
                EvalPolyLocallyRing + PrincipalIdealRing + Domain + SelfIso,
        {
            type Output = El<<P::Type as RingExtension>::BaseRing>;

            fn execute(self) -> Self::Output
            where
                <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing,
            {
                let new_poly_ring = DensePolyRing::new(
                    AsField::from(AsFieldBase::promise_is_perfect_field(self.ring.base_ring())),
                    "X",
                );
                let hom = new_poly_ring.lifted_hom(
                    &self.ring,
                    WrapHom::to_delegate_ring(new_poly_ring.base_ring().get_ring()),
                );
                let result = resultant_finite_field(&new_poly_ring, hom.map(self.f), hom.map(self.g));
                return UnwrapHom::from_delegate_ring(new_poly_ring.base_ring().get_ring()).map(result);
            }

            fn fallback(self) -> Self::Output {
                let ring = self.ring;
                let f = self.f;
                let g = self.g;
                let base_ring = ring.base_ring();
                if ring.is_zero(&f) || ring.is_zero(&g) {
                    return base_ring.zero();
                }
                let n = ring.degree(&f).unwrap() + ring.degree(&g).unwrap();
                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();
                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.0;

                let work_locally = base_ring.get_ring().local_computation(ln_max_norm);
                let mut resultants = Vec::new();
                for i in 0..base_ring.get_ring().local_ring_count(&work_locally) {
                    let embedding = EvaluatePolyLocallyReductionMap::new(base_ring.get_ring(), &work_locally, i);
                    let new_poly_ring = DensePolyRing::new(embedding.codomain(), "X");
                    let poly_ring_embedding = new_poly_ring.lifted_hom(ring, &embedding);
                    let local_f = poly_ring_embedding.map_ref(&f);
                    let local_g = poly_ring_embedding.map_ref(&g);
                    let local_resultant = <_ as ComputeResultantRing>::resultant(&new_poly_ring, local_f, local_g);
                    resultants.push(local_resultant);
                }
                return base_ring.get_ring().lift_combine(&work_locally, &resultants);
            }
        }
        R::specialize(ComputeResultant { ring, f, g })
    }
}

impl<I> ComputeResultantRing for RationalFieldBase<I>
where
    I: RingStore,
    I::Type: IntegerRing,
{
    fn resultant<P>(ring: P, f: El<P>, g: El<P>) -> El<<P::Type as RingExtension>::BaseRing>
    where
        P: RingStore + Copy,
        P::Type: PolyRing,
        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
    {
        if ring.is_zero(&g) || ring.is_zero(&f) {
            return ring.base_ring().zero();
        }
        let QQ = ring.base_ring();
        let ZZ = QQ.base_ring();
        let f_deg = ring.degree(&f).unwrap();
        let g_deg = ring.degree(&g).unwrap();
        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));
        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));
        let ZZX = DensePolyRing::new(ZZ, "X");
        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)
        }));
        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)
        }));
        let result_unscaled = <_ as ComputeResultantRing>::resultant(&ZZX, f_int, g_int);
        return QQ.from_fraction(
            result_unscaled,
            ZZ.mul(ZZ.pow(f_den_lcm, g_deg), ZZ.pow(g_den_lcm, f_deg)),
        );
    }
}

#[cfg(test)]
use crate::algorithms::buchberger::buchberger_simple;
#[cfg(test)]
use crate::algorithms::poly_gcd::PolyTFracGCDRing;
#[cfg(test)]
use crate::integer::BigIntRing;
#[cfg(test)]
use crate::primitive_int::StaticRing;
#[cfg(test)]
use crate::rings::multivariate::multivariate_impl::MultivariatePolyRingImpl;
#[cfg(test)]
use crate::rings::multivariate::*;
#[cfg(test)]
use crate::rings::zn::ZnRingStore;
#[cfg(test)]
use crate::rings::zn::zn_64::Zn;

#[test]
#[allow(deprecated)]
fn test_resultant_global() {
    let ZZ = StaticRing::<i64>::RING;
    let ZZX = DensePolyRing::new(ZZ, "X");

    // a quadratic polynomial and its derivative - the resultant should be the discriminant
    let f = ZZX.from_terms([(3, 0), (-5, 1), (1, 2)].into_iter());
    let g = ZZX.from_terms([(-5, 0), (2, 1)].into_iter());

    assert_el_eq!(ZZ, -13, resultant_global(&ZZX, ZZX.clone_el(&f), ZZX.clone_el(&g)));
    assert_el_eq!(ZZ, -13, resultant_global(&ZZX, g, f));

    // if f and g have common factors, this should be zero
    let f = ZZX.from_terms([(1, 0), (-2, 1), (1, 2)].into_iter());
    let g = ZZX.from_terms([(-1, 0), (1, 2)].into_iter());
    assert_el_eq!(ZZ, 0, resultant_global(&ZZX, f, g));

    // a slightly larger example
    let f = ZZX.from_terms([(5, 0), (-1, 1), (3, 2), (1, 4)].into_iter());
    let g = ZZX.from_terms([(-1, 0), (-1, 2), (1, 3), (4, 5)].into_iter());
    assert_el_eq!(ZZ, 642632, resultant_global(&ZZX, f, g));

    let Fp = Zn::new(512409557603043077).as_field().ok().unwrap();
    let FpX = DensePolyRing::new(Fp, "X");
    let f = FpX.from_terms([(Fp.one(), 4), (Fp.int_hom().map(-2), 1), (Fp.int_hom().map(2), 0)].into_iter());
    let g = FpX.from_terms([(Fp.one(), 64), (Fp.one(), 0)].into_iter());
    assert_el_eq!(
        Fp,
        Fp.coerce(&StaticRing::<i64>::RING, 148105674794572153),
        resultant_global(&FpX, f, g)
    );
}

#[test]
#[allow(deprecated)]
fn test_resultant_finite_field() {
    let Fp = Zn::new(17).as_field().ok().unwrap();
    let FpX = DensePolyRing::new(Fp, "X");

    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,
        ]
    });
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(8),
        resultant_finite_field(&FpX, FpX.clone_el(&f), FpX.clone_el(&g))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(9),
        resultant_finite_field(&FpX, FpX.clone_el(&g), FpX.clone_el(&f))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(8),
        resultant_global(&FpX, FpX.clone_el(&f), FpX.clone_el(&g))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(9),
        resultant_global(&FpX, FpX.clone_el(&g), FpX.clone_el(&f))
    );

    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,
        ]
    });
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(5),
        resultant_finite_field(&FpX, FpX.clone_el(&f), FpX.clone_el(&g))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(5),
        resultant_finite_field(&FpX, FpX.clone_el(&g), FpX.clone_el(&f))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(5),
        resultant_global(&FpX, FpX.clone_el(&f), FpX.clone_el(&g))
    );
    assert_el_eq!(
        Fp,
        Fp.int_hom().map(5),
        resultant_global(&FpX, FpX.clone_el(&g), FpX.clone_el(&f))
    );
}

#[test]
#[allow(deprecated)]
fn test_resultant_local_polynomial() {
    let ZZ = BigIntRing::RING;
    let QQ = RationalField::new(ZZ);
    static_assert_impls!(RationalFieldBase<BigIntRing>: PolyTFracGCDRing);
    // we eliminate `Y`, so add it as the outer indeterminate
    let QQX = DensePolyRing::new(&QQ, "X");
    let QQXY = DensePolyRing::new(&QQX, "Y");
    let ZZ_to_QQ = QQ.int_hom();

    // 1 + X^2 + 2 Y + (1 + X) Y^2
    let f = QQXY.from_terms(
        [(vec![(1, 0), (1, 2)], 0), (vec![(2, 0)], 1), (vec![(1, 0), (1, 1)], 2)]
            .into_iter()
            .map(|(v, i)| (QQX.from_terms(v.into_iter().map(|(c, j)| (ZZ_to_QQ.map(c), j))), i)),
    );

    // 3 + X + (2 + X) Y + (1 + X + X^2) Y^2
    let g = QQXY.from_terms(
        [
            (vec![(3, 0), (1, 1)], 0),
            (vec![(2, 0), (1, 1)], 1),
            (vec![(1, 0), (1, 1), (1, 2)], 2),
        ]
        .into_iter()
        .map(|(v, i)| (QQX.from_terms(v.into_iter().map(|(c, j)| (ZZ_to_QQ.map(c), j))), i)),
    );

    let actual = QQX.normalize(resultant_global(&QQXY, QQXY.clone_el(&f), QQXY.clone_el(&g)));
    let actual_local = <_ as ComputeResultantRing>::resultant(&QQXY, QQXY.clone_el(&f), QQXY.clone_el(&g));
    assert_el_eq!(&QQX, &actual, &actual_local);

    let QQYX = MultivariatePolyRingImpl::new(&QQ, 2);
    // reverse the order of indeterminates, so that we indeed eliminate `Y`
    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),
        ]
    });

    let gb = buchberger_simple::<_, _>(&QQYX, vec![f, g], Lex);
    let expected = gb
        .into_iter()
        .filter(|poly| QQYX.appearing_indeterminates(&poly).len() == 1)
        .collect::<Vec<_>>();
    assert!(expected.len() == 1);
    let expected = QQX.normalize(
        QQX.from_terms(
            QQYX.terms(&expected[0])
                .map(|(c, m)| (QQ.clone_el(c), QQYX.exponent_at(m, 1))),
        ),
    );

    assert_el_eq!(QQX, expected, actual);
}

#[test]
fn test_resultant_local_integer() {
    let ZZ = BigIntRing::RING;
    let ZZX = DensePolyRing::new(ZZ, "X");

    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(32) + 1, X.pow_ref(2) - X - 1]);
    assert_el_eq!(
        ZZ,
        ZZ.int_hom().map(4870849),
        <_ as ComputeResultantRing>::resultant(&ZZX, f, g)
    );

    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 2 * X + 2, X.pow_ref(64) + 1]);
    assert_el_eq!(
        ZZ,
        ZZ.parse("381380816531458621441", 10).unwrap(),
        <_ as ComputeResultantRing>::resultant(&ZZX, f, g)
    );

    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(32) - 2 * X.pow_ref(8) + 2, X.pow_ref(512) + 1]);
    assert_el_eq!(
        ZZ,
        ZZ.pow(ZZ.parse("381380816531458621441", 10).unwrap(), 8),
        <_ as ComputeResultantRing>::resultant(&ZZX, f, g)
    );
}

#[test]
#[ignore]
fn test_resultant_large() {
    let ZZ = BigIntRing::RING;
    let ZZX = DensePolyRing::new(ZZ, "X");

    let [f] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 2 * X + 2]);
    let g = ZZX.from_terms([(ZZ.one(), 1 << 14), (ZZ.one(), 0)]);
    println!("start");
    let result = BigIntRingBase::resultant(&ZZX, f, g);
    println!("{}", ZZ.format(&result))
}