Ideal saturation: Resolving old question found in code (#4238)

* Ideal saturation: Resolving old question found in code

* Improvement, backwards compatibility

* correction

docs/src/CommutativeAlgebra/ideals.md

@@ -134,7 +134,7 @@ Given two ideals $I, J$ of a ring $R$, the saturation of $I$ with respect to $J$
 $I:J^{\infty} = \bigl\{ f \in R \:\big|\: f J^k \!\subset I {\text{ for some }}k\geq 1 \bigr\} = \textstyle{\bigcup\limits_{k=1}^{\infty} (I:J^k)}.$
 
 ```@docs
-saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T
+saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T}; iteration::Bool=false) where T
 saturation_with_index(I::MPolyIdeal{T}, J::MPolyIdeal{T}) where T
 ```
 

src/Rings/mpoly-ideals.jl

@@ -302,11 +302,16 @@ end
 
 # saturation #######################################################
 @doc raw"""
-    saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T
+    saturation(I::MPolyIdeal{T}, 
+               J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I))); 
+               iteration::Bool=false) 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.
 
+If `iteration` is set to `true`, the saturation is done by carrying out successive ideal quotient computations as in 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.
+
 # Examples
 ```jldoctest
 julia> R, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z])
@@ -338,27 +343,33 @@ 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))
+function saturation(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I))); iteration::Bool=false) where T
+  if iteration
+    K, _ = Singular.saturation(singular_generators(I), singular_generators(J))
+  else
+    K, _ = Singular.saturation2(singular_generators(I), singular_generators(J))
+  end
   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
+_saturation2(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T = saturation(I, J)
+# kept for backwards compatibility
 
 #######################################################
 @doc raw"""
-    saturation_with_index(I::MPolyIdeal{T}, J::MPolyIdeal{T} = ideal(base_ring(I), gens(base_ring(I)))) where T
+    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}$ (saturation index).
 
-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.
 
+!!! note
+    If the saturation index is not needed, we recommend to use `saturation(I, J)` which is typically faster.
+
 # Examples
 ```jldoctest
 julia> R, (x, y, z) = polynomial_ring(QQ, [:x, :y, :z])
@@ -386,6 +397,7 @@ julia> K, m = saturation_with_index(I, J)
 
 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