Saturation for modules

docs/src/CommutativeAlgebra/ModulesOverMultivariateRings/subquotients.md

@@ -363,9 +363,18 @@ quotient(M::SubquoModule{T}, N::SubquoModule{T}) where T
 ```
 
 ```@docs
-quotient(M::SubquoModule{T}, J::MPolyIdeal{T}) where T
+quotient(M::SubquoModule, J::Ideal)
 ```
 
+```@docs
+saturation(M:: SubquoModule, J::Ideal)
+```
+
+```@docs
+saturation_with_index(M:: SubquoModule, J::Ideal)
+```
+
+
 ## Submodules and Quotients
 
 ```@docs

docs/src/CommutativeAlgebra/affine_algebras.md

@@ -296,7 +296,6 @@ minimal_generating_set(I::MPolyQuoIdeal{<:MPolyDecRingElem})
 
 ```@docs
 intersect(a::MPolyQuoIdeal{T}, bs::MPolyQuoIdeal{T}...) where T
-intersect(V::Vector{MPolyQuoIdeal{T}}) where T
 ```
 
 #### Ideal Quotients

src/Modules/UngradedModules/SubquoModule.jl

@@ -1779,8 +1779,7 @@
 (::Colon)(M::SubquoModule{T}, N::SubquoModule{T}) where T = quotient(M, N)
 
 @doc raw"""
-     quotient(M::SubquoModule{T}, J::MPolyIdeal{T}) where T
-
+     quotient(M::SubquoModule, J::Ideal)
 
 Return the quotient of `M` by `J`.
 
@@ -1819,12 +1818,10 @@
   z^2
 
 julia> L = quotient(M, J)
-Subquotient of Submodule with 5 generators
+Subquotient of Submodule with 3 generators
 1 -> x*e[1]
-2 -> y^3*e[1]
-3 -> z^4*e[1]
-4 -> y*z^3*e[1]
-5 -> y^2*z^2*e[1]
+2 -> y*z^3*e[1]
+3 -> y^2*z^2*e[1]
 by Submodule with 3 generators
 1 -> x^2*e[1]
 2 -> y^3*e[1]
@@ -1833,32 +1830,286 @@
 julia> ambient_free_module(L) == ambient_free_module(M)
 true
 
+```
+
+```jldoctest
+julia> S, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z]);
+
+julia> R, _ = quo(S, ideal(S, [x+y+z]));
+
+julia> F = free_module(R, 1);
+
+julia> AM = R[x;];
+
+julia> BM = R[x^2; y^3; z^4];
+
+julia> M = SubquoModule(F, AM, BM)
+Subquotient of Submodule with 1 generator
+1 -> (-y - z)*e[1]
+by Submodule with 3 generators
+1 -> (y^2 + 2*y*z + z^2)*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+julia> J = ideal(R, [x, y, z])^2
+Ideal generated by
+  x^2
+  x*y
+  x*z
+  y^2
+  y*z
+  z^2
+
+julia> quotient(M, J)
+Subquotient of Submodule with 2 generators
+1 -> z*e[1]
+2 -> y*e[1]
+by Submodule with 3 generators
+1 -> (y^2 + 2*y*z + z^2)*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
 ```
 """
-function quotient(M::SubquoModule{T}, J::MPolyIdeal{T}) where T
+function quotient(M::SubquoModule, J::Ideal)
+  @assert base_ring(M) == base_ring(J)
+  M_quo = isdefined(M, :quo) ? M.quo : SubModuleOfFreeModule(ambient_free_module(M), Vector{elem_type(ambient_free_module(M))}())
+  U = M.sub+M.quo
+  UF = _quotient(U, J)
+  res = SubquoModule(UF)
+  res.quo = M.quo
+  return simplify_light(res)[1]
+end
+
+(::Colon)(M::SubquoModule, J::Ideal) = quotient(M, N)
+
+function _quotient(U::SubModuleOfFreeModule, J::Ideal) ### TODO Replace by generic method
   error("not implemented for the given types of modules.")
 end
 
-(::Colon)(M::SubquoModule{T}, J::MPolyIdeal{T}) where T = quotient(M, N)
+function _quotient(U::SubModuleOfFreeModule, J::Ideal{T}) where T <: Union{MPolyRingElem, MPolyQuoRingElem}
+  F = ambient_free_module(U)
+  SgU = singular_generators(U.gens)
+  SgJ = singular_generators(J.gens)
+  SQ = Singular.quotient(SgU, SgJ)
+  MG = ModuleGens(F, SQ)
+  return SubModuleOfFreeModule(F, MG)
+end
 
-function quotient(M::SubquoModule{T}, J::MPolyIdeal{T}) where T <: MPolyRingElem
+########################################
+### saturation for modules
+########################################
+
+@doc raw"""
+     saturation(M::SubquoModule,
+               J::Ideal = ideal(base_ring(M), gens(base_ring(M)));
+               iteration::Bool = false)
+
+Return the saturation $M:J^{\infty}$ of `M` with respect to `J`.
+
+If the ideal `J` is not given, the ideal generated by the generators (variables) of `base_ring(M)` is used.
+
+Setting `iteration` to `true` only has an effect over rings of type `MPolyRing` or `MPolyQuoRing`. Over such rings,
+if `iteration` is set to `true`, the saturation is done by carrying out successive ideal quotient computations as
+suggested by the definition of saturation. Otherwise, a more sophisticated Gröbner basis approach is used which is typically
+faster. Applying the two approaches may lead to different generating sets of the saturation.
+
+!!! note
+    By definition,
+    $M:J^{\infty} = \{ a \in A \mid J^ka \subset M \text{ for some } k\geq 1 \} = \bigcup\limits_{k=1}^{\infty} (M:J^k).$
+    Here, $A$ is the ambient module of $M$.
+
+# Examples
+```jldoctest
+julia> R, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z]);
+
+julia> F = free_module(R, 1);
+
+julia> AM = R[x;];
+
+julia> BM = R[x^2; y^3; z^4];
+
+julia> M = SubquoModule(F, AM, BM)
+Subquotient of Submodule with 1 generator
+1 -> x*e[1]
+by Submodule with 3 generators
+1 -> x^2*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+julia> J = ideal(R, [x, y, z])^2
+Ideal generated by
+  x^2
+  x*y
+  x*z
+  y^2
+  y*z
+  z^2
+
+julia> saturation(M, J)
+Subquotient of Submodule with 1 generator
+1 -> e[1]
+by Submodule with 3 generators
+1 -> x^2*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+```
+"""
+function saturation(M::SubquoModule, J::Ideal = ideal(base_ring(M), gens(base_ring(M))); iteration::Bool = false)
   @assert base_ring(M) == base_ring(J)
-  F = ambient_free_module(M)
+  M_quo = isdefined(M, :quo) ? M.quo : SubModuleOfFreeModule(ambient_free_module(M), Vector{elem_type(ambient_free_module(M))}())
   U = M.sub+M.quo
+  UF = _saturation(U, J; iteration = iteration)
+  res = SubquoModule(UF)
+  res.quo = M.quo
+  return simplify_light(res)[1]
+end
+
+function _saturation(U::SubModuleOfFreeModule, J::Ideal) ### TODO Replace by generic method
+  error("not implemented for the given types of modules.")
+end
+
+function _saturation(U::SubModuleOfFreeModule, J::Ideal{T}; iteration::Bool = false) where T <: Union{MPolyRingElem, MPolyQuoRingElem}
+  F = ambient_free_module(U)
   SgU = singular_generators(U.gens)
-  singular_assure(J)
-  SQ = Singular.quotient(SgU, J.gens.S)
+  SgJ = singular_generators(J.gens)
+  SQ, _ = Singular.saturation(SgU, SgJ)
   MG = ModuleGens(F, SQ)
-  UF = SubModuleOfFreeModule(F, MG)
+  return SubModuleOfFreeModule(F, MG)
+end
+
+@doc raw"""
+     saturation_with_index(M::SubquoModule,
+               J::MPolyIdeal = ideal(base_ring(M), gens(base_ring(M)));
+               iteration::Bool = false)
+
+
+Return the saturation $M:J^{\infty}$ of $M$ with respect to $J$ and the smallest integer $k$ such that $I:J^k = I:J^{\infty}$ (saturation index).
+
+If the ideal `J` is not given, the ideal generated by the generators (variables) of `base_ring(M)` is used.
+
+Setting `iteration` to `true` only has an effect over rings of type `MPolyRing` or `MPolyQuoRing`. Over such rings,
+if `iteration` is set to `true`, the saturation is done by carrying out successive module quotient computations as
+suggested by the definition of saturation. Otherwise, a more sophisticated Gröbner basis approach is used which is typically
+faster. Applying the two approaches may lead to different generating sets of the saturation.
+
+!!! note
+    By definition,
+    $M:J^{\infty} = \{ a \in A \mid J^ka \subset M \text{ for some } k\geq 1 \} = \bigcup\limits_{k=1}^{\infty} (M:J^k).$
+    Here, $A$ is the ambient module of $M$.
+
+# Examples
+```jldoctest
+julia> R, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z]);
+
+julia> F = free_module(R, 1);
+
+julia> AM = R[x;];
+
+julia> BM = R[x^2; y^3; z^4];
+
+julia> M = SubquoModule(F, AM, BM)
+Subquotient of Submodule with 1 generator
+1 -> x*e[1]
+by Submodule with 3 generators
+1 -> x^2*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+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 = saturation_with_index(M, J);
+
+julia> L[1]
+Subquotient of Submodule with 1 generator
+1 -> e[1]
+by Submodule with 3 generators
+1 -> x^2*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+julia> L[2]
+3
+
+```
+
+```jldoctest
+julia> S, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z]);
+
+julia> R, _ = quo(S, ideal(S, [x+y+z]));
+
+julia> F = free_module(R, 1);
+
+julia> AM = R[x;];
+
+julia> BM = R[x^2; y^3; z^4];
+
+julia> M = SubquoModule(F, AM, BM)
+Subquotient of Submodule with 1 generator
+1 -> (-y - z)*e[1]
+by Submodule with 3 generators
+1 -> (y^2 + 2*y*z + z^2)*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+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 = saturation_with_index(M, J);
+
+julia> L[1]
+Subquotient of Submodule with 1 generator
+1 -> e[1]
+by Submodule with 3 generators
+1 -> (y^2 + 2*y*z + z^2)*e[1]
+2 -> y^3*e[1]
+3 -> z^4*e[1]
+
+julia> L[2]
+2
+
+```
+"""
+function saturation_with_index(M::SubquoModule, J::Ideal = ideal(base_ring(M), gens(base_ring(M))); iteration::Bool = false)
+@assert base_ring(M) == base_ring(J)
+  M_quo = isdefined(M, :quo) ? M.quo : SubModuleOfFreeModule(ambient_free_module(M), Vector{elem_type(ambient_free_module(M))}())
+  U = M.sub+M.quo
+  UF, k = _saturation_with_index(U, J; iteration = iteration)
   res = SubquoModule(UF)
   res.quo = M.quo
-  return res
+  return simplify_light(res)[1], k
 end
 
+function _saturation_with_index(U::SubModuleOfFreeModule, J::Ideal) ### TODO Replace by generic method
+  error("not implemented for the given types of modules.")
+end
 
+function _saturation_with_index(U::SubModuleOfFreeModule, J::Ideal{T}; iteration::Bool = false) where T <: Union{MPolyRingElem, MPolyQuoRingElem}
+  F = ambient_free_module(U)
+  SgU = singular_generators(U.gens)
+  SgJ = singular_generators(J.gens)
+  SQ, k = Singular.saturation(SgU, SgJ)
+  MG = ModuleGens(F, SQ)
+  return SubModuleOfFreeModule(F, MG), k
+end
 
 ########################################
 
+
 @doc raw"""
     represents_element(a::FreeModElem, SQ::SubquoModule)