feanor-math 3.5.18

use crate::MAX_PROBABILISTIC_REPETITIONS;
use crate::algorithms::poly_factor::FactorPolyField;
use crate::algorithms::poly_gcd::PolyTFracGCDRing;
use crate::algorithms::resultant::ComputeResultantRing;
use crate::computation::*;
use crate::field::*;
use crate::homomorphism::*;
use crate::integer::*;
use crate::ordered::OrderedRingStore;
use crate::pid::EuclideanRing;
use crate::primitive_int::StaticRing;
use crate::reduce_lift::poly_eval::InterpolationBaseRing;
use crate::ring::*;
use crate::rings::extension::{FreeAlgebra, FreeAlgebraStore};
use crate::rings::poly::dense_poly::DensePolyRing;
use crate::rings::poly::{PolyRing, PolyRingStore};
use crate::specialization::FiniteRingSpecializable;

#[stability::unstable(feature = "enable")]
pub struct ProbablyNotSquarefree;

/// Factors a polynomial with coefficients in a field `K` that is a simple, finite-degree
/// field extension of a base field that supports polynomial factorization.
///
/// If the given polynomial is square-free, this will succeed, except with probability
/// `2^-attempts`. If it is not square-free, this will always fail (i.e. return `Err`).
#[stability::unstable(feature = "enable")]
pub fn poly_factor_squarefree_extension<P, Controller>(
    LX: P,
    f: &El<P>,
    attempts: usize,
    controller: Controller,
) -> Result<Vec<El<P>>, ProbablyNotSquarefree>
where
    P: RingStore,
    P::Type: PolyRing + EuclideanRing,
    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field + FreeAlgebra + PolyTFracGCDRing,
    <<<<P::Type as RingExtension>::BaseRing as RingStore>::Type as RingExtension>::BaseRing as RingStore>::Type:
        PerfectField + PolyTFracGCDRing + FactorPolyField + InterpolationBaseRing + FiniteRingSpecializable + SelfIso,
    Controller: ComputationController,
{
    controller.run_computation(
        format_args!(
            "factor_ring_ext(pdeg={}, extdeg={})",
            LX.degree(f).unwrap(),
            LX.base_ring().rank()
        ),
        |controller| {
            let L = LX.base_ring();
            let K = L.base_ring();

            assert!(LX.base_ring().is_one(LX.lc(f).unwrap()));

            let KX = DensePolyRing::new(K, "X");
            // Y will take the role of the ring generator `theta`, and `X` remains the indeterminate
            let KXY = DensePolyRing::new(KX.clone(), "Y");

            let Norm = |f: El<P>| {
                let f_over_KY = <_ as RingStore>::sum(
                    &KXY,
                    LX.terms(&f).map(|(c, i)| {
                        let mut result = L.poly_repr(&KXY, c, KX.inclusion());
                        KXY.inclusion()
                            .mul_assign_map(&mut result, KX.pow(KX.indeterminate(), i));
                        result
                    }),
                );
                let gen_poly = L.generating_poly(&KXY, KX.inclusion());

                return <_ as ComputeResultantRing>::resultant(&KXY, f_over_KY, gen_poly);
            };

            // we want to find `k` such that `N(f(X + k theta))` remains square-free, where `theta`
            // generates `L`; there are only finitely many `k` that don't work, by the
            // following idea: Write `f = f1 ... fr` for the factorization of `f`, where
            // `r <= d`. Then `N(f(X + k theta))` is not square-free, iff there exist `i
            // != j` and `sigma` such that `fi(X + k theta) = sigma(fj(X + k theta)) = sigma(fj)(X +
            // k sigma(theta))`. Now note that for given `i, j, sigma`, `fi(X + k theta)
            // = sigma(fj)(X + k sigma(theta))` for at most one `k`, as otherwise
            // `fi(a (k1 - k2) theta) sigma(theta)) = sigma(fj)(a (k1 - k2) sigma(theta))` for all
            // `a in Z` (evaluation is ring hom).
            //
            // Now note that `fi(X) - sigma(fj(X))` is a linear combination of `X^m` and
            // `sigma(X^m)`, i.e. of `deg(fi) + deg(fj) <= d` basic functions. The
            // vectors `(a1 theta)^m, ..., (al theta)^m` for `m <= deg(fi)` and `sigma(a1 theta)^m,
            // ..., sigma(al theta)^m` for `m <= deg(fj)` are all linearly independent
            // (assuming `a1, ..., al` distinct and `l >= deg(fi) + deg(fj) + 2`),
            // thus `char(K) >= d + 2` (or `char(K) = 0`) implies that `fi = fj = 0`, a
            // contradiction.
            //
            // Thus we find that there are at most `d(d + 1)/2 * [L : K]` many `k` such that `N(f(X
            // + k theta))` is not squarefree. This would even work if `k` is not an
            // integer, but any element of `K`. Note that we still require `char(K) >= d + 2` for
            // the previous paragraph to work.

            let ZZbig = BigIntRing::RING;
            let characteristic = K.characteristic(&ZZbig).unwrap();
            // choose bound about twice as large as necessary, so the probability of succeeding is
            // almost 1/2
            let bound = LX.degree(f).unwrap() * LX.degree(f).unwrap() * L.rank();
            assert!(
                ZZbig.is_zero(&characteristic)
                    || ZZbig.is_geq(
                        &characteristic,
                        &int_cast(bound.try_into().unwrap(), ZZbig, StaticRing::<i64>::RING)
                    )
            );

            let mut rng = oorandom::Rand64::new(1);

            for _ in 0..attempts {
                let k = StaticRing::<i32>::RING.get_uniformly_random(&bound.try_into().unwrap(), || rng.rand_u64());
                let lin_transform = LX.from_terms([(L.mul(L.canonical_gen(), L.int_hom().map(k)), 0), (L.one(), 1)]);
                let f_transformed = LX.evaluate(f, &lin_transform, LX.inclusion());

                let norm_f_transformed = Norm(LX.clone_el(&f_transformed));
                log_progress!(controller, "(norm)");
                let degree = KX.degree(&norm_f_transformed).unwrap();
                let squarefree_part = <_ as PolyTFracGCDRing>::squarefree_part(&KX, &norm_f_transformed);
                log_progress!(controller, "(squarefree)");

                if KX.degree(&squarefree_part).unwrap() == degree {
                    let lin_transform_rev =
                        LX.from_terms([(L.mul(L.canonical_gen(), L.int_hom().map(-k)), 0), (L.one(), 1)]);
                    let (factorization, _unit) =
                        <_ as FactorPolyField>::factor_poly_with_controller(&KX, &squarefree_part, controller.clone());
                    log_progress!(controller, "(factored)");

                    return Ok(factorization
                        .into_iter()
                        .map(|(factor, e)| {
                            assert!(e == 1);
                            let f_factor = LX.normalize(<_ as PolyTFracGCDRing>::gcd(
                                &LX,
                                &f_transformed,
                                &LX.lifted_hom(&KX, L.inclusion()).map(factor),
                            ));
                            return LX.evaluate(&f_factor, &lin_transform_rev, LX.inclusion());
                        })
                        .collect());
                }
            }
            return Err(ProbablyNotSquarefree);
        },
    )
}

/// Factors a polynomial with coefficients in a field `K` that is a simple, finite-degree
/// field extension of a base field that supports polynomial factorization.
#[stability::unstable(feature = "enable")]
pub fn poly_factor_extension<P, Controller>(
    poly_ring: P,
    f: &El<P>,
    controller: Controller,
) -> (Vec<(El<P>, usize)>, El<<P::Type as RingExtension>::BaseRing>)
where
    P: RingStore,
    P::Type: PolyRing + EuclideanRing,
    <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
        FreeAlgebra + PerfectField + FiniteRingSpecializable + PolyTFracGCDRing,
    <<<<P::Type as RingExtension>::BaseRing as RingStore>::Type as RingExtension>::BaseRing as RingStore>::Type:
        PerfectField + PolyTFracGCDRing + FactorPolyField + InterpolationBaseRing + FiniteRingSpecializable + SelfIso,
    Controller: ComputationController,
{
    let KX = &poly_ring;
    let K = KX.base_ring();

    // We use the approach outlined in Cohen's "a course in computational algebraic number theory".
    //  - Use square-free reduction to assume wlog that `f` is square-free
    //  - Observe that the factorization of `f` is the product over `gcd(f, g)` where `g` runs through
    //    the factors of `N(f)` over `QQ[X]` - assuming that `N(f)` is square-free! Here `N(f)` is the
    //    "norm" of `f`, i.e. the product `prod_sigma sigma(f)` where `sigma` runs through the
    //    embeddings `K -> CC`.
    //  - It is now left to actually compute `N(f)`, which is not so simple as we do not known the
    //    `sigma`. As it turns out, this is the resultant of `f` and `MiPo(theta)` where `theta`
    //    generates `K`

    assert!(!KX.is_zero(f));
    let mut result: Vec<(El<P>, usize)> = Vec::new();
    for (non_irred_factor, k) in <_ as PolyTFracGCDRing>::power_decomposition(KX, f) {
        for factor in
            poly_factor_squarefree_extension(KX, &non_irred_factor, MAX_PROBABILISTIC_REPETITIONS, controller.clone())
                .ok()
                .unwrap()
        {
            if let Some((i, _)) = result.iter().enumerate().find(|(_, f)| KX.eq_el(&f.0, &factor)) {
                result[i].1 += k;
            } else {
                result.push((factor, k));
            }
        }
    }
    return (result, K.clone_el(KX.coefficient_at(f, 0)));
}