feanor_math/algorithms/poly_factor/
extension.rs

1use crate::MAX_PROBABILISTIC_REPETITIONS;
2use crate::algorithms::poly_factor::FactorPolyField;
3use crate::algorithms::poly_gcd::PolyTFracGCDRing;
4use crate::algorithms::resultant::ComputeResultantRing;
5use crate::computation::*;
6use crate::field::*;
7use crate::homomorphism::*;
8use crate::integer::*;
9use crate::ordered::OrderedRingStore;
10use crate::pid::EuclideanRing;
11use crate::primitive_int::StaticRing;
12use crate::reduce_lift::poly_eval::InterpolationBaseRing;
13use crate::ring::*;
14use crate::rings::extension::{FreeAlgebra, FreeAlgebraStore};
15use crate::rings::poly::dense_poly::DensePolyRing;
16use crate::rings::poly::{PolyRing, PolyRingStore};
17use crate::specialization::FiniteRingSpecializable;
18
19#[stability::unstable(feature = "enable")]
20pub struct ProbablyNotSquarefree;
21
22/// Factors a polynomial with coefficients in a field `K` that is a simple, finite-degree
23/// field extension of a base field that supports polynomial factorization.
24///
25/// If the given polynomial is square-free, this will succeed, except with probability
26/// `2^-attempts`. If it is not square-free, this will always fail (i.e. return `Err`).
27#[stability::unstable(feature = "enable")]
28pub fn poly_factor_squarefree_extension<P, Controller>(
29    LX: P,
30    f: &El<P>,
31    attempts: usize,
32    controller: Controller,
33) -> Result<Vec<El<P>>, ProbablyNotSquarefree>
34where
35    P: RingStore,
36    P::Type: PolyRing + EuclideanRing,
37    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field + FreeAlgebra + PolyTFracGCDRing,
38    <<<<P::Type as RingExtension>::BaseRing as RingStore>::Type as RingExtension>::BaseRing as RingStore>::Type:
39        PerfectField + PolyTFracGCDRing + FactorPolyField + InterpolationBaseRing + FiniteRingSpecializable + SelfIso,
40    Controller: ComputationController,
41{
42    controller.run_computation(
43        format_args!(
44            "factor_ring_ext(pdeg={}, extdeg={})",
45            LX.degree(f).unwrap(),
46            LX.base_ring().rank()
47        ),
48        |controller| {
49            let L = LX.base_ring();
50            let K = L.base_ring();
51
52            assert!(LX.base_ring().is_one(LX.lc(f).unwrap()));
53
54            let KX = DensePolyRing::new(K, "X");
55            // Y will take the role of the ring generator `theta`, and `X` remains the indeterminate
56            let KXY = DensePolyRing::new(KX.clone(), "Y");
57
58            let Norm = |f: El<P>| {
59                let f_over_KY = <_ as RingStore>::sum(
60                    &KXY,
61                    LX.terms(&f).map(|(c, i)| {
62                        let mut result = L.poly_repr(&KXY, c, KX.inclusion());
63                        KXY.inclusion()
64                            .mul_assign_map(&mut result, KX.pow(KX.indeterminate(), i));
65                        result
66                    }),
67                );
68                let gen_poly = L.generating_poly(&KXY, KX.inclusion());
69
70                return <_ as ComputeResultantRing>::resultant(&KXY, f_over_KY, gen_poly);
71            };
72
73            // we want to find `k` such that `N(f(X + k theta))` remains square-free, where `theta`
74            // generates `L`; there are only finitely many `k` that don't work, by the
75            // following idea: Write `f = f1 ... fr` for the factorization of `f`, where
76            // `r <= d`. Then `N(f(X + k theta))` is not square-free, iff there exist `i
77            // != j` and `sigma` such that `fi(X + k theta) = sigma(fj(X + k theta)) = sigma(fj)(X +
78            // k sigma(theta))`. Now note that for given `i, j, sigma`, `fi(X + k theta)
79            // = sigma(fj)(X + k sigma(theta))` for at most one `k`, as otherwise
80            // `fi(a (k1 - k2) theta) sigma(theta)) = sigma(fj)(a (k1 - k2) sigma(theta))` for all
81            // `a in Z` (evaluation is ring hom).
82            //
83            // Now note that `fi(X) - sigma(fj(X))` is a linear combination of `X^m` and
84            // `sigma(X^m)`, i.e. of `deg(fi) + deg(fj) <= d` basic functions. The
85            // vectors `(a1 theta)^m, ..., (al theta)^m` for `m <= deg(fi)` and `sigma(a1 theta)^m,
86            // ..., sigma(al theta)^m` for `m <= deg(fj)` are all linearly independent
87            // (assuming `a1, ..., al` distinct and `l >= deg(fi) + deg(fj) + 2`),
88            // thus `char(K) >= d + 2` (or `char(K) = 0`) implies that `fi = fj = 0`, a
89            // contradiction.
90            //
91            // Thus we find that there are at most `d(d + 1)/2 * [L : K]` many `k` such that `N(f(X
92            // + k theta))` is not squarefree. This would even work if `k` is not an
93            // integer, but any element of `K`. Note that we still require `char(K) >= d + 2` for
94            // the previous paragraph to work.
95
96            let ZZbig = BigIntRing::RING;
97            let characteristic = K.characteristic(&ZZbig).unwrap();
98            // choose bound about twice as large as necessary, so the probability of succeeding is
99            // almost 1/2
100            let bound = LX.degree(f).unwrap() * LX.degree(f).unwrap() * L.rank();
101            assert!(
102                ZZbig.is_zero(&characteristic)
103                    || ZZbig.is_geq(
104                        &characteristic,
105                        &int_cast(bound.try_into().unwrap(), ZZbig, StaticRing::<i64>::RING)
106                    )
107            );
108
109            let mut rng = oorandom::Rand64::new(1);
110
111            for _ in 0..attempts {
112                let k = StaticRing::<i32>::RING.get_uniformly_random(&bound.try_into().unwrap(), || rng.rand_u64());
113                let lin_transform = LX.from_terms([(L.mul(L.canonical_gen(), L.int_hom().map(k)), 0), (L.one(), 1)]);
114                let f_transformed = LX.evaluate(f, &lin_transform, LX.inclusion());
115
116                let norm_f_transformed = Norm(LX.clone_el(&f_transformed));
117                log_progress!(controller, "(norm)");
118                let degree = KX.degree(&norm_f_transformed).unwrap();
119                let squarefree_part = <_ as PolyTFracGCDRing>::squarefree_part(&KX, &norm_f_transformed);
120                log_progress!(controller, "(squarefree)");
121
122                if KX.degree(&squarefree_part).unwrap() == degree {
123                    let lin_transform_rev =
124                        LX.from_terms([(L.mul(L.canonical_gen(), L.int_hom().map(-k)), 0), (L.one(), 1)]);
125                    let (factorization, _unit) =
126                        <_ as FactorPolyField>::factor_poly_with_controller(&KX, &squarefree_part, controller.clone());
127                    log_progress!(controller, "(factored)");
128
129                    return Ok(factorization
130                        .into_iter()
131                        .map(|(factor, e)| {
132                            assert!(e == 1);
133                            let f_factor = LX.normalize(<_ as PolyTFracGCDRing>::gcd(
134                                &LX,
135                                &f_transformed,
136                                &LX.lifted_hom(&KX, L.inclusion()).map(factor),
137                            ));
138                            return LX.evaluate(&f_factor, &lin_transform_rev, LX.inclusion());
139                        })
140                        .collect());
141                }
142            }
143            return Err(ProbablyNotSquarefree);
144        },
145    )
146}
147
148/// Factors a polynomial with coefficients in a field `K` that is a simple, finite-degree
149/// field extension of a base field that supports polynomial factorization.
150#[stability::unstable(feature = "enable")]
151pub fn poly_factor_extension<P, Controller>(
152    poly_ring: P,
153    f: &El<P>,
154    controller: Controller,
155) -> (Vec<(El<P>, usize)>, El<<P::Type as RingExtension>::BaseRing>)
156where
157    P: RingStore,
158    P::Type: PolyRing + EuclideanRing,
159    <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
160        FreeAlgebra + PerfectField + FiniteRingSpecializable + PolyTFracGCDRing,
161    <<<<P::Type as RingExtension>::BaseRing as RingStore>::Type as RingExtension>::BaseRing as RingStore>::Type:
162        PerfectField + PolyTFracGCDRing + FactorPolyField + InterpolationBaseRing + FiniteRingSpecializable + SelfIso,
163    Controller: ComputationController,
164{
165    let KX = &poly_ring;
166    let K = KX.base_ring();
167
168    // We use the approach outlined in Cohen's "a course in computational algebraic number theory".
169    //  - Use square-free reduction to assume wlog that `f` is square-free
170    //  - Observe that the factorization of `f` is the product over `gcd(f, g)` where `g` runs through
171    //    the factors of `N(f)` over `QQ[X]` - assuming that `N(f)` is square-free! Here `N(f)` is the
172    //    "norm" of `f`, i.e. the product `prod_sigma sigma(f)` where `sigma` runs through the
173    //    embeddings `K -> CC`.
174    //  - It is now left to actually compute `N(f)`, which is not so simple as we do not known the
175    //    `sigma`. As it turns out, this is the resultant of `f` and `MiPo(theta)` where `theta`
176    //    generates `K`
177
178    assert!(!KX.is_zero(f));
179    let mut result: Vec<(El<P>, usize)> = Vec::new();
180    for (non_irred_factor, k) in <_ as PolyTFracGCDRing>::power_decomposition(KX, f) {
181        for factor in
182            poly_factor_squarefree_extension(KX, &non_irred_factor, MAX_PROBABILISTIC_REPETITIONS, controller.clone())
183                .ok()
184                .unwrap()
185        {
186            if let Some((i, _)) = result.iter().enumerate().find(|(_, f)| KX.eq_el(&f.0, &factor)) {
187                result[i].1 += k;
188            } else {
189                result.push((factor, k));
190            }
191        }
192    }
193    return (result, K.clone_el(KX.coefficient_at(f, 0)));
194}
feanor_math/algorithms/poly_factor/extension.rs

feanor_math/algorithms/poly_factor/
extension.rs
