Primary decomposition rerouting

src/Rings/MPolyMap/MPolyAnyMap.jl

@@ -250,7 +250,14 @@
 # Julia functions in both maps
 function compose(F::MPolyAnyMap{D, C, <: Function}, G::MPolyAnyMap{C, E, <: Function}) where {D, C, E}
   @req codomain(F) === domain(G) "Incompatible (co)domain in composition"
-  return hom(domain(F), codomain(G), x -> coefficient_map(G)(coefficient_map(F)(x)), G.(_images(F)), check=false)
+  b = coefficient_map(F)(one(coefficient_ring(domain(F))))
+  if parent(b) === domain(G)
+    return hom(domain(F), codomain(G), x -> G(coefficient_map(F)(x)), G.(_images(F)), check=false)
+  elseif parent(b) === coefficient_ring(domain(G))
+    return hom(domain(F), codomain(G), x -> coefficient_map(G)(coefficient_map(F)(x)), G.(_images(F)), check=false)
+  else
+    error("coefficient map is not admissible")
+  end
 end
 
 # Now compose with arbitrary maps

src/Rings/Rings.jl

@@ -35,4 +35,4 @@ include("PBWAlgebra.jl")
 include("PBWAlgebraQuo.jl")
 include("FreeAssAlgIdeal.jl")
 include("hilbert.jl")
-
+include("primary_decomposition_helpers.jl")

src/Rings/mpoly-ideals.jl

@@ -531,6 +531,23 @@ function primary_decomposition(I::T; algorithm::Symbol=:GTZ, cache::Bool=true) w
   end::Vector{Tuple{T,T}}
 end
 
+function primary_decomposition(
+    I::MPolyIdeal{T}; 
+    algorithm::Symbol=:GTZ, cache::Bool=true
+  ) where {U<:Union{nf_elem, <:Hecke.NfRelElem}, T<:MPolyRingElem{U}}
+  return get_attribute!(I, :primary_decomposition) do
+    R = base_ring(I)
+    R_flat, iso, iso_inv = _flatten_to_QQ(R)
+    I_flat = ideal(R_flat, iso_inv.(gens(I)))
+    dec = primary_decomposition(I_flat; algorithm)
+    result = Vector{Tuple{typeof(I), typeof(I)}}()
+    for (P, Q) in dec
+      push!(result, (ideal(R, iso.(gens(P))), ideal(R, iso.(gens(Q)))))
+    end
+    result
+  end::Vector{Tuple{typeof(I), typeof(I)}}
+end
+
 function _compute_primary_decomposition(I::MPolyIdeal; algorithm::Symbol=:GTZ)
   R = base_ring(I)
   if isa(base_ring(R), NumField) && !isa(base_ring(R), AnticNumberField)

src/Rings/mpolyquo-localizations.jl

@@ -1884,14 +1884,17 @@ end
 ############################################################################# 
 
 @doc raw"""
-     primary_decomposition(I::Union{<:MPolyQuoIdeal, <:MPolyQuoLocalizedIdeal, <:MPolyLocalizedIdeal})
+    primary_decomposition(I::Union{<:MPolyQuoIdeal, <:MPolyQuoLocalizedIdeal, <:MPolyLocalizedIdeal})
 
 Return the primary decomposition of ``I``.
 """
-function primary_decomposition(I::Union{<:MPolyQuoIdeal, <:MPolyQuoLocalizedIdeal, <:MPolyLocalizedIdeal})
+function primary_decomposition(
+    I::Union{<:MPolyQuoIdeal, <:MPolyQuoLocalizedIdeal, <:MPolyLocalizedIdeal};
+    algorithm::Symbol=:GTZ, cache::Bool=true
+  )
   Q = base_ring(I)
   R = base_ring(Q)
-  decomp = primary_decomposition(saturated_ideal(I))
+  decomp = primary_decomposition(saturated_ideal(I); algorithm, cache)
   result = [(ideal(Q, Q.(gens(a))), ideal(Q, Q.(gens(b)))) for (a, b) in decomp]
   erase = Int[]
   for i in 1:length(result)

src/Rings/primary_decomposition_helpers.jl

@@ -0,0 +1,60 @@
+# Since these need the declaration of MPolyQuoRing, they have to come later here 
+# in an extra file.
+function _flatten_ring(R::MPolyRing{T}) where {T<:QQFieldElem}
+  return R, identity_map(R), identity_map(R)
+end
+
+function _flatten_ring(Q::MPolyQuoRing{<:MPolyRingElem{T}}) where {T<:QQFieldElem}
+  return Q, identity_map(Q), identity_map(Q)
+end
+
+function _flatten_ring(R::MPolyRing{T}) where {T<:Union{nf_elem, <:Hecke.NfRelElem}}
+  K = coefficient_ring(R)
+  alpha = first(gens(K))
+  kk = base_field(K)
+  P, _ = polynomial_ring(kk, vcat([:theta], symbols(R)), cached=false)
+  theta = first(gens(P))
+  f = defining_polynomial(K)
+  d = degree(f)
+  powers = elem_type(P)[theta^k for k in 0:d-1]
+  function help_map(a)
+    return sum(c*powers[i] for (i, c) in enumerate(coefficients(a)); init=zero(P))
+  end
+  R_flat, pr = quo(P, ideal(P, evaluate(f, theta)))
+  to_R_flat = hom(R, R_flat, x->pr(help_map(x)), gens(R_flat)[2:end])
+  to_R = hom(R_flat, R, vcat([R(alpha)], gens(R)))
+  return R_flat, to_R, to_R_flat
+end
+
+function _flatten_ring(Q::MPolyQuoRing{<:MPolyRingElem{T}}) where {T<:Union{nf_elem, <:Hecke.NfRelElem}}
+  R = base_ring(Q)
+  R_flat, iso, iso_inv = _flatten_ring(R)
+  I = modulus(Q)
+  I_flat = ideal(R_flat, iso_inv.(gens(I)))
+  Q_flat, pr = quo(R_flat, I_flat)
+  alpha = first(gens(coefficient_ring(Q)))
+  to_Q = hom(Q_flat, Q, vcat([Q(alpha)], gens(Q)))
+  theta = Q_flat[1]
+  f = defining_polynomial(coefficient_ring(Q))
+  d = degree(f)
+  powers = [theta^k for k in 0:d-1]
+  function help_map(a)
+    return sum(c*powers[i] for (i, c) in enumerate(coefficients(a)); init=zero(Q_flat))
+  end
+  to_Q_flat = hom(Q, Q_flat, help_map, gens(Q_flat)[2:end])
+  return Q_flat, to_Q, to_Q_flat
+end
+
+function _flatten_to_QQ(R::Union{<:MPolyRing, <:MPolyQuoRing})
+  R_flat, to_R, to_R_flat = _flatten_ring(R)
+  res, a, b = _flatten_to_QQ(R_flat)
+  return res, compose(a, to_R), compose(to_R_flat, b)
+end
+
+function _flatten_to_QQ(R::MPolyRing{T}) where {T<:QQFieldElem}
+  return _flatten_ring(R)
+end
+
+function _flatten_to_QQ(R::MPolyQuoRing{<:MPolyRingElem{T}}) where {T<:QQFieldElem}
+  return _flatten_ring(R)
+end

test/Rings/mpoly.jl

@@ -331,3 +331,40 @@ end
   @test hessian(f) == -24*x*z + 24*y^2*z
 end
 
+@testset "rerouting of primary decomposition over number fields" begin
+  P, t = QQ[:t]
+  kk = splitting_field(t^3 - 1)
+  kks, s = kk[:s]
+  kk2 = splitting_field(s^2 + 1)
+
+  R, (x, y) = kk[:x, :y]
+  alpha = first(gens(kk))
+  R_flat, iso, iso_inv = Oscar._flatten_ring(R)
+  theta = first(gens(R_flat))
+  @test iso_inv(R(alpha)) == theta
+  @test iso(theta) == alpha
+  id1 = compose(iso, iso_inv)
+  @test all(x->x==id1(x), gens(R_flat))
+  id2 = compose(iso_inv, iso)
+  @test all(x->x==id2(x), gens(R))
+
+  S, (u, v) = kk2[:u, :v]
+  S_flat, iso, iso_inv = Oscar._flatten_ring(S)
+  S_ff, iso2, iso_inv2 = Oscar._flatten_ring(S_flat)
+
+  iso_tot = compose(iso2, iso)
+  iso_inv_tot = compose(iso_inv, iso_inv2)
+
+  beta = first(gens(kk2))
+  @test iso_inv_tot(S(beta)) == S_ff[2]
+  @test iso_inv_tot(S(alpha)) == S_ff[1]
+
+  Sff, iso2, iso_inv2 = Oscar._flatten_to_QQ(S)
+  @test iso_inv_tot(S(beta)) == S_ff[2]
+  @test iso_inv_tot(S(alpha)) == S_ff[1]
+
+  I = ideal(S, [(u^2 + 1)^3])
+  dec = primary_decomposition(I)
+  @test length(dec) == 2
+end
+