Speed up polynomial mappings (oscar-system#4124)

*  speed up mapping of polynomials.

src/AlgebraicGeometry/Schemes/AffineSchemeOpenSubscheme/Rings/Methods.jl

@@ -506,3 +506,7 @@ function Base.show(io::IO, ::MIME"text/plain", a::AffineSchemeOpenSubschemeRingE
   end
 end
 
+# overwrite a method for mapping of rings which would throw otherwise
+function _allunique(lst::Vector{T}) where {T<:MPolyQuoRingElem{<:MPolyRingElem{<:AffineSchemeOpenSubschemeRingElem}}}
+  return false
+end

src/AlgebraicGeometry/Schemes/AffineSchemes/Morphisms/Attributes.jl

@@ -46,7 +46,7 @@ Affine scheme morphism
   from [x1, x2, x3]  scheme(x1)
   to   [x1, x2, x3]  affine 3-space over QQ
 given by the pullback function
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 
@@ -95,7 +95,7 @@ Affine scheme morphism
   from [x1, x2, x3]  scheme(x1)
   to   [x1, x2, x3]  affine 3-space over QQ
 given by the pullback function
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 
@@ -143,7 +143,7 @@ Ring homomorphism
   from multivariate polynomial ring in 3 variables over QQ
   to quotient of multivariate polynomial ring by ideal (x1)
 defined by
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 ```
@@ -254,7 +254,7 @@ Affine scheme morphism
   from [x1, x2, x3]  scheme(x1)
   to   [x1, x2, x3]  affine 3-space over QQ
 given by the pullback function
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 

src/AlgebraicGeometry/Schemes/AffineSchemes/Morphisms/Constructors.jl

@@ -132,7 +132,7 @@ Affine scheme morphism
   from [x1, x2, x3]  scheme(x1)
   to   [x1, x2, x3]  affine 3-space over QQ
 given by the pullback function
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 
@@ -181,7 +181,7 @@ Affine scheme morphism
   from [x1, x2, x3]  scheme(x1)
   to   [x1, x2, x3]  affine 3-space over QQ
 given by the pullback function
-  x1 -> 0
+  x1 -> x1
   x2 -> x2
   x3 -> x3
 

src/AlgebraicGeometry/Schemes/AffineSchemes/Objects/Attributes.jl

@@ -234,7 +234,7 @@ Affine scheme morphism
   from [x, y]  scheme(x)
   to   [x, y]  affine 2-space over QQ
 given by the pullback function
-  x -> 0
+  x -> x
   y -> y
 
 julia> inc == inclusion_morphism(Y, X)

src/Rings/MPolyMap/MPolyAnyMap.jl

@@ -42,16 +42,32 @@ const _DomainTypes = Union{MPolyRing, MPolyQuoRing}
   coeff_map::U
   img_gens::Vector{V}
   temp_ring           # temporary ring used when evaluating maps
+  variable_indices::Vector{Int} # a table where the i-th entry contains the 
+                                # index of the variable where it is mapped to
+                                # in case the mapping takes such a particularly
+                                # simple form.
 
   function MPolyAnyMap{D, C, U, V}(domain::D,
-                                codomain::C,
-                                coeff_map::U,
-                                img_gens::Vector{V}) where {D, C, U, V}
-      @assert V === elem_type(C)
-      for g in img_gens
-        @assert parent(g) === codomain "elements does not have the correct parent"
-      end
-    return new{D, C, U, V}(domain, codomain, coeff_map, img_gens)
+                                   codomain::C,
+                                   coeff_map::U,
+                                   img_gens::Vector{V};
+                                   check_for_mapping_of_vars::Bool=true
+                                  ) where {D, C, U, V}
+    @assert V === elem_type(C)
+    for g in img_gens
+      @assert parent(g) === codomain "elements does not have the correct parent"
+    end
+    result = new{D, C, U, V}(domain, codomain, coeff_map, img_gens)
+    # If it ever turned out that doing the checks within the following if-block
+    # is a bottleneck, consider passing on the `check_for_mapping_of_vars` kw 
+    # argument to the outer constructors or make the outer constructors 
+    # call the inner one with this argument set to `false`. This way the check 
+    # can safely be disabled.
+    if check_for_mapping_of_vars && all(_is_gen, img_gens) && _allunique(img_gens)
+      gens_codomain = gens(codomain)
+      result.variable_indices = [findfirst(==(x), gens_codomain) for x in img_gens]
+    end
+    return result
   end
 end
 

src/Rings/MPolyMap/MPolyRing.jl

@@ -132,17 +132,111 @@ end
 #
 ################################################################################
 
+# Some additional methods needed for the test in the constructor for MPolyAnyMap
+_is_gen(x::MPolyQuoRingElem) = _is_gen(lift(x))
+_is_gen(x::MPolyDecRingElem) = is_gen(forget_grading(x))
+_is_gen(x::MPolyRingElem) = is_gen(x)
+# default method; overwrite if you want this to work for your rings.
+_is_gen(x::NCRingElem) = false
+
+# In case there is a type of ring elements for which hashing is correctly implemented 
+# and does not throw an error, this gives the opportunity to overwrite the `allunique`
+# to be used within the constructor for maps. 
+function _allunique(lst::Vector{T}) where {T<:MPolyRingElem}
+  return allunique(lst)
+end
+
+# We have a lot of rings which do/can not implement correct hashing. 
+# So we make the following the default.
+function _allunique(lst::Vector{T}) where {T<:RingElem}
+  return all(!(x in lst[i+1:end]) for (i, x) in enumerate(lst))
+end
+
+function _build_poly(u::MPolyRingElem, indices::Vector{Int}, S::MPolyRing)
+  kk = coefficient_ring(S)
+  r = ngens(S)
+  ctx = MPolyBuildCtx(S)
+  for (c, e) in zip(AbstractAlgebra.coefficients(u), AbstractAlgebra.exponent_vectors(u))
+    ee = [0 for _ in 1:r]
+    for (i, k) in enumerate(e)
+      ee[indices[i]] = k
+    end
+    push_term!(ctx, kk(c), ee)
+  end
+  return finish(ctx)
+end
+
+function _evaluate_plain(F::MPolyAnyMap{<:MPolyRing, <:MPolyRing}, u)
+  isdefined(F, :variable_indices) && return _build_poly(u, F.variable_indices, codomain(F))
+  return evaluate(u, F.img_gens)
+end
+
 function _evaluate_plain(F::MPolyAnyMap{<: MPolyRing}, u)
   return evaluate(u, F.img_gens)
 end
 
 # See the comment in MPolyQuo.jl
 function _evaluate_plain(F::MPolyAnyMap{<:MPolyRing, <:MPolyQuoRing}, u)
   A = codomain(F)
+  R = base_ring(A)
+  isdefined(F, :variable_indices) && return A(_build_poly(u, F.variable_indices, R))
   v = evaluate(lift(u), lift.(_images(F)))
   return simplify(A(v))
 end
 
+# The following assumes `p` to be in `S[x₁,…,xₙ]` where `S` is the 
+# actual codomain of the map.
+function _evaluate_with_build_ctx(
+    p::MPolyRingElem, ind::Vector{Int}, cod_ring::MPolyRing
+  )
+  @assert cod_ring === coefficient_ring(parent(p))
+  r = ngens(cod_ring)
+  kk = coefficient_ring(cod_ring)
+  ctx = MPolyBuildCtx(cod_ring)
+  for (q, e) in zip(AbstractAlgebra.coefficients(p), AbstractAlgebra.exponent_vectors(p))
+    ee = [0 for _ in 1:r]
+    for (i, k) in enumerate(e)
+      ee[ind[i]] = k
+    end
+    for (c, d) in zip(AbstractAlgebra.coefficients(q), AbstractAlgebra.exponent_vectors(q))
+      push_term!(ctx, kk(c), ee+d)
+    end
+  end
+  return finish(ctx)
+end
+
+function _evaluate_general(F::MPolyAnyMap{<:MPolyRing, <:MPolyRing}, u)
+  if domain(F) === codomain(F) && coefficient_map(F) === nothing
+    return evaluate(map_coefficients(coefficient_map(F), u,
+                                     parent = domain(F)), F.img_gens)
+  else
+    S = temp_ring(F)
+    if S !== nothing
+      if !isdefined(F, :variable_indices) || coefficient_ring(S) !== codomain(F)
+        return evaluate(map_coefficients(coefficient_map(F), u,
+                                         parent = S), F.img_gens)
+      else
+        tmp_poly = map_coefficients(coefficient_map(F), u, parent = S)
+        return _evaluate_with_build_ctx(tmp_poly, F.variable_indices, codomain(F))
+      end
+    else
+      if !isdefined(F, :variable_indices)
+        return evaluate(map_coefficients(coefficient_map(F), u), F.img_gens)
+      else
+        # For the case where we can recycle the method above, do so.
+        tmp_poly = map_coefficients(coefficient_map(F), u)
+        coefficient_ring(parent(tmp_poly)) === codomain(F) && return _evaluate_with_build_ctx(
+                   tmp_poly,
+                   F.variable_indices,
+                   codomain(F)
+                 )
+        # Otherwise default to the standard evaluation for the time being.
+        return evaluate(tmp_poly, F.img_gens)
+      end
+    end
+  end
+end
+
 function _evaluate_general(F::MPolyAnyMap{<: MPolyRing}, u)
   if domain(F) === codomain(F) && coefficient_map(F) === nothing
     return evaluate(map_coefficients(coefficient_map(F), u,
@@ -151,7 +245,7 @@ function _evaluate_general(F::MPolyAnyMap{<: MPolyRing}, u)
     S = temp_ring(F)
     if S !== nothing
       return evaluate(map_coefficients(coefficient_map(F), u,
-                                parent = S), F.img_gens)
+                                       parent = S), F.img_gens)
     else
       return evaluate(map_coefficients(coefficient_map(F), u), F.img_gens)
     end

src/Rings/mpolyquo-localizations.jl

@@ -620,7 +620,7 @@ function is_zero_divisor(f::MPolyQuoLocRingElem{<:Field})
   # once more functionality is working, it will actually do stuff and
   # the above signature can be widened.
   if is_constant(lifted_numerator(f)) && is_constant(lifted_denominator(f))
-    c = first(coefficients(lift(numerator(f))))
+    c = first(AbstractAlgebra.coefficients(lift(numerator(f))))
     return is_zero_divisor(c)
   end
   return !is_zero(quotient(ideal(parent(f), zero(f)), ideal(parent(f), f)))
@@ -1378,7 +1378,7 @@ end
   gb = groebner_basis(J, ordering=oo)
 
   # TODO: Speed up and use build context.
-  res_gens = elem_type(A)[f for f in gb if all(e -> is_zero(view(e, 1:(n+r))), exponents(f))]
+  res_gens = elem_type(A)[f for f in gb if all(e -> is_zero(view(e, 1:(n+r))), AbstractAlgebra.exponent_vectors(f))]
   img_gens2 = vcat([zero(R) for i in 1:(n+r)], gens(R))
   result = ideal(R, elem_type(R)[evaluate(g, img_gens2) for g in res_gens])
   return result
@@ -2717,3 +2717,90 @@ function dim(R::MPolyQuoLocRing)
   return dim(modulus(R))
 end
 
+### extra methods for speedup of mappings
+# See `src/Rings/MPolyMap/MPolyRing.jl` for the original implementation and 
+# the rationale. These methods make speedup available for quotient rings 
+# and localizations of polynomial rings.
+function _allunique(lst::Vector{T}) where {T<:MPolyLocRingElem}
+  return _allunique(numerator.(lst))
+end
+
+function _allunique(lst::Vector{T}) where {T<:MPolyQuoLocRingElem}
+  return _allunique(lifted_numerator.(lst))
+end
+
+_is_gen(x::MPolyLocRingElem) = is_one(denominator(x)) && _is_gen(numerator(x))
+_is_gen(x::MPolyQuoLocRingElem) = is_one(lifted_denominator(x)) && _is_gen(lifted_numerator(x))
+
+function _evaluate_plain(
+    F::MPolyAnyMap{<:MPolyRing, CT}, u
+  ) where {CT <: Union{<:MPolyLocRing, <:MPolyQuoLocRing}}
+  W = codomain(F)::MPolyLocRing
+  S = base_ring(W)
+  isdefined(F, :variable_indices) && return W(_build_poly(u, F.variable_indices, S))
+  return evaluate(u, F.img_gens)
+end
+
+function _evaluate_general(
+    F::MPolyAnyMap{<:MPolyRing, CT}, u
+  ) where {CT <: Union{<:MPolyQuoRing, <:MPolyLocRing, <:MPolyQuoLocRing}}
+  S = temp_ring(F)
+  if S !== nothing
+    if !isdefined(F, :variable_indices) || coefficient_ring(S) !== codomain(F)
+      return evaluate(map_coefficients(coefficient_map(F), u,
+                                       parent = S), F.img_gens)
+    else
+      tmp_poly = map_coefficients(coefficient_map(F), u, parent = S)
+      return _evaluate_with_build_ctx(tmp_poly, F.variable_indices, codomain(F))
+    end
+  else
+    if !isdefined(F, :variable_indices)
+      return evaluate(map_coefficients(coefficient_map(F), u), F.img_gens)
+    else
+      # For the case where we can recycle the method above, do so.
+      tmp_poly = map_coefficients(coefficient_map(F), u)
+      coefficient_ring(parent(tmp_poly)) === codomain(F) && return _evaluate_with_build_ctx(
+                                                                                            tmp_poly,
+                                                                                            F.variable_indices,
+                                                                                            codomain(F)
+                                                                                           )
+      # Otherwise default to the standard evaluation for the time being.
+      return evaluate(tmp_poly, F.img_gens)
+    end
+  end
+end
+
+function _evaluate_with_build_ctx(
+    p::MPolyRingElem, ind::Vector{Int}, 
+    Q::Union{<:MPolyQuoRing, <:MPolyLocRing, <:MPolyQuoLocRing}
+  )
+  @assert Q === coefficient_ring(parent(p))
+  cod_ring = base_ring(Q)
+  r = ngens(cod_ring)
+  kk = coefficient_ring(cod_ring)
+  ctx = MPolyBuildCtx(cod_ring)
+  for (q, e) in zip(AbstractAlgebra.coefficients(p), AbstractAlgebra.exponent_vectors(p))
+    ee = [0 for _ in 1:r]
+    for (i, k) in enumerate(e)
+      ee[ind[i]] = k
+    end
+    for (c, d) in zip(_coefficients(q), _exponents(q))
+      push_term!(ctx, kk(c), ee+d)
+    end
+  end
+  return Q(finish(ctx))
+end
+
+# The following methods are only safe, because they are called 
+# exclusively in a setup where variables map to variables 
+# and no denominators are introduced. 
+_coefficients(x::MPolyRingElem) = AbstractAlgebra.coefficients(x)
+_coefficients(x::MPolyQuoRingElem) = AbstractAlgebra.coefficients(lift(x))
+_coefficients(x::MPolyLocRingElem) = AbstractAlgebra.coefficients(numerator(x))
+_coefficients(x::MPolyQuoLocRingElem) = AbstractAlgebra.coefficients(lifted_numerator(x))
+
+_exponents(x::MPolyRingElem) = AbstractAlgebra.exponent_vectors(x)
+_exponents(x::MPolyQuoRingElem) = AbstractAlgebra.exponent_vectors(lift(x))
+_exponents(x::MPolyLocRingElem) = AbstractAlgebra.exponent_vectors(numerator(x))
+_exponents(x::MPolyQuoLocRingElem) = AbstractAlgebra.exponent_vectors(lifted_numerator(x))
+