mpoly-ideals.jl

# constructors #######################################################

@doc raw"""
    ideal(R::MPolyRing, V::Vector)

Given a vector `V` of polynomials in `R`, return the ideal of `R` generated by these polynomials.

!!! note
    In the graded case, the entries of `V` must be homogeneous.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x*y-3*x,y^3-2*x^2*y])
Ideal generated by
  x*y - 3*x
  -2*x^2*y + y^3

julia> typeof(I)
MPolyIdeal{QQMPolyRingElem}

julia> S, (x, y) = graded_polynomial_ring(QQ, ["x", "y"],  [1, 2])
(Graded multivariate polynomial ring in 2 variables over QQ, MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[x, y])

julia> J = ideal(S, [(x^2+y)^2])
Ideal generated by
  x^4 + 2*x^2*y + y^2

julia> typeof(J)
MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}
```
"""
function ideal(R::MPolyRing, g::Vector)
  h = elem_type(R)[R(f) for f = g]
  #if isempty(g)
  #  push!(h, R())
  #end
  return MPolyIdeal(R, h)
end

function ideal(I::IdealGens{T}) where {T <: MPolyRingElem}
  return MPolyIdeal(I)
end

function ideal(x::MPolyRingElem{T}) where T <: RingElem
  return ideal([x])
end

function ideal(Qxy::MPolyRing{T}, x::MPolyRingElem{T}) where T <: RingElem
  return ideal(Qxy, [x])
end

function ideal(g::Vector{T}) where {T <: MPolyRingElem}
  @req length(g) > 0 "List of elements must be non-empty"
  return ideal(parent(g[1]), g)
end

# Coerce an ungraded ideal in a graded ring
function ideal(S::MPolyDecRing, I::MPolyIdeal)
  @req base_ring(I) === forget_grading(S) "Rings do not coincide"
  return ideal(S, [ S(f) for f in gens(I) ])
end

function is_graded(I::MPolyIdeal)
  return is_graded(Hecke.ring(I))
end

# elementary operations #######################################################
@doc raw"""
    check_base_rings(I::MPolyIdeal, J::MPolyIdeal)

Throws an error if the base rings of the ideals `I` and `J` do not coincide.
"""
function check_base_rings(I::MPolyIdeal, J::MPolyIdeal)
  if !isequal(base_ring(I), base_ring(J))
    error("Base rings must coincide.")
  end
end

@doc raw"""
    ^(I::MPolyIdeal, m::Int)

Return the `m`-th power of `I`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [x, y])
Ideal generated by
  x
  y

julia> I^3
Ideal generated by
  x^3
  x^2*y
  x*y^2
  y^3
```
"""
function Base.:^(I::MPolyIdeal, m::Int)
  return MPolyIdeal(base_ring(I), (singular_generators(I))^m)
end

@doc raw"""
    +(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T

Return the sum of `I` and `J`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [x, y])
Ideal generated by
  x
  y

julia> J = ideal(R, [z^2])
Ideal generated by
  z^2

julia> I+J
Ideal generated by
  x
  y
  z^2
```
"""
function Base.:+(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T
  oscar_assure(I)
  oscar_assure(J)
  check_base_rings(I, J)
  newgens = filter!(!iszero, unique!(vcat(I.gens.O,J.gens.O)))
  return ideal(base_ring(I), newgens)
end
Base.:-(I::MPolyIdeal, J::MPolyIdeal) = I+J

@doc raw"""
    *(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T

Return the product of `I` and `J`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [x, y])
Ideal generated by
  x
  y

julia> J = ideal(R, [z^2])
Ideal generated by
  z^2

julia> I*J
Ideal generated by
  x*z^2
  y*z^2
```
"""
function Base.:*(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T
  oscar_assure(I)
  oscar_assure(J)
  check_base_rings(I, J)
  newgens = elem_type(base_ring(I))[]
  for g in I.gens.O
    for h in J.gens.O
      gh = g * h
      if !iszero(gh)
        push!(newgens, gh)
      end
    end
  end
  return ideal(base_ring(I), unique!(newgens))
end

#######################################################

# ideal intersection #######################################################
@doc raw"""
    intersect(I::MPolyIdeal{T}, Js::MPolyIdeal{T}...) where T
    intersect(V::Vector{MPolyIdeal{T}}) where T

Return the intersection of two or more ideals.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2;

julia> J = ideal(R, [y^2-x^3+x]);

julia> intersect(I, J)
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> intersect([I, J])
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2
```
"""
function Base.intersect(I::MPolyIdeal{T}, Js::MPolyIdeal{T}...) where T
  si = singular_generators(I)
  si = Singular.intersection(si, [singular_generators(J) for J in Js]...)
  return MPolyIdeal(base_ring(I), si)
end

function Base.intersect(V::Vector{MPolyIdeal{T}}) where T
  @assert length(V) != 0
  length(V) == 1 && return V[1]

  return Base.intersect(V[1], V[2:end]...)
end


#######################################################

@doc raw"""
    quotient(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T

Return the ideal quotient of `I` by `J`. Alternatively, use `I:J`.

    quotient(I::MPolyIdeal{T}, f::MPolyRingElem{T}) where T

Return the ideal quotient of `I` by the ideal generated by `f`. Alternatively, use `I:f`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [x^4+x^2*y*z+y^3*z, y^4+x^3*z+x*y^2*z, x^3*y+x*y^3])
Ideal generated by
  x^4 + x^2*y*z + y^3*z
  x^3*z + x*y^2*z + y^4
  x^3*y + x*y^3

julia> J = ideal(R, [x, y, z])^2
Ideal generated by
  x^2
  x*y
  x*z
  y^2
  y*z
  z^2

julia> L = quotient(I, J)
Ideal generated by
  x^3*z + x*y^2*z + y^4
  x^3*y + x*y^3
  x^4 + x^2*y*z + y^3*z
  x^3*z^2 - x^2*y*z^2 + x*y^2*z^2 - y^3*z^2
  x^2*y^2*z - x^2*y*z^2 - y^3*z^2
  x^3*z^2 + x^2*y^3 - x^2*y^2*z + x*y^2*z^2

julia> I:J
Ideal generated by
  x^3*z + x*y^2*z + y^4
  x^3*y + x*y^3
  x^4 + x^2*y*z + y^3*z
  x^3*z^2 - x^2*y*z^2 + x*y^2*z^2 - y^3*z^2
  x^2*y^2*z - x^2*y*z^2 - y^3*z^2
  x^3*z^2 + x^2*y^3 - x^2*y^2*z + x*y^2*z^2

julia> I:x
Ideal generated by
  x^2*y + y^3
  x^3*z + x*y^2*z + y^4
  x^2*z^2 + x*y^3 - x*y^2*z + y^2*z^2
  x^4
  x^3*z^2 - x^2*z^3 + 2*x*y^2*z^2 - y^2*z^3
  -x^2*z^4 + x*y^2*z^3 - y^2*z^4
```
"""
function quotient(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T
  @assert base_ring(I) == base_ring(J)
  return MPolyIdeal(base_ring(I), Singular.quotient(singular_generators(I), singular_generators(J)))
end

function quotient(I::MPolyIdeal{T}, f::T) where T
  R = base_ring(I)
  @assert R == parent(f)
  return quotient(I, ideal(R, [f]))
end
(::Colon)(I::MPolyIdeal, J::MPolyIdeal) = quotient(I, J)
(::Colon)(I::MPolyIdeal, f::MPolyRingElem) = quotient(I, f)

#######################################################

# saturation #######################################################
@doc raw"""
    saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T

Return the saturation of `I` with respect to `J`.
If the second ideal `J` is not given, the ideal generated by the generators (variables) of `base_ring(I)` is used.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [z^3, y*z^2, x*z^2, y^2*z, x*y*z, x^2*z, x*y^2, x^2*y])
Ideal generated by
  z^3
  y*z^2
  x*z^2
  y^2*z
  x*y*z
  x^2*z
  x*y^2
  x^2*y

julia> J = ideal(R, [x, y, z])
Ideal generated by
  x
  y
  z

julia> K = saturation(I, J)
Ideal generated by
  z
  x*y

julia> K = saturation(I)
Ideal generated by
  z
  x*y
```
"""
function saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T
  K, _ = Singular.saturation(singular_generators(I), singular_generators(J))
  return MPolyIdeal(base_ring(I), K)
end

# the following is corresponding to saturation2 from Singular
# TODO: think about how to use use this properly/automatically
function _saturation2(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T
  K, _ = Singular.saturation2(singular_generators(I), singular_generators(J))
  return MPolyIdeal(base_ring(I), K)
end

#######################################################
@doc raw"""
    saturation_with_index(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T

Return $I:J^{\infty}$ together with the smallest integer $m$ such that $I:J^m = I:J^{\infty}$.
If the second ideal `J` is not given, the ideal generated by the generators (variables) of `base_ring(I)` is used.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [z^3, y*z^2, x*z^2, y^2*z, x*y*z, x^2*z, x*y^2, x^2*y])
Ideal generated by
  z^3
  y*z^2
  x*z^2
  y^2*z
  x*y*z
  x^2*z
  x*y^2
  x^2*y

julia> J = ideal(R, [x, y, z])
Ideal generated by
  x
  y
  z

julia> K, m = saturation_with_index(I, J)
(Ideal (z, x*y), 2)

julia> K, m = saturation_with_index(I)
(Ideal (z, x*y), 2)
```
"""
function saturation_with_index(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T
  K, k = Singular.saturation(singular_generators(I), singular_generators(J))
  return (MPolyIdeal(base_ring(I), K), k)
 end

# elimination #######################################################
@doc raw"""
    eliminate(I::MPolyIdeal{T}, V::Vector{T}) where T <: MPolyRingElem

Given a vector `V` of polynomials which are variables, these variables are eliminated from `I`.
That is, return the ideal generated by all polynomials in `I` which only involve the remaining variables.

    eliminate(I::MPolyIdeal, V::AbstractVector{Int})

Given a vector `V` of indices which specify variables, these variables are eliminated from `I`.
That is, return the ideal generated by all polynomials in `I` which only involve the remaining variables.

!!! note
    The return value is an ideal of the original ring.

# Examples
```jldoctest
julia> R, (t, x, y, z) = polynomial_ring(QQ, ["t", "x", "y", "z"])
(Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[t, x, y, z])

julia> I = ideal(R, [t-x, t^2-y, t^3-z])
Ideal generated by
  t - x
  t^2 - y
  t^3 - z

julia> A = [t]
1-element Vector{QQMPolyRingElem}:
 t

julia> TC = eliminate(I, A)
Ideal generated by
  -x*z + y^2
  x*y - z
  x^2 - y

julia> A = [1]
1-element Vector{Int64}:
 1

julia> TC = eliminate(I, A)
Ideal generated by
  -x*z + y^2
  x*y - z
  x^2 - y

julia> base_ring(TC)
Multivariate polynomial ring in 4 variables t, x, y, z
  over rational field
```
"""
function eliminate(I::MPolyIdeal{T}, l::Vector{T}) where T <: MPolyRingElem
  S = singular_polynomial_ring(I)
  s = Singular.eliminate(singular_generators(I), [S(x) for x = l]...)
  return MPolyIdeal(base_ring(I), s)
end
function eliminate(I::MPolyIdeal, l::AbstractVector{Int})
  R = base_ring(I)
  return eliminate(I, [gen(R, i) for i=l])
end

### todo: wenn schon GB bzgl. richtiger eliminationsordnung bekannt ...
### Frage: return MPolyIdeal(base_ring(I), s) ???

###################################################

# primary decomposition #######################################################

#######################################################
@doc raw"""
    radical(I::MPolyIdeal)

Return the radical of `I`.

# Implemented Algorithms

If the base ring of `I` is a polynomial
ring over a field, a combination of the algorithms of Krick and Logar
(with modifications by Laplagne) and Kemper is used. For polynomial
rings over the integers, the algorithm proceeds as suggested by
Pfister, Sadiq, and Steidel. See [KL91](@cite),
[Kem02](@cite), and [PSS11](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> RI = radical(I)
Ideal generated by
  x^4 - x^3*y - x^3 - x^2 - x*y^2 + x*y + x + y^3 + y^2
```
```jldoctest
julia> R, (a, b, c, d) = polynomial_ring(ZZ, ["a", "b", "c", "d"])
(Multivariate polynomial ring in 4 variables over ZZ, ZZMPolyRingElem[a, b, c, d])

julia> I = intersect(ideal(R, [9,a,b]), ideal(R, [3,c]))
Ideal generated by
  9
  3*b
  3*a
  b*c
  a*c

julia> I = intersect(I, ideal(R, [11,2a,7b]))
Ideal generated by
  99
  3*b
  3*a
  b*c
  a*c

julia> I = intersect(I, ideal(R, [13a^2,17b^4]))
Ideal generated by
  39*a^2
  13*a^2*c
  51*b^4
  17*b^4*c
  3*a^2*b^4
  a^2*b^4*c

julia> I = intersect(I, ideal(R, [9c^5,6d^5]))
Ideal generated by
  78*a^2*d^5
  117*a^2*c^5
  102*b^4*d^5
  153*b^4*c^5
  6*a^2*b^4*d^5
  9*a^2*b^4*c^5
  39*a^2*c^5*d^5
  51*b^4*c^5*d^5
  3*a^2*b^4*c^5*d^5

julia> I = intersect(I, ideal(R, [17,a^15,b^15,c^15,d^15]))
Ideal generated by
  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

julia> RI = radical(I)
Ideal generated by
  102*b*d
  78*a*d
  51*b*c
  39*a*c
  6*a*b*d
  3*a*b*c
```
"""
@attr T function radical(I::T) where {T <: MPolyIdeal}
  R = base_ring(I)
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    r = radical(map_coefficients(pseudo_inv(mA), I))
    Irad = map_coefficients(mA, r, parent = R)
    set_attribute!(Irad, :is_radical => true)
    return Irad
  elseif elem_type(base_ring(R)) <: FieldElement
    J = Singular.LibPrimdec.radical(singular_polynomial_ring(I), singular_generators(I))
  elseif base_ring(singular_polynomial_ring(I)) isa Singular.Integers
    J = Singular.LibPrimdecint.radicalZ(singular_polynomial_ring(I), singular_generators(I))
  else
    error("not implemented for base ring")
  end
  Irad = ideal(R, J)
  set_attribute!(Irad, :is_radical => true)
  return Irad
end

# Rerouting via expansion of the coefficient field
function radical(
    I::MPolyIdeal{T};
    factor_generators::Bool=true
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  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
    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)
    Irad = iso(I_flat_rad)
    set_attribute!(Irad, :is_radical => true)
    @hassert :IdealSheaves 2 !is_one(Irad)
    Irad
  end::MPolyIdeal{T}
end

function map_coefficients(mp, I::MPolyIdeal; parent = nothing)
  if parent === nothing
    parent = Oscar.parent(map_coefficients(mp, gen(I, 1)))
  end
  return ideal(parent, [map_coefficients(mp, g, parent = parent) for g = gens(I)])
end

@doc raw"""
    is_radical(I::MPolyIdeal)

Return whether `I` is a radical ideal.

Computes the radical.
"""
@attr Bool function is_radical(I::MPolyIdeal)
  if has_attribute(I, :is_prime) && is_prime(I)
    return true
  end

  return I == radical(I)
end
#######################################################
@doc raw"""
    primary_decomposition(I::MPolyIdeal; algorithm = :GTZ, cache=true)

Return a minimal primary decomposition of `I`.

The decomposition is returned as a vector of tuples $(Q_1, P_1), \dots, (Q_t, P_t)$, say,
where each $Q_i$ is a primary ideal with associated prime $P_i$, and where the intersection of
the $Q_i$ is `I`.

# Implemented Algorithms

If the base ring of `I` is a polynomial ring over a field, the algorithm of Gianni, Trager, and Zacharias
is used by default (`algorithm = :GTZ`). Alternatively, the algorithm by Shimoyama and Yokoyama can be used
by specifying `algorithm = :SY`.  For polynomial rings over the integers, the algorithm proceeds as suggested by
Pfister, Sadiq, and Steidel. See [GTZ88](@cite), [SY96](@cite), and [PSS11](@cite).

!!! warning
    The algorithm of Gianni, Trager, and Zacharias may not terminate over a small finite field. If it terminates, the result is correct.

!!! warning
    If computations are done in a ring over a number field, then the output may contain redundant components.

If `cache=false` is set, the primary decomposition is recomputed and not cached.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = primary_decomposition(I)
3-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}}:
 (Ideal (x^3 - x - y^2), Ideal (x^3 - x - y^2))
 (Ideal (x^2 - 2*x*y - 2*x + y^2 + 2*y + 1), Ideal (x - y - 1))
 (Ideal (y, x^2), Ideal (x, y))

julia> L = primary_decomposition(I, algorithm = :SY, cache=false)
3-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}}:
 (Ideal (x^3 - x - y^2), Ideal (x^3 - x - y^2))
 (Ideal (x^2 - 2*x*y - 2*x + y^2 + 2*y + 1), Ideal (x - y - 1))
 (Ideal (y, x^2), Ideal (y, x))
```
```jldoctest
julia> R, (a, b, c, d) = polynomial_ring(ZZ, ["a", "b", "c", "d"])
(Multivariate polynomial ring in 4 variables over ZZ, ZZMPolyRingElem[a, b, c, d])

julia> I = ideal(R, [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])
Ideal generated by
  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

julia> L = primary_decomposition(I)
8-element Vector{Tuple{MPolyIdeal{ZZMPolyRingElem}, MPolyIdeal{ZZMPolyRingElem}}}:
 (Ideal (d^5, c^5), Ideal (d, c))
 (Ideal (a^2, b^4), Ideal (b, a))
 (Ideal (2, c^5), Ideal (2, c))
 (Ideal (3), Ideal (3))
 (Ideal (13, b^4), Ideal (13, b))
 (Ideal (17, a^2), Ideal (17, a))
 (Ideal (17, d^15, c^15, b^15, a^15), Ideal (17, d, c, b, a))
 (Ideal (9, 3*d^5, d^10), Ideal (3, d))
```
"""
function primary_decomposition(I::T; algorithm::Symbol=:GTZ, cache::Bool=true) where {T<:MPolyIdeal}
  !cache && return _compute_primary_decomposition(I, algorithm=algorithm)
  return get_attribute!(I, :primary_decomposition) do
    return _compute_primary_decomposition(I, algorithm=algorithm)
  end::Vector{Tuple{T,T}}
end

function primary_decomposition(
    I::MPolyIdeal{T}; 
    algorithm::Symbol=:GTZ, cache::Bool=true
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  if has_attribute(I, :primary_decomposition)
    return get_attribute(I, :primary_decomposition)::Vector{Tuple{typeof(I), typeof(I)}}
  end
  R = base_ring(I)
  R_flat, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  I_flat = ideal(R_flat, iso_inv.(gens(I)))
  dec = primary_decomposition(I_flat; algorithm, cache)
  result = Vector{Tuple{typeof(I), typeof(I)}}()
  for (P, Q) in dec
    push!(result, (ideal(R, unique!([x for x in iso.(gens(P)) if !iszero(x)])), 
                   ideal(R, unique!([x for x in iso.(gens(Q)) if !iszero(x)]))))
  end

  for (Q,P) in result
    set_attribute!(P, :is_prime=>true)
    set_attribute!(Q, :is_primary=>true)
  end

  cache && set_attribute!(I, :primary_decomposition=>result)
  return result
end

function _compute_primary_decomposition(I::MPolyIdeal; algorithm::Symbol=:GTZ)
  R = base_ring(I)
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    pd = primary_decomposition(map_coefficients(pseudo_inv(mA), I), algorithm = algorithm, cache = false)
    if isempty(pd)
      return Tuple{typeof(I), typeof(I)}[]
    end
    return Tuple{typeof(I), typeof(I)}[(map_coefficients(mA, x[1], parent = R), map_coefficients(mA, x[2], parent = R)) for x = pd]
  end
  if elem_type(base_ring(R)) <: FieldElement
    if algorithm == :GTZ
      L = Singular.LibPrimdec.primdecGTZ(singular_polynomial_ring(I), singular_generators(I))
    elseif algorithm == :SY
      L = Singular.LibPrimdec.primdecSY(singular_polynomial_ring(I), singular_generators(I))
    else
      error("algorithm invalid")
    end
  elseif base_ring(singular_polynomial_ring(I)) isa Singular.Integers
    L = Singular.LibPrimdecint.primdecZ(singular_polynomial_ring(I), singular_generators(I))
  else
    error("base ring not implemented")
  end
  V = [(ideal(R, q[1]), ideal(R, q[2])) for q in L]
  if length(V) == 1 && is_one(gen(V[1][1], 1))
    return Tuple{typeof(I), typeof(I)}[]
  end
  return V
end

########################################################
@doc raw"""
    absolute_primary_decomposition(I::MPolyIdeal{<:MPolyRingElem{QQFieldElem}})

Given an ideal `I` in a multivariate polynomial ring over the rationals, return an absolute minimal primary decomposition of `I`.

Return the decomposition as a vector of tuples $(Q_i, P_i, P_{ij}, d_{ij})$, say,
where $(Q_i, P_i)$ is a (primary, prime) tuple as returned by `primary_decomposition(I)`,
and $P_{ij}$ represents a corresponding class of conjugated absolute associated primes
defined over a number field of degree $d_{ij}$ whose generator prints as `_a`.

# Implemented Algorithms

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"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[y, z])

julia> p = z^2+1
z^2 + 1

julia> q = z^3+2
z^3 + 2

julia> I = ideal(R, [p*q^2, y-z^2])
Ideal generated by
  z^8 + z^6 + 4*z^5 + 4*z^3 + 4*z^2 + 4
  y - z^2

julia> L = primary_decomposition(I)
2-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}}}:
 (Ideal (z^2 + 1, y - z^2), Ideal (z^2 + 1, y - z^2))
 (Ideal (z^6 + 4*z^3 + 4, y - z^2), Ideal (z^3 + 2, y - z^2))

julia> AL = absolute_primary_decomposition(I)
2-element Vector{Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}, Int64}}:
 (Ideal (z^2 + 1, y + 1), Ideal (z^2 + 1, y + 1), Ideal (z - _a, y + 1), 2)
 (Ideal (z^6 + 4*z^3 + 4, y - z^2), Ideal (z^3 + 2, y - z^2), Ideal (z - _a, y - _a^2), 3)

julia> AP = AL[1][3]
Ideal generated by
  z - _a
  y + 1

julia> RAP = base_ring(AP)
Multivariate polynomial ring in 2 variables y, z
  over number field of degree 2 over QQ

julia> NF = coefficient_ring(RAP)
Number field with defining polynomial x^2 + 1
  over rational field

julia> a = gen(NF)
_a

julia> minpoly(a)
x^2 + 1
```


```jldoctest
julia> R, (x, y) = graded_polynomial_ring(QQ, ["x", "y"])
(Graded multivariate polynomial ring in 2 variables over QQ, MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[x, y])

julia> I = ideal(R, [x^2+y^2])
Ideal generated by
  x^2 + y^2

julia> AL = absolute_primary_decomposition(I)
1-element Vector{Tuple{MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}, MPolyIdeal{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}, MPolyIdeal{MPolyDecRingElem{AbsSimpleNumFieldElem, AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}}, Int64}}:
 (Ideal (x^2 + y^2), Ideal (x^2 + y^2), Ideal (x + _a*y), 2)

julia> AP = AL[1][3]
Ideal generated by
  x + _a*y

julia> RAP = base_ring(AP)
Multivariate polynomial ring in 2 variables over number field graded by 
  x -> [1]
  y -> [1]
```
"""
@attr function absolute_primary_decomposition(I::MPolyIdeal{<:MPolyRingElem{QQFieldElem}})
  R = base_ring(I)
  if is_zero(I)
     return [(ideal(R, zero(R)), ideal(R, zero(R)), ideal(R, zero(R)), 1)]
  end
  (S, d) = Singular.LibPrimdec.absPrimdecGTZ(singular_polynomial_ring(I), singular_generators(I))
  decomp = d[:primary_decomp]
  absprimes = d[:absolute_primes]
  @assert length(decomp) == length(absprimes)
  V =  [(_map_last_var(R, decomp[i][1], 1, one(QQ))) for i in 1:length(decomp)]
  if length(V) == 1 && is_one(gen(V[1], 1))
    return Tuple{MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{QQMPolyRingElem}, MPolyIdeal{AbstractAlgebra.Generic.MPoly{AbsSimpleNumFieldElem}}, Int64}[]
  end 
  return [(V[i], _map_last_var(R, decomp[i][2], 1, one(QQ)),
         _map_to_ext(R, absprimes[i][1]),
         absprimes[i][2]::Int)
         for i in 1:length(decomp)]
end

@attr function absolute_primary_decomposition(
    I::MPolyIdeal{T}
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  R = base_ring(I)
  kk = coefficient_ring(R)
  R_exp, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  iso_inv(one(R))
  I_exp = ideal(R_exp, iso_inv.(gens(I)))
  res = absolute_primary_decomposition(I_exp)
  full_res = []
  for (P_ext, Q_ext, P_prime, d) in res
    @assert base_ring(P_ext) === R_exp
    P = ideal(R, unique!(iso.(gens(P_ext))))
    Q = ideal(R, unique!(iso.(gens(Q_ext))))
    RR = base_ring(P_prime)
    L = coefficient_ring(RR)
    f = defining_polynomial(L)
    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)
    LoR, to_LoR = change_base_ring(kk_ext, R)
    help_map = hom(br, LoR, iso_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
  return full_res
end

# the ideals in QQbar[x] come back in QQ[x,a] with an extra variable a added
# and the minpoly of a prepended to the ideal generator list
function _map_to_ext(Qx::MPolyRing, I::Oscar.Singular.sideal)
  Qxa = base_ring(I)
  @assert nvars(Qxa) == nvars(Qx) + 1
  p = I[1]
  minpoly = zero(Hecke.Globals.Qx)
  for (c, e) in zip(AbstractAlgebra.coefficients(p), AbstractAlgebra.exponent_vectors(p))
    setcoeff!(minpoly, e[nvars(Qxa)], QQ(c))
  end
  R, a = number_field(minpoly)
  if is_graded(Qx)
    Rx, _ = graded_polynomial_ring(R, symbols(Qx), [degree(x) for x = gens(Qx)])
  else
    Rx, _ = polynomial_ring(R, symbols(Qx))
  end
  return _map_last_var(Rx, I, 2, a)
end

# the ideals in QQ[x] also come back in QQ[x,a]
function _map_last_var(Qx::MPolyRing, I::Singular.sideal, start, a)
  newgens = elem_type(Qx)[]
  for i in start:ngens(I)
    p = I[i]
    g = MPolyBuildCtx(Qx)
    for (c, e) in zip(AbstractAlgebra.coefficients(p), AbstractAlgebra.exponent_vectors(p))
      ca = QQ(c)*a^pop!(e)
      push_term!(g, ca, e)
    end
    push!(newgens, finish(g))
  end
  return ideal(Qx, newgens)
end

#######################################################
@doc raw"""
    minimal_primes(I::MPolyIdeal; algorithm::Symbol = :GTZ)

Return a vector containing the minimal associated prime ideals of `I`.

# Implemented Algorithms

If the base ring of `I` is a polynomial ring over a field, the algorithm of
Gianni, Trager, and Zacharias is used by default (`algorithm = :GTZ`). Alternatively, characteristic sets can be
used by specifying `algorithm = :charSets`. For polynomial rings over the integers,
the algorithm proceeds as suggested by Pfister, Sadiq, and Steidel.
See [GTZ88](@cite) and [PSS11](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = minimal_primes(I)
2-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (x - y - 1)
 Ideal (x^3 - x - y^2)

julia> L = minimal_primes(I, algorithm = :charSets)
2-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (x - y - 1)
 Ideal (x^3 - x - y^2)
```
```jldoctest
julia> R, (a, b, c, d) = polynomial_ring(ZZ, ["a", "b", "c", "d"])
(Multivariate polynomial ring in 4 variables over ZZ, ZZMPolyRingElem[a, b, c, d])

julia> I = ideal(R, [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])
Ideal generated by
  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

julia> L = minimal_primes(I)
6-element Vector{MPolyIdeal{ZZMPolyRingElem}}:
 Ideal (d, c)
 Ideal (b, a)
 Ideal (2, c)
 Ideal (3)
 Ideal (13, b)
 Ideal (17, a)
```
"""
function minimal_primes(I::MPolyIdeal; algorithm::Symbol = :GTZ, cache::Bool=true)
  has_attribute(I, :minimal_primes) && return get_attribute(I, :minimal_primes)::Vector{typeof(I)}
  R = base_ring(I)
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    mp = minimal_primes(map_coefficients(pseudo_inv(mA), I); algorithm = algorithm)
    result = typeof(I)[map_coefficients(mA, x) for x = mp]
    cache && set_attribute!(I, :minimal_primes=>result)
    return result
  end
  if elem_type(base_ring(R)) <: FieldElement
    if algorithm == :GTZ
      l = Singular.LibPrimdec.minAssGTZ(singular_polynomial_ring(I), singular_generators(I))
    elseif algorithm == :charSets
      l = Singular.LibPrimdec.minAssChar(singular_polynomial_ring(I), singular_generators(I))
    else
      error("algorithm invalid")
    end
  elseif base_ring(singular_polynomial_ring(I)) isa Singular.Integers
    l = Singular.LibPrimdecint.minAssZ(singular_polynomial_ring(I), singular_generators(I))
  else
    error("base ring not implemented")
  end
  V = [ideal(R, i) for i in l]
  if length(V) == 1 && is_one(gen(V[1], 1))
    result = typeof(I)[]
    cache && set_attribute!(I, :minimal_primes=>result)
    return result
  end

  for Ptemp in V
    set_attribute!(Ptemp, :is_prime=>true)
  end

  return V
end

# rerouting the procedure for minimal primes this way leads to 
# much longer computations compared to the flattening of the coefficient
# field implemented above.
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}}
  has_attribute(I, :minimal_primes) && return get_attribute(I, :minimal_primes)::Vector{typeof(I)}

  R = base_ring(I)
  is_one(I) && return typeof(I)[]

  # Try to eliminate variables first. This will often speed up computations significantly.
  if simplify_ring
    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)
    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 = typeof(I)[ideal(R, elem_type(R)[])]
    for g in gens(I)
      K = typeof(I)[]
      is_zero(g) && continue
      for (b, k) in factor(g)
        for j in J
          push!(K, j + ideal(R, b))
        end
      end
      J = K
    end
    result = unique!(filter!(!is_one, vcat([minimal_primes(j; algorithm, factor_generators=false) for j in J]...)))
    # 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)
  result = Vector{typeof(I)}()
  for Q in dec
    push!(result, ideal(R, iso.(gens(Q))))
  end

  cache && set_attribute!(I, :minimal_primes=>result)
  return result
end

#######################################################
@doc raw"""
    equidimensional_decomposition_weak(I::MPolyIdeal)

Return a vector of equidimensional ideals where the last entry is the
equidimensional hull of `I`, that is, the intersection of the primary
components of `I` of maximal dimension. Each of the previous entries
is an ideal of lower dimension whose associated primes are exactly the associated
primes of `I` of that dimension.

# Implemented Algorithms

The implementation relies on ideas of Eisenbud, Huneke, and Vasconcelos. See [EHV92](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = equidimensional_decomposition_weak(I)
2-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (y, x)
 Ideal with 1 generator
```
"""
@attr function equidimensional_decomposition_weak(I::MPolyIdeal)
  R = base_ring(I)
  iszero(I) && return [I]
  @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field"
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    eq = equidimensional_decomposition_weak(map_coefficients(pseudo_inv(mA), I))
    retval =  typeof(I)[map_coefficients(mA, x, parent = R) for x = eq]
    for Itemp in retval
      set_attribute!(Itemp, :is_equidimensional=>true)
    end
    return retval
  end
  l = Singular.LibPrimdec.equidim(singular_polynomial_ring(I), singular_generators(I))
  V = [ideal(R, i) for i in l]
  if length(V) == 1 && is_one(gen(V[1], 1))
    return typeof(I)[]
  end

  for Itemp in V
    set_attribute!(Itemp, :is_equidimensional=>true)
  end

  return V
end

@attr function equidimensional_decomposition_weak(
    I::MPolyIdeal{T}
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  R = base_ring(I)
  R_ext, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  I_ext = ideal(R_ext, iso_inv.(gens(I)))
  res = equidimensional_decomposition_weak(I_ext)
  retval =  typeof(I)[ideal(R, unique!([x for x in iso.(gens(I)) if !iszero(x)])) for I in res]
  for Itemp in retval
    set_attribute!(Itemp, :is_equidimensional=>true)
  end
  return retval
end



@doc raw"""
    equidimensional_decomposition_radical(I::MPolyIdeal)

Return a vector of equidimensional radical ideals increasingly ordered by dimension.
For each dimension, the returned radical ideal is the intersection of the associated primes
of `I` of that dimension. 

# Implemented Algorithms

The implementation combines the algorithms of Krick and Logar (with modifications by Laplagne) and Kemper. See [KL91](@cite) and [Kem02](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = equidimensional_decomposition_radical(I)
2-element Vector{MPolyIdeal{QQMPolyRingElem}}:
 Ideal (y, x)
 Ideal (x^4 - x^3*y - x^3 - x^2 - x*y^2 + x*y + x + y^3 + y^2)
```
"""
@attr function equidimensional_decomposition_radical(I::MPolyIdeal)
  R = base_ring(I)
  @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field"
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    eq = equidimensional_decomposition_radical(map_coefficients(pseudo_inv(mA), I))
    retval = typeof(I)[map_coefficients(mA, x) for x = eq]
    for Itemp in retval
      set_attribute!(Itemp, :is_radical=>true)
      set_attribute!(Itemp, :is_equidimensional=>true)
    end
    return retval
  end
  l = Singular.LibPrimdec.prepareAss(singular_polynomial_ring(I), singular_generators(I))
  V = [ideal(R, i) for i in l]
  if length(V) == 1 && is_one(gen(V[1], 1))
    return typeof(I)[]
  end

  for Ptemp in V
    set_attribute!(Ptemp, :is_radical=>true)
    set_attribute!(Ptemp, :is_equidimensional=>true)
  end

  return V
end

@attr function equidimensional_decomposition_radical(
    I::MPolyIdeal{T}
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  R = base_ring(I)
  R_ext, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  I_ext = ideal(R_ext, iso_inv.(gens(I)))
  res = equidimensional_decomposition_weak(I_ext)
  retval = [ideal(R, unique!(iso.(gens(I)))) for I in res]
  for Itemp in retval
    set_attribute!(Itemp, :is_equidimensional=>true)
    set_attribute!(Itemp, :is_radical=>true)
  end
  return retval
end

#######################################################
@doc raw"""
    equidimensional_hull(I::MPolyIdeal)

If the base ring of `I` is a polynomial ring over a field, return the intersection
of the primary components of `I` of maximal dimension. In the case of polynomials
over the integers, return the intersection of the primary components of I of
minimal height.  If `I` is the unit ideal, return `[ideal(1)]`.

# Implemented Algorithms

For polynomial rings over a field, the implementation relies on ideas as used by
Gianni, Trager, and Zacharias or Krick and Logar. For polynomial rings over the integers,
the algorithm proceeds as suggested by Pfister, Sadiq, and Steidel. See [GTZ88](@cite),
[KL91](@cite),  and [PSS11](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = equidimensional_hull(I)
Ideal generated by
  x^5 - 2*x^4*y - 2*x^4 + x^3*y^2 + 2*x^3*y - x^2*y^2 + 2*x^2*y + 2*x^2 + 2*x*y^3 + x*y^2 - 2*x*y - x - y^4 - 2*y^3 - y^2
```
```jldoctest
julia> R, (a, b, c, d) = polynomial_ring(ZZ, ["a", "b", "c", "d"])
(Multivariate polynomial ring in 4 variables over ZZ, ZZMPolyRingElem[a, b, c, d])

julia> I = ideal(R, [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])
Ideal generated by
  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

julia> L = equidimensional_hull(I)
Ideal generated by
  3
```
"""
function equidimensional_hull(I::MPolyIdeal)
  R = base_ring(I)
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    eq = equidimensional_hull(map_coefficients(pseudo_inv(mA), I))
    return map_coefficients(mA, eq)
  end
  if elem_type(base_ring(R)) <: FieldElement
    i = Singular.LibPrimdec.equidimMax(singular_polynomial_ring(I), singular_generators(I))
  elseif base_ring(singular_polynomial_ring(I)) isa Singular.Integers
    i = Singular.LibPrimdecint.equidimZ(singular_polynomial_ring(I), singular_generators(I))
  else
    error("base ring not implemented")
  end
  return ideal(R, i)
end

function equidimensional_hull(
    I::MPolyIdeal{T}
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  R = base_ring(I)
  R_ext, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  I_ext = ideal(R_ext, iso_inv.(gens(I)))
  res = equidimensional_hull(I_ext)
  return ideal(R, unique!(iso.(gens(res))))
end

#######################################################
@doc raw"""
    equidimensional_hull_radical(I::MPolyIdeal)

Return the intersection of the associated primes of `I` of maximal dimension.
If `I` is the unit ideal, return `[ideal(1)]`.

# Implemented Algorithms

The implementation relies on a combination of the algorithms of Krick and Logar
(with modifications by Laplagne) and Kemper. See [KL91](@cite) and [Kem02](@cite).

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = intersect(ideal(R, [x, y])^2, ideal(R, [y^2-x^3+x]))
Ideal generated by
  x^3*y - x*y - y^3
  x^4 - x^2 - x*y^2

julia> I = intersect(I, ideal(R, [x-y-1])^2)
Ideal generated by
  x^5*y - 2*x^4*y^2 - 2*x^4*y + x^3*y^3 + 2*x^3*y^2 - x^2*y^3 + 2*x^2*y^2 + 2*x^2*y + 2*x*y^4 + x*y^3 - 2*x*y^2 - x*y - y^5 - 2*y^4 - y^3
  x^6 - 2*x^5 - 3*x^4*y^2 - 2*x^4*y + 2*x^3*y^3 + 3*x^3*y^2 + 2*x^3*y + 2*x^3 + 5*x^2*y^2 + 2*x^2*y - x^2 + 3*x*y^4 - 5*x*y^2 - 2*x*y - 2*y^5 - 4*y^4 - 2*y^3

julia> L = equidimensional_hull_radical(I)
Ideal generated by
  x^4 - x^3*y - x^3 - x^2 - x*y^2 + x*y + x + y^3 + y^2
```
"""
function equidimensional_hull_radical(I::MPolyIdeal)
  R = base_ring(I)
  @req coefficient_ring(R) isa AbstractAlgebra.Field "The coefficient ring must be a field"
  if isa(base_ring(R), NumField) && !isa(base_ring(R), AbsSimpleNumField)
    A, mA = absolute_simple_field(base_ring(R))
    eq = equidimensional_hull_radical(map_coefficients(pseudo_inv(mA), I))
    return map_coefficients(mA, eq)
  end
  i = Singular.LibPrimdec.equiRadical(singular_polynomial_ring(I), singular_generators(I))
  return ideal(R, i)
end

function equidimensional_hull_radical(
    I::MPolyIdeal{T}
  ) where {U<:Union{AbsSimpleNumFieldElem, <:Hecke.RelSimpleNumFieldElem}, T<:MPolyRingElem{U}}
  R = base_ring(I)
  R_ext, iso, iso_inv = _expand_coefficient_field_to_QQ(R)
  I_ext = ideal(R_ext, iso_inv.(gens(I)))
  res = equidimensional_hull_radical(I_ext)
  return ideal(R, unique!(iso.(gens(res))))
end

#######################################################
@doc raw"""
    ==(I::MPolyIdeal, J::MPolyIdeal)

Return `true` if `I` is equal to `J`, `false` otherwise.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x^2])
Ideal generated by
  x^2

julia> J = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> I == J
false
```
"""
function ==(I::MPolyIdeal, J::MPolyIdeal)
  check_base_rings(I, J)
  I === J && return true
  gens(I) == gens(J) && return true
  return issubset(I, J) && issubset(J, I)
end

### todo: wenn schon GB's  bekannt ...

#######################################################
@doc raw"""
    is_subset(I::MPolyIdeal, J::MPolyIdeal)

Return `true` if `I` is contained in `J`, `false` otherwise.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x^2])
Ideal generated by
  x^2

julia> J = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> is_subset(I, J)
true
```
"""
function is_subset(I::MPolyIdeal, J::MPolyIdeal)
  check_base_rings(I, J)
  return Singular.iszero(Singular.reduce(singular_generators(I), singular_groebner_generators(J)))
end

### todo: wenn schon GB's  bekannt ...

#######################################################

@doc raw"""
    ideal_membership(f::T, I::MPolyIdeal{T}) where T

Return `true` if `f` is contained in `I`, `false` otherwise. Alternatively, use `f in I`.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> f = x^2
x^2

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> ideal_membership(f, I)
true

julia> g = x
x

julia> g in I
false
```
"""
function ideal_membership(f::T, I::MPolyIdeal{T}; ordering::MonomialOrdering = default_ordering(base_ring(I))) where T
  Sx = singular_polynomial_ring(I, ordering)
  return Singular.iszero(Singular.reduce(Sx(f), singular_groebner_generators(I, ordering)))
end
Base.:in(f::MPolyRingElem, I::MPolyIdeal) = ideal_membership(f,I)
#######################################################
@doc raw"""
    radical_membership(f::T, I::MPolyIdeal{T}) where T

Return `true` if `f` is contained in the radical of `I`, `false` otherwise.
Alternatively, use `inradical(f, I)`.

# Examples
```jldoctest
julia> R, (x,) = polynomial_ring(QQ, ["x"])
(Multivariate polynomial ring in 1 variable over QQ, QQMPolyRingElem[x])

julia> f = x
x

julia> I = ideal(R,  [x^2])
Ideal generated by
  x^2

julia> radical_membership(f, I)
true

julia> g = x+1
x + 1

julia> inradical(g, I)
false
```
"""
function radical_membership(f::T, I::MPolyIdeal{T}) where T
  Sx = singular_polynomial_ring(I)
  return Singular.LibPolylib.rad_con(Sx(f), singular_generators(I)) == 1
end
inradical(f::MPolyRingElem, I::MPolyIdeal) = radical_membership(f,I)
################################################################################
@doc raw"""
    is_prime(I::MPolyIdeal)

Return `true` if `I` is prime, `false` otherwise.

!!! warning
    The function computes the minimal associated primes of `I`. This may take some time.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> is_prime(I)
false
```
"""
@attr Bool function is_prime(I::MPolyIdeal)
  D = minimal_primes(I)
  return length(D) == 1 && issubset(D[1], I)
end

################################################################################
@doc raw"""
    is_primary(I::MPolyIdeal)

Return `true` if `I` is primary, `false` otherwise.

!!! warning
    The function computes a primary decomposition of `I`. This may take some time.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> is_primary(I)
true
```
"""
@attr Bool function is_primary(I::MPolyIdeal)
  if has_attribute(I, :is_prime) && is_prime(I)
    return true
  end

  D = primary_decomposition(I)
  return length(D) == 1
end

#######################################################
@doc raw"""
    base_ring(I::MPolyIdeal)

Return the ambient ring of `I`.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> base_ring(I)
Multivariate polynomial ring in 2 variables x, y
  over rational field
```
"""
function base_ring(I::MPolyIdeal{S}) where {S}
  return I.gens.Ox::parent_type(S)
end

base_ring_type(::Type{MPolyIdeal{S}}) where {S} = parent_type(S)

#######################################################
@doc raw"""
    coefficient_ring(I::MPolyIdeal)

Return the coefficient ring of `I`, which is the coefficient ring of the
polynomial ring containing the ideal.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> coefficient_ring(I)
Rational field
```
"""
coefficient_ring(I::MPolyIdeal) = coefficient_ring(base_ring(I))


#######################################################
@doc raw"""
    number_of_generators(I::MPolyIdeal)

Return the number of generators of `I`.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> number_of_generators(I)
3
```
"""
function number_of_generators(I::MPolyIdeal)
  return length(I.gens)
end

#######################################################
@doc raw"""
    gens(I::MPolyIdeal)

Return the generators of `I`.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, y])^2
Ideal generated by
  x^2
  x*y
  y^2

julia> gens(I)
3-element Vector{QQMPolyRingElem}:
 x^2
 x*y
 y^2
```
"""
function gens(I::MPolyIdeal)
  return [I.gens[Val(:O), i] for i=1:ngens(I)]
end

gen(I::MPolyIdeal, i::Int) = I.gens[Val(:O), i]

#######################################################
@doc raw"""
    dim(I::MPolyIdeal)

Return the Krull dimension of `I`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [y-x^2, x-z^3])
Ideal generated by
  -x^2 + y
  x - z^3

julia> dim(I)
1
```
"""
@attr Int function dim(I::MPolyIdeal)
  if I.dim > -1
    return I.dim
  end
  is_zero(ngens(base_ring(I))) && return 0 # Catch a boundary case
  I.dim = Singular.dimension(singular_groebner_generators(I, false, true))
  return I.dim
end

#######################################################
#######################################################
@doc raw"""
    codim(I::MPolyIdeal)

Return the codimension of `I`.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"])
(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x, y, z])

julia> I = ideal(R, [y-x^2, x-z^3])
Ideal generated by
  -x^2 + y
  x - z^3

julia> codim(I)
2
```
"""
codim(I::MPolyIdeal{T}) where {T<:MPolyRingElem{<:FieldElem}}= nvars(base_ring(I)) - dim(I)
codim(I::MPolyIdeal) = dim(base_ring(I)) - dim(I)

# Some fixes which were necessary for the above
dim(R::MPolyRing) = dim(base_ring(R)) + nvars(R)
dim(R::ZZRing) = 1


################################################################################
#
#  iszero and isone functions
#
################################################################################

@doc raw"""
    is_zero(I::MPolyIdeal)

Return `true` if `I` is the zero ideal, `false` otherwise.
# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, y-x^2)
Ideal generated by
  -x^2 + y

julia> is_zero(I)
false
```
"""
@attr Bool function is_zero(I::MPolyIdeal)
  lg = gens(I)
  return isempty(lg) || all(iszero, lg)
end

@doc raw"""
    is_one(I::MPolyIdeal)

Return `true` if `I` is generated by `1`, `false` otherwise.
# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> I = ideal(R, [x, x + y, y - 1])
Ideal generated by
  x
  x + y
  y - 1

julia> is_one(I)
true
```
"""
@attr Bool function is_one(I::MPolyIdeal)
  R = base_ring(I)
  if iszero(I)
      return false
  end
  if any(x -> (!is_zero(x) && is_constant(x) && is_unit(first(AbstractAlgebra.coefficients(x)))), gens(I))
    return true
  end
  gb = groebner_basis(I, complete_reduction = true)
  return is_constant(gb[1]) && is_unit(first(AbstractAlgebra.coefficients(gb[1])))
end

################################################################################
#
#  is_monomial functions
#
################################################################################

@doc raw"""
    is_monomial(f::MPolyRingElem)

Return `true` if `f` is a monomial, `false` otherwise.

    is_monomial(I::MPolyIdeal)

Return `true` if `I` can be generated by monomials, `false` otherwise.

# Examples
```jldoctest
julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> f = 2*x+y
2*x + y

julia> g = y
y

julia> is_monomial(f)
false

julia> is_monomial(g)
true

julia> is_monomial(ideal(R, [f, g]))
true
```
"""
function is_monomial(f::MPolyRingElem)
   if (length(f) == 1 && coeff(f, 1) == 1)
      return true
   else
      return false
   end
end

function _ismonomial(V::Vector{<: MPolyRingElem})
  return all(is_monomial, V)
end

@attr Bool function is_monomial(I::MPolyIdeal)
  if _ismonomial(gens(I))
    return true
  end
  GB = gens(groebner_basis(I, complete_reduction = true))
  if _ismonomial(GB)
    return true
  end
  return false
end

################################################################################
#
# Small generating set
#
################################################################################

@doc raw"""
    small_generating_set(I::MPolyIdeal)

Given a ideal `I` in a multivariate polynomial ring over a field,
return an array containing a set of generators of `I`, which is usually
smaller than the original one.
If `I` is the zero ideal an empty list is returned.

!!! note
   Minimal generating sets exist only in the local and the homogeneous case. Beyond these cases, the best one can hope for is some small set of generators.

# Examples
```jldoctest
julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);

julia> V = [x, z^2, x^3+y^3, y^4, y*z^5];

julia> I = ideal(R, V)
Ideal generated by
  x
  z^2
  x^3 + y^3
  y^4
  y*z^5

julia> small_generating_set(I)
3-element Vector{QQMPolyRingElem}:
 x
 z^2
 y^3

julia> WW = [x-y^3,z-y^5,x^2-z^3, y^6-y^15]  # certainly not homogeneous
4-element Vector{QQMPolyRingElem}:
 x - y^3
 -y^5 + z
 x^2 - z^3
 -y^15 + y^6

julia> J = ideal(R,WW)
Ideal generated by
  x - y^3
  -y^5 + z
  x^2 - z^3
  -y^15 + y^6

julia> small_generating_set(J)
3-element Vector{QQMPolyRingElem}:
 -x^2 + z^3
 -x + y^3
 x*y^2 - z
```
"""
function small_generating_set(
    I::MPolyIdeal{T};
    algorithm::Symbol=:simple
  ) where {T<:MPolyRingElem{<:FieldElem}}
  # For non-homogeneous ideals, we do not have a notion of minimal generating
  # set, but Singular.mstd still provides a good heuristic to find a small
  # generating set.

  R = base_ring(I)

  # in the ungraded case, mstd's heuristic returns smaller gens when recomputing gb
  sing_gb, sing_min = Singular.mstd(singular_generators(I))
  ring = I.gens.Ox
  computed_gb = IdealGens(ring, sing_gb, true)
  if !haskey(I.gb,computed_gb.ord)
  # if not yet present, store gb for later use
    I.gb[computed_gb.ord] = computed_gb
  end

  # we do not have a notion of minimal generating set in this context!
  # If we are unlucky, mstd can even produce a larger generating set
  # than the original one!!!
  return_value = filter(!iszero, (R).(gens(sing_min)))

  # The following is a common phenomenon which we can not fully explain yet. So far nothing but a 
  # restart really seems to help, unfortunately.
  if is_zero(length(return_value))
    !is_zero(I) && error("singular crashed in the background; please restart your session!")
  end
  if length(return_value) <= ngens(I)
    return return_value
  else
    return gens(I)
  end
end

# in graded rings, reusing a cached gb makes sense, so use minimal_generating_set there
small_generating_set(I::MPolyIdeal{<:MPolyDecRingElem}; algorithm::Symbol=:simple) = minimal_generating_set(I)

################################################################################
#
# Grassmann Pluecker Ideal
#
################################################################################
#returns Pluecker ideal in ring with standard grading
function grassmann_pluecker_ideal(subspace_dimension::Int, ambient_dimension::Int)
  I = convert(MPolyIdeal{QQMPolyRingElem},
              Polymake.ideal.pluecker_ideal(subspace_dimension, ambient_dimension))
  base, _ = grade(base_ring(I))
  o = degrevlex(base)

  return ideal(IdealGens(base,base.(gens(I)), o;
                   keep_ordering=true,
                   isReduced=true,
                   isGB=true))
end

@doc raw"""
    grassmann_pluecker_ideal([ring::MPolyRing,] subspace_dimension::Int, ambient_dimension::Int)

Given a (possibly graded) ring, an ambient dimension, and a subspace dimension, return the ideal in the ring
generated by the Plücker relations. If the input ring is not graded, return the ideal in the ring with the standard grading.
If the ring is not specified return the ideal in a multivariate polynomial ring over the rationals with the standard grading.

The Grassmann-Plücker ideal is the homogeneous ideal generated by the relations defined by 
the Plücker Embedding of the Grassmannian. That is given Gr$(k, n)$ the Moduli 
space of all $k$-dimensional subspaces of an $n$-dimensional vector space, the relations
are given by all $d \times d$ minors of a $d \times n$ matrix. For the algorithm see
[Stu93](@cite).

# Examples
```jldoctest
julia> grassmann_pluecker_ideal(2, 4)
Ideal generated by
  x[1]*x[6] - x[2]*x[5] + x[3]*x[4]

julia> R, x = polynomial_ring(residue_ring(ZZ, 7)[1], "x" => (1:2, 1:3))
(Multivariate polynomial ring in 6 variables over ZZ/(7), zzModMPolyRingElem[x[1, 1] x[1, 2] x[1, 3]; x[2, 1] x[2, 2] x[2, 3]])

julia> grassmann_pluecker_ideal(R, 2, 4)
Ideal generated by
  x[1, 1]*x[2, 3] + 6*x[2, 1]*x[1, 3] + x[1, 2]*x[2, 2]
```
"""
function grassmann_pluecker_ideal(ring::MPolyRing,
                                 subspace_dimension::Int,
                                 ambient_dimension::Int) 
  I = grassmann_pluecker_ideal(subspace_dimension, ambient_dimension)
  coeff_ring = base_ring(ring)
  o = degrevlex(ring)
  coeffmap = c -> begin
    @assert is_integral(c)
    return coeff_ring(numerator(c))
  end
  if !is_graded(ring)
	 ring, _ = grade(ring)
  end
  h = hom(base_ring(I), ring, coeffmap, gens(ring))
  converted_generators = elem_type(ring)[h(g) for g in groebner_basis(I; ordering = degrevlex(base_ring(I)))]
  ideal(IdealGens(ring, converted_generators, o;
                   keep_ordering=true,
                   isReduced=true,
                   isGB=true))
end

################################################################################
#
#  Printing
#
################################################################################

function AbstractAlgebra.expressify(a::MPolyIdeal; context = nothing)
  return Expr(:call, :ideal, [expressify(g, context = context) for g in collect(a.gens)]...)
end

# See Rings/mpolyquo-localizations.jl for show functions

################################################################################
#
#  Pluecker ideal
#
################################################################################

_pluecker_sgn(a::Vector{Int}, b::Vector{Int}, t::Int)::Int =
 iseven(count(z -> z > t, a) + count(z -> z < t, b)) ? 1 : -1
@doc raw"""
    flag_pluecker_ideal(F::Union{Field, MPolyRing}, dimensions::Vector{Int}, n::Int; minimal::Bool=true)

Returns the generators of the defining ideal for the complete flag variety 
$\text{Fl}(\mathbb{F}, (d_1,\dots,d_k), n)$, where $(d_1,\dots,d_k)
=$`dimensions`, with $d_j\leq n-1$,  denotes the rank.  That is, the vanishing
set of this ideal corresponds to the space of $k$-step flags of linear
subspaces $V_1\subset\dots\subset V_k$ in $\mathbb{F}^n$, where
$\text{dim}(V_j) = d_{j}$.  You can obtain the generators for the
*complete flag variety* of $\mathbb{F}^{n}$ by taking `dimensions`
$=(1,\dots,n-1)$ and `n`$=n$.  We remark that evaluating for `F = QQ` yields
the same set of generators as any field of characteristic $0$.

The first parameter can either be $\mathbb{F}$, or a polynomial ring over $\mathbb{F}$,
with $\sum^{k}_{j=1}{n\choose d_j}$ variables.
The parameter `dimensions` needs to be a vector of distinct increasing entries.
Evaluating this function with the parameter `minimal = true` returns the reduced Gröbner basis for
the flag Plücker ideal with respect to the degree reverse lexicographical order. For more details, see Theorem 14.6 [MS05](@cite)
# Examples

Complete flag variety $\text{Fl}(\mathbb{Q}, (1,2,3), 4)$. 
```jldoctest
julia> flag_pluecker_ideal(QQ,[1,2,3],4)
Ideal generated by
  x[[2, 4]]*x[[1, 2, 3]] - x[[2, 3]]*x[[1, 2, 4]] + x[[1, 2]]*x[[2, 3, 4]]
  x[[1, 4]]*x[[1, 2, 3]] - x[[1, 3]]*x[[1, 2, 4]] + x[[1, 2]]*x[[1, 3, 4]]
  -x[[3, 4]]*x[[1, 2, 3]] - x[[1, 3]]*x[[2, 3, 4]] + x[[2, 3]]*x[[1, 3, 4]]
  -x[[1, 4]]*x[[2, 3, 4]] + x[[2, 4]]*x[[1, 3, 4]] - x[[3, 4]]*x[[1, 2, 4]]
  -x[[1]]*x[[2, 3, 4]] + x[[2]]*x[[1, 3, 4]] - x[[3]]*x[[1, 2, 4]] + x[[4]]*x[[1, 2, 3]]
  -x[[1, 4]]*x[[2, 3]] + x[[2, 4]]*x[[1, 3]] - x[[3, 4]]*x[[1, 2]]
  -x[[1]]*x[[2, 3]] + x[[2]]*x[[1, 3]] - x[[3]]*x[[1, 2]]
  -x[[2]]*x[[3, 4]] + x[[3]]*x[[2, 4]] - x[[4]]*x[[2, 3]]
  -x[[1]]*x[[3, 4]] + x[[3]]*x[[1, 4]] - x[[4]]*x[[1, 3]]
  -x[[1]]*x[[2, 4]] + x[[2]]*x[[1, 4]] - x[[4]]*x[[1, 2]]

```
Flag variety $\text{Fl}(\mathbb{Q},(1,3),4)$.
```jldoctest
julia> flag_pluecker_ideal(QQ,[1,3],4)
Ideal generated by
  -x[[1]]*x[[2, 3, 4]] + x[[2]]*x[[1, 3, 4]] - x[[3]]*x[[1, 2, 4]] + x[[4]]*x[[1, 2, 3]]
```
An example with a custom ring as input.
```jldoctest
julia> R, _ = polynomial_ring(QQ, 8)
(Multivariate polynomial ring in 8 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4, x5, x6, x7, x8])

julia> flag_pluecker_ideal(R, [1,3], 4; minimal=false)
Ideal generated by
  x1*x6 - x2*x5 + x3*x7 - x4*x8
```
"""
function flag_pluecker_ideal(F::Field, dimensions::Vector{Int}, n::Int; minimal::Bool=true)
  L = reduce(vcat, [subsets(n, d) for d in dimensions])
  ring, _ = polynomial_ring(F, "x" => L; cached=false)

  return flag_pluecker_ideal(ring, dimensions, n; minimal=minimal)

end

flag_pluecker_ideal(dimensions::Vector{Int}, n::Int; minimal::Bool=true) = flag_pluecker_ideal(QQ,dimensions,n, minimal=minimal)

function flag_pluecker_ideal(ring::MPolyRing{<: FieldElem}, dimensions::Vector{Int}, n::Int; minimal::Bool=true)
  dimensions = unique(dimensions)
   
  @req maximum(dimensions) <= n-1 "The dimensions must be at most n-1"
  @req sum(binomial(n, d) for d in dimensions) == ngens(ring) "The ring doesn't have the right number of variables"
  @req issorted(dimensions) "The dimensions must be increasing"

  x = gens(ring)
  L = reduce(vcat, [subsets(n, d) for d in dimensions])

  xdict = Dict{Vector{Int},MPolyRingElem}([L[i] => x[i] for i in 1:length(L)])

  L1 = reduce(vcat, [subsets(n, d - 1) for d in dimensions])
  L2 = reduce(vcat, [subsets(n, d + 1) for d in dimensions])

  X = [(a, b) for a in L1, b in L2 if length(b) - length(a) >= 2]

  T = [
    [
      (_pluecker_sgn(a, b, t)) * xdict[sort(union(a, t))] * xdict[setdiff(b, t)] for
      t in setdiff(b, a)
    ] for (a, b) in X
  ]

  genz = [sum(t) for t in T]
  filter!(!iszero, genz)
  gg = unique(x ->  divexact(x, canonical_unit(x)), genz)
  I  = ideal(gg)

  !minimal && return I 

  o = degrevlex(ring)
  converted_generators = collect(groebner_basis(I; ordering = o))
  return ideal(IdealGens(ring, converted_generators, o;
                   keep_ordering=true,
                   isReduced=true,
                   isGB=true))
end