Overhaul preprocessing for radical computations.

experimental/GaloisGrp/src/Solve.jl

@@ -236,14 +236,14 @@ function _fixed_field(C::GaloisCtx, S::SubField, U::PermGroup; invar=nothing, ma
   return SS
 end
 
-function Oscar.extension_field(f::AbstractAlgebra.Generic.Poly{QQPolyRingElem}; cached::Bool, check::Bool)
+function Oscar.extension_field(f::AbstractAlgebra.Generic.Poly{QQPolyRingElem}; cached::Bool=false, check::Bool=true)
   C = base_ring(f)
   Qt, t = rational_function_field(QQ, symbols(C)[1], cached = false)
   ff = map_coefficients(x->x(t), f)
   return extension_field(ff, cached = cached, check = check)
 end
 
-function Oscar.extension_field(f::AbstractAlgebra.Generic.Poly{<:NumFieldElem}; cached::Bool, check::Bool)
+function Oscar.extension_field(f::AbstractAlgebra.Generic.Poly{<:NumFieldElem}; cached::Bool=false, check::Bool=true)
   return number_field(f; cached, check)
 end
 

src/AlgebraicGeometry/Schemes/AffineSchemes/Objects/Attributes.jl

@@ -615,7 +615,7 @@ julia> singular_locus(A3)
 (scheme(1), Hom: scheme(1) -> affine 3-space)
 
 julia> singular_locus(X)
-(scheme(x^2 - y^2 + z^2, 2*x, -2*y, 2*z), Hom: scheme(x^2 - y^2 + z^2, 2*x, -2*y, 2*z) -> scheme(x^2 - y^2 + z^2))
+(scheme(x^2 - y^2 + z^2, x, y, z), Hom: scheme(x^2 - y^2 + z^2, x, y, z) -> scheme(x^2 - y^2 + z^2))
 
 julia> U = complement_of_point_ideal(R, [0,0,0])
 Complement

src/AlgebraicGeometry/Schemes/CoveredSchemes/Morphisms/Attributes.jl

@@ -45,7 +45,7 @@ julia> I, s = singular_locus(Ycov)
 julia> covering_morphism(s)
 Covering morphism
   from covering with 1 patch
-    1a: [(x//z), (y//z)]   scheme((x//z)^3 - (y//z)^2, (y//z), (x//z))
+    1a: [(x//z), (y//z)]   scheme((x//z)^3 - (y//z)^2, (x//z), (y//z))
   to covering with 3 patches
     1b: [(y//x), (z//x)]   scheme(-(y//x)^2*(z//x) + 1)
     2b: [(x//y), (z//y)]   scheme((x//y)^3 - (z//y))

src/AlgebraicGeometry/Schemes/CoveredSchemes/Objects/Attributes.jl

@@ -238,14 +238,14 @@ Scheme
   over rational field
 with default covering
   described by patches
-    1: scheme((x//z)^3 - (y//z)^2, (y//z), (x//z))
+    1: scheme((x//z)^3 - (y//z)^2, (x//z), (y//z))
   in the coordinate(s)
     1: [(x//z), (y//z)]
 
 julia> s
 Covered scheme morphism
   from scheme over QQ covered with 1 patch
-    1a: [(x//z), (y//z)]   scheme((x//z)^3 - (y//z)^2, (y//z), (x//z))
+    1a: [(x//z), (y//z)]   scheme((x//z)^3 - (y//z)^2, (x//z), (y//z))
   to scheme over QQ covered with 3 patches
     1b: [(y//x), (z//x)]   scheme(-(y//x)^2*(z//x) + 1)
     2b: [(x//y), (z//y)]   scheme((x//y)^3 - (z//y))

src/AlgebraicGeometry/Schemes/Sheaves/Types.jl

@@ -775,8 +775,9 @@ Ideal generated by
 
 julia> JJ(U)
 Ideal generated by
-  (z//x)
   (y//x)
+  (y//x)*(z//x)
+  (z//x)
 
 ```
 """
@@ -899,7 +900,7 @@ Sheaf of ideals
 with restrictions
   1: Ideal (1)
   2: Ideal (1)
-  3: Ideal (x_2_3, x_1_3)
+  3: Ideal (x_1_3, x_2_3)
 
 julia> typeof(JJ)
 Oscar.SingularLocusIdealSheaf{CoveredScheme{QQField}, AbsAffineScheme, Ideal, Map}

src/Rings/MPolyQuo.jl

@@ -554,11 +554,15 @@ end
   return is_prime(saturated_ideal(I))
 end
 
-function radical(I::MPolyQuoIdeal; eliminate_variables::Bool=true)
+function radical(
+    I::MPolyQuoIdeal; 
+    eliminate_variables::Bool=true, 
+    use_squarefree_parts_of_generators::Bool=true
+  )
   get_attribute!(I, :radical) do
     R = base_ring(I)
     J = saturated_ideal(I)
-    return ideal(R, [g for g in R.(gens(radical(J; eliminate_variables))) if !iszero(g)])
+    return ideal(R, [g for g in R.(gens(radical(J; eliminate_variables, use_squarefree_parts_of_generators))) if !iszero(g)])
   end::typeof(I)
 end
 

src/Rings/mpoly-graded.jl

@@ -879,17 +879,16 @@ end
 
 function factor(x::MPolyDecRingElem)
   R = parent(x)
-  D = Dict{elem_type(R), Int64}()
   F = factor(forget_decoration(x))
-  n=length(F.fac)
-  #if n == 1
-  #  return Fac(R(F.unit), D)
-  #else
-    for i in keys(F.fac)
-     push!(D, R(i) => Int64(F[i]))
-    end
-  return Fac(R(F.unit), D)
-  #end
+  D = Dict{elem_type(R), Int}(R(i) => e for (i, e) in F)
+  return Fac(R(unit(F)), D)
+end
+
+function factor_squarefree(x::MPolyDecRingElem)
+  R = parent(x)
+  F = factor_squarefree(forget_decoration(x))
+  D = Dict{elem_type(R), Int}(R(i) => e for (i, e) in F)
+  return Fac(R(unit(F)), D)
 end
 
 function gcd(x::MPolyDecRingElem, y::MPolyDecRingElem)

src/Rings/mpoly-ideals.jl

@@ -475,7 +475,7 @@ end
 
 #######################################################
 @doc raw"""
-    radical(I::MPolyIdeal{T}; eliminate_variables::Bool = true) where {T <: MPolyRingElem}
+    radical(I::MPolyIdeal{T}; eliminate_variables::Bool = true, use_squarefree_parts_of_generators::Bool = true) where {T <: MPolyRingElem}
 
 Return the radical of `I`.
 
@@ -487,6 +487,8 @@ Return the radical of `I`.
 
 When `eliminate_variables` is set to `true`, a preprocessing heuristic is applied in order to find variables which can be eliminated using `I`.
 
+When `use_squarefree_parts_of_generators` is set to `true`, the squarefree part of all generators is computed before passing to the actual radical computation. 
+
 # Examples
 ```jldoctest
 julia> R, (x, y) = polynomial_ring(QQ, [:x, :y])
@@ -575,42 +577,31 @@ Ideal generated by
 
 julia> RI = radical(I)
 Ideal generated by
+  1326*a*d
+  663*a*c
   102*b*d
-  78*a*d
   51*b*c
+  663*a*c*d
+  51*b*c*d
+  78*a*d
   39*a*c
   6*a*b*d
   3*a*b*c
-  1326*a^2*d^5
-  1989*a^2*c^5
-  102*b^4*d^5
-  153*b^4*c^5
-  663*a^2*c^5*d^5
-  51*b^4*c^5*d^5
-  78*a^2*d^15
-  117*a^2*c^15
-  78*a^15*d^5
-  117*a^15*c^5
-  6*a^2*b^4*d^15
-  9*a^2*b^4*c^15
-  39*a^2*c^5*d^15
-  39*a^2*c^15*d^5
-  6*a^2*b^15*d^5
-  9*a^2*b^15*c^5
-  6*a^15*b^4*d^5
-  9*a^15*b^4*c^5
-  39*a^15*c^5*d^5
-  3*a^2*b^4*c^5*d^15
-  3*a^2*b^4*c^15*d^5
-  3*a^2*b^15*c^5*d^5
-  3*a^15*b^4*c^5*d^5
+  39*a*c*d
+  3*a*b*c*d
 ```
 """
 function radical(
     I::MPolyIdeal{T}; 
-    eliminate_variables::Bool=true
+    eliminate_variables::Bool=true,
+    use_squarefree_parts_of_generators::Bool=true
   ) where {T <: MPolyRingElem}
   get_attribute!(I, :radical) do
+    # We first check for the squarefree flag. 
+    # Elimination of variables will probably affect how the generators factor. 
+    if use_squarefree_parts_of_generators
+      return radical(_squarefree_generator_ideal(I); eliminate_variables, use_squarefree_parts_of_generators=false)
+    end
     if eliminate_variables
       is_known_to_be_radical(I) && return I
       # Calling `elimpart` (within `simplify`) turns out to significantly speed things up in many cases.
@@ -620,7 +611,7 @@ function radical(
       is_zero(ngens(S)) && return I # This only happens when all variables can be eliminated. 
       # Then the ideal defines a reduced point.
       J = modulus(S)
-      pre_res = _compute_radical(J)
+      pre_res = radical(J; eliminate_variables=false, use_squarefree_parts_of_generators)
       pre_res = ideal(R, elem_type(R)[lift(iso_inv(S(g))) for g in gens(pre_res)])
       res = ideal(R, small_generating_set(pre_res + I))
       set_attribute!(res, :is_radical=>true)
@@ -650,39 +641,43 @@ function _compute_radical(I::T) where {T <: MPolyIdeal}
   return Irad
 end
 
+# Go through the generators of I and replace each generator by its squarefree part. 
+# This is a heuristic preprocessing method for radical computations.
+function _squarefree_generator_ideal(I::MPolyIdeal)
+  R = base_ring(I)
+  res_gens = elem_type(R)[]
+  for g in gens(I)
+    is_zero(g) && continue
+    # In characteristic zero, there are special methods, but not in pos. char. 
+    # See e.g. Geddes et. al. "Algorithms for Computer Algebra", Springer (1992)
+    fact = is_zero(characteristic(R)) ? factor_squarefree(g) : factor(g)
+    is_empty(fact) && return ideal(R, one(R))
+    h = prod(x for (x, _) in fact; init=one(R))
+    push!(res_gens, h)
+  end
+  return ideal(R, res_gens)
+end
+
 # Rerouting via expansion of the coefficient field
 function radical(
     I::MPolyIdeal{T};
-    factor_generators::Bool=true, 
+    use_squarefree_parts_of_generators::Bool=true, 
     eliminate_variables::Bool=true
   ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
   is_known_to_be_radical(I) && return I
   get_attribute!(I, :radical) do
     is_one(I) && return I
     R = base_ring(I)
-    J = ideal(R, zero(R))
-    if factor_generators
-      # In practice this will often lead to significant speedup due to reduction of degrees.
-      # TODO:  is the following faster? radical(ab,c) = intersect(radical(a,c),radical(b,c)) 
-      for g in gens(I)
-        is_zero(g) && continue
-        fact = factor(g)
-        is_empty(fact) && continue
-        h = one(g)
-        for (x, k) in fact
-          h = h*x
-        end
-        J = J + ideal(R, h)
-      end
-    else
-      J = I
-    end
+
+    use_squarefree_parts_of_generators && return radical(_squarefree_generator_ideal(I); eliminate_variables, use_squarefree_parts_of_generators=false)
+
     R_flat, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
-    I_flat = ideal(R_flat, iso_inv.(gens(J)))
-    I_flat_rad = radical(I_flat; eliminate_variables)
+    I_flat = ideal(R_flat, iso_inv.(gens(I)))
+    I_flat_rad = radical(I_flat; 
+                         eliminate_variables, 
+                         use_squarefree_parts_of_generators=false) # This only pays off on the top level. 
     Irad = iso(I_flat_rad)
     set_attribute!(Irad, :is_radical => true)
-    @hassert :IdealSheaves 2 !is_one(Irad)
     Irad
   end::MPolyIdeal{T}
 end
@@ -891,9 +886,6 @@ defined over a number field of degree $d_{ij}$ whose generator prints as `_a`.
 The implementation combines the algorithm of Gianni, Trager, and Zacharias for primary
 decomposition with absolute polynomial factorization.
 
-!!! warning
-    Over number fields this proceduce might return redundant output.
-
 # Examples
 ```jldoctest
 julia> R, (y, z) = polynomial_ring(QQ, [:y, :z])
@@ -1007,10 +999,9 @@ end
     f_kk = map_coefficients(kk, f)
     h = first([h for (h, _) in factor(f_kk)])
     kk_ext, zeta = extension_field(h)
-    iso_kk_ext = hom(L, kk_ext, zeta)
-    br = base_ring(P_prime)
+    inc_kk_ext = hom(L, kk_ext, zeta)
     LoR, to_LoR = change_base_ring(kk_ext, R)
-    help_map = hom(br, LoR, iso_kk_ext, to_LoR.(iso.(gens(R_exp))))
+    help_map = hom(RR, LoR, inc_kk_ext, to_LoR.(iso.(gens(R_exp))))
     P_prime_ext = ideal(LoR, help_map.(gens(P_prime)))
     push!(full_res, (P, Q, P_prime_ext, degree(h)))
   end
@@ -1178,7 +1169,6 @@ end
 function minimal_primes(
     I::MPolyIdeal{T}; 
     algorithm::Symbol=:GTZ, 
-    factor_generators::Bool=true,
     simplify_ring::Bool=true,
     cache::Bool=true
   ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
@@ -1193,57 +1183,14 @@ function minimal_primes(
     W, id, id_inv = simplify(Q)
     @assert domain(id) === codomain(id_inv) === Q
     @assert codomain(id) === domain(id_inv) === W
-    res_simp = minimal_primes(modulus(W); algorithm, factor_generators, simplify_ring=false)
+    res_simp = minimal_primes(modulus(W); algorithm, simplify_ring=false)
     result = [I + ideal(R, lift.(id_inv.(W.(gens(j))))) for j in res_simp]
     for p in result
       set_attribute!(p, :is_prime=>true)
     end
     return result
   end
 
-  # This will in many cases lead to an easy simplification of the problem
-  if factor_generators
-    J = [ideal(R, gens(I))] # A copy of I as initialization
-    for g in gens(I)
-      K = typeof(I)[]
-      is_zero(g) && continue
-      for (b, k) in factor(g)
-        # Split the already collected components with b
-        for j in J
-          push!(K, j + ideal(R, b))
-        end
-      end
-      J = K
-    end
-
-    unique_comp = typeof(I)[]
-    for q in J
-      is_one(q) && continue
-      q in unique_comp && continue
-      push!(unique_comp, q)
-    end
-    J = unique_comp
-
-    # unique! seems to fail here. We have to do it manually.
-    pre_result = filter!(!is_one, vcat([minimal_primes(j; algorithm, factor_generators=false) for j in J]...))
-    result = typeof(I)[]
-    for p in pre_result
-      p in result && continue
-      push!(result, p)
-    end
-
-    # The list might not consist of minimal primes only. We have to discard the embedded ones
-    final_list = typeof(I)[]
-    for p in result
-      any((q !== p && is_subset(q, p)) for q in result) && continue
-      push!(final_list, p)
-    end
-    for p in final_list
-      set_attribute!(p, :is_prime=>true)
-    end
-    return final_list
-  end
-
   R_flat, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
   I_flat = ideal(R_flat, elem_type(R_flat)[g for g in iso_inv.(gens(I))])
   dec = minimal_primes(I_flat; algorithm)

src/Rings/primary_decomposition_helpers.jl

@@ -169,13 +169,22 @@ end
   J = saturated_ideal(I)
   res = absolute_primary_decomposition(J)
   result = []
+  if is_empty(res)
+    U = ideal(A, one(A))
+    return [(U, U, U, 0)]
+  end
+
+  # Create the ring for the return values
   for (P, Q, P_prime, d) in res
+    R_prime = base_ring(P_prime) # the new polynomial ring
+    L = coefficient_ring(R_prime) # the new field for the result
+    A_ext, ext_map = change_base_ring(L, A) # recreate the quo-ring over that field
+    @assert coefficient_ring(base_ring(A_ext)) === L
+    help_map = hom(R_prime, A_ext, gens(A_ext); check=false)
     PP = ideal(A, A.(gens(P)))
     QQ = ideal(A, A.(gens(Q)))
-    R_prime = base_ring(P_prime)
-    L = coefficient_ring(R_prime)
-    A_ext, ext_map = change_base_ring(L, R_prime)
-    PP_prime = ideal(A_ext, ext_map.(gens(P_prime)))
+    trans_gens = help_map.(gens(P_prime))
+    PP_prime = ideal(A_ext, trans_gens)
     push!(result, (PP, QQ, PP_prime, d))
   end
   return result