Overhaul preprocessing for radical computations.

experimental/GaloisGrp/src/Solve.jl

@@ -236,14 +236,14 @@
   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

@@ -613,7 +613,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,10 +554,15 @@ end
   return is_prime(saturated_ideal(I))
 end
 
-@attr typeof(I) function radical(I::MPolyQuoIdeal; eliminate_variables::Bool=true)
+@attr typeof(I) function radical(                                                   
+    I::MPolyQuoIdeal; 
+    preprocessing::Bool=true,
+    eliminate_variables::Bool=true, 
+    factor_generators::Bool=true
+  )
   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; preprocessing, eliminate_variables, factor_generators))) if !iszero(g)])
 end
 
 # The following is to streamline the programmer's

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

@@ -490,7 +490,7 @@ end
 
 #######################################################
 @doc raw"""
-    radical(I::MPolyIdeal{T}; eliminate_variables::Bool = true) where {T <: MPolyRingElem}
+    radical(I::MPolyIdeal{T}; preprocessing::Bool=true) where {T <: MPolyRingElem}
 
 Return the radical of `I`.
 
@@ -500,7 +500,7 @@ Return the radical of `I`.
   * If the `base_ring` of `I` is a number field, then we first expand the minimal polynomials to reduce to a computation over the rationals.
   * For polynomial rings over the integers, the algorithm proceeds as suggested by Pfister, Sadiq, and Steidel. See [KL91](@cite), [Kem02](@cite), and [PSS11](@cite).
 
-When `eliminate_variables` is set to `true`, a preprocessing heuristic is applied in order to find variables which can be eliminated using `I`.
+When `preprocessing` is set to `true`, several preprocessing heuristics are applied. See the source code for more detailed information and options.
 
 # Examples
 ```jldoctest
@@ -590,48 +590,44 @@ 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
 ```
 """
-@attr MPolyIdeal{T} function radical(I::MPolyIdeal{T}; eliminate_variables::Bool=true) where {T <: MPolyRingElem}
-  if eliminate_variables
+@attr MPolyIdeal{T} function radical(
+    I::MPolyIdeal{T}; 
+    preprocessing::Bool=true,
+    eliminate_variables::Bool=true,
+    factor_generators::Bool=true
+  ) where {T <: MPolyRingElem}
+  # We first check for the squarefree flag. 
+  # Elimination of variables will probably affect how the generators factor. 
+  if preprocessing && factor_generators
+    # Replace every generator by its squarefree part to get closer to the actual radical 
+    # even before passing computations on to Singular.
+    return radical(_squarefree_generator_ideal(I); preprocessing, eliminate_variables, factor_generators=false)
+  end
+  if preprocessing && eliminate_variables
     is_known_to_be_radical(I) && return I
-    # Calling `elimpart` (within `simplify`) turns out to significantly speed things up in many cases.
+    # Call `elimpart` (within `simplify`) in order to substitute variables. 
+    # This turns out to significantly speed things up in many cases.
     R = base_ring(I)
     Q, pr = quo(R, I)
     S, iso, iso_inv = simplify(Q)
     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; preprocessing, eliminate_variables=false, factor_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)
@@ -660,38 +656,40 @@ 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
+    fact = factor_squarefree(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
 @attr MPolyIdeal{T} function radical(
     I::MPolyIdeal{T};
+    preprocessing::Bool=true,
     factor_generators::Bool=true, 
     eliminate_variables::Bool=true
-  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
+  ) where {U<:Hecke.RelSimpleNumFieldElem, T<:MPolyRingElem{U}}
   is_known_to_be_radical(I) && return I
   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
+  preprocessing && factor_generators && return radical(_squarefree_generator_ideal(I); preprocessing, eliminate_variables, factor_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; 
+                       preprocessing,
+                       eliminate_variables, 
+                       factor_generators=false) # Taking the squarefree part of generators only pays off on the top level. 
   Irad = iso(I_flat_rad)
   set_attribute!(Irad, :is_radical => true)
-  @hassert :IdealSheaves 2 !is_one(Irad)
   return Irad
 end
 
@@ -830,7 +828,7 @@ end
 function primary_decomposition(
     I::MPolyIdeal{T}; 
     algorithm::Symbol=:GTZ, cache::Bool=true
-  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
+  ) where {U<:Hecke.RelSimpleNumFieldElem, T<:MPolyRingElem{U}}
   if has_attribute(I, :primary_decomposition)
     return get_attribute(I, :primary_decomposition)::Vector{Tuple{typeof(I), typeof(I)}}
   end
@@ -899,9 +897,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])
@@ -1015,10 +1010,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
@@ -1186,8 +1180,7 @@ end
 function minimal_primes(
     I::MPolyIdeal{T}; 
     algorithm::Symbol=:GTZ, 
-    factor_generators::Bool=true,
-    simplify_ring::Bool=true,
+    eliminate_variables::Bool=true,
     cache::Bool=true
   ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
   has_attribute(I, :minimal_primes) && return get_attribute(I, :minimal_primes)::Vector{typeof(I)}
@@ -1196,62 +1189,19 @@ function minimal_primes(
   is_one(I) && return typeof(I)[]
 
   # Try to eliminate variables first. This will often speed up computations significantly.
-  if simplify_ring
+  if eliminate_variables
     Q, pr = quo(R, I)
     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, eliminate_variables=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)