diff --git a/docs/src/AlgebraicGeometry/Schemes/Cycles.md b/docs/src/AlgebraicGeometry/Schemes/Cycles.md
index 88d4efd83c5..192aa1b213b 100644
--- a/docs/src/AlgebraicGeometry/Schemes/Cycles.md
+++ b/docs/src/AlgebraicGeometry/Schemes/Cycles.md
@@ -16,7 +16,7 @@ algebraic_cycle(I::AbsIdealSheaf)
 ```
 ### Properties
 ```@docs
-scheme(D::AbsAlgebraicCycle)
+ambient_scheme(D::AbsAlgebraicCycle)
 components(D::AbsAlgebraicCycle)
 dim(D::AbsAlgebraicCycle)
 irreducible_decomposition(D::AbsAlgebraicCycle)
@@ -46,7 +46,6 @@ weil_divisor(I::AbsIdealSheaf, R::Ring; check::Bool=true)
 Besides the methods for [`AbsAlgebraicCycle`](@ref)
 the following are available.
 ```@docs
-colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
 is_in_linear_system(f::VarietyFunctionFieldElem, D::AbsWeilDivisor; regular_on_complement::Bool=false, check::Bool=true)
 order_of_vanishing(f::VarietyFunctionFieldElem, D::AbsWeilDivisor; check::Bool=true)
 intersect(D::AbsWeilDivisor, E::AbsWeilDivisor; covering::Covering=default_covering(scheme(D)))
@@ -76,8 +75,8 @@ cartier_divisor(IP::AbsProjectiveScheme, f::Union{MPolyDecRingElem, MPolyQuoRing
 ### Attributes
 ```@docs
 ideal_sheaf(C::EffectiveCartierDivisor)
-scheme(C::EffectiveCartierDivisor)
-scheme(C::CartierDivisor)
+ambient_scheme(C::EffectiveCartierDivisor)
+ambient_scheme(C::CartierDivisor)
 coefficient_ring(C::CartierDivisor)
 components(C::CartierDivisor)
 trivializing_covering(C::EffectiveCartierDivisor)
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl b/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
index 3ff9ee18238..5239e0a16fa 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
@@ -23,11 +23,11 @@ coefficient_ring_elem_type(::Type{AbsAlgebraicCycle{S, U}}) where {S, U} = elem_
 ### essential getters and functionality
 
 @doc raw"""
-    scheme(D::AbsAlgebraicCycle)
+    ambient_scheme(D::AbsAlgebraicCycle)
 
 Return the `CoveredScheme` ``X`` on which `D` is defined.
 """
-scheme(D::AbsAlgebraicCycle) = scheme(underlying_cycle(D))
+ambient_scheme(D::AbsAlgebraicCycle) = ambient_scheme(underlying_cycle(D))
 
 # For an element `I` of `components(D)`, this returns the coefficient 
 # of `I` in the formal sum for `D`.
@@ -139,7 +139,7 @@ function underlying_cycle(D::AbsAlgebraicCycle)
 end
 
 ### implementation of the essential functionality
-scheme(D::AlgebraicCycle) = D.X
+ambient_scheme(D::AlgebraicCycle) = D.X
 getindex(D::AlgebraicCycle, I::AbsIdealSheaf) = (D.coefficients)[I]
 
 components(D::AlgebraicCycle) = collect(keys(D.coefficients))
@@ -168,7 +168,7 @@ function AlgebraicCycle(X::AbsCoveredScheme, R::Ring; check::Bool=true)
 end
 
 function zero(D::AbsAlgebraicCycle; check::Bool=true) 
-  return AlgebraicCycle(scheme(D), coefficient_ring(D); check)
+  return AlgebraicCycle(ambient_scheme(D), coefficient_ring(D); check)
 end
 
 # provide non-camelcase methods
@@ -284,7 +284,7 @@ function copy(D::AlgebraicCycle)
   for I in keys(coefficient_dict(D))
     new_dict[I] = D[I]
   end
-  return AlgebraicCycle(scheme(D), coefficient_ring(D), new_dict)
+  return AlgebraicCycle(ambient_scheme(D), coefficient_ring(D), new_dict)
 end
 
 ###############################################################################
@@ -301,7 +301,7 @@ end
 # needs to take care about some left offsets.
 function Base.show(io::IO, ::MIME"text/plain", D::AlgebraicCycle)
   io = pretty(io)
-  X = scheme(D)
+  X = ambient_scheme(D)
   # If the IO context knows about a covering to be used, we use this one.
   # Otherwise, we check whether X has a simplified covering. If not, we use the
   # default covering of X
@@ -348,7 +348,7 @@ end
 
 function Base.show(io::IO, D::AlgebraicCycle)
   io = pretty(io)
-  X = scheme(D)
+  X = ambient_scheme(D)
   eff = all(i >= 0 for i in collect(values(D.coefficients)))
   if length(components(D)) == 1
     prim = D[components(D)[1]] == 1 ? true : false
@@ -360,15 +360,15 @@ function Base.show(io::IO, D::AlgebraicCycle)
   elseif is_terse(io)
     print(io, "Algebraic cycle")
   elseif length(components(D)) == 0
-    print(io, "Zero algebraic cycle on ", Lowercase(), scheme(D))
+    print(io, "Zero algebraic cycle on ", Lowercase(), ambient_scheme(D))
   elseif eff
     if prim
-      print(io, "Irreducible algebraic cycle on ", Lowercase(), scheme(D))
+      print(io, "Irreducible algebraic cycle on ", Lowercase(), ambient_scheme(D))
     else
-      print(io, "Effective algebraic cycle on ", Lowercase(), scheme(D))
+      print(io, "Effective algebraic cycle on ", Lowercase(), ambient_scheme(D))
     end
   else
-    print(io, "Algebraic cycle on ", Lowercase(), scheme(D))
+    print(io, "Algebraic cycle on ", Lowercase(), ambient_scheme(D))
   end
 end
 
@@ -391,8 +391,8 @@ end
 # so the implementation can not be truly generic. 
 
 function +(D::T, E::T) where {T<:AbsAlgebraicCycle}
-  X = scheme(D)
-  X === scheme(E) || error("divisors do not live on the same scheme")
+  X = ambient_scheme(D)
+  X === ambient_scheme(E) || error("divisors do not live on the same scheme")
   R = coefficient_ring(D)
   R === coefficient_ring(E) || error("coefficient rings do not coincide")
   dict = IdDict{AbsIdealSheaf, elem_type(R)}()
@@ -419,7 +419,7 @@ function -(D::T) where {T<:AbsAlgebraicCycle}
   for I in keys(coefficient_dict(D))
     dict[I] = -D[I]
   end
-  return AlgebraicCycle(scheme(D), coefficient_ring(D), dict, check=false)
+  return AlgebraicCycle(ambient_scheme(D), coefficient_ring(D), dict, check=false)
 end
 
 -(D::T, E::T) where {T<:AbsAlgebraicCycle} = D + (-E)
@@ -434,7 +434,7 @@ function *(a::RingElem, E::AbsAlgebraicCycle)
       dict[I] = c
     end
   end
-  return AlgebraicCycle(scheme(E), coefficient_ring(E), dict, check=false)
+  return AlgebraicCycle(ambient_scheme(E), coefficient_ring(E), dict, check=false)
 end
 
 *(a::Int, E::AbsAlgebraicCycle) = coefficient_ring(E)(a)*E
@@ -458,7 +458,7 @@ function irreducible_decomposition(D::AbsAlgebraicCycle)
       k = _colength_in_localization(I, P)
       next_dict[P] = coefficient_ring(D)(k)
     end
-    result = result + a * AlgebraicCycle(scheme(D), coefficient_ring(D), next_dict, check=false)
+    result = result + a * AlgebraicCycle(ambient_scheme(D), coefficient_ring(D), next_dict, check=false)
   end
   return result
 end
@@ -519,7 +519,7 @@ the lengths of all the components of dimension `0` of ``W``.
 """
 function integral(W::AbsAlgebraicCycle; check::Bool=true)
   result = zero(coefficient_ring(W))
-  X = scheme(W)
+  X = ambient_scheme(W)
   for I in components(W)
     @check begin
       dim(I) == 0 || continue
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl b/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
index 5a58eb8f967..254c9bcc1f8 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
@@ -2,6 +2,8 @@ function (C::EffectiveCartierDivisor)(U::AbsAffineScheme)
   return gens(C.I(U))
 end
 
+iszero(C::EffectiveCartierDivisor) = isone(ideal_sheaf(C))
+
 @doc raw"""
     ideal_sheaf(C::EffectiveCartierDivisor)
     
@@ -10,11 +12,11 @@ Return the sheaf of ideals $\mathcal{I}_C \subseteq \mathcal{O}_X$ representing
 ideal_sheaf(C::EffectiveCartierDivisor) = C.I
 
 @doc raw"""
-    scheme(C::EffectiveCartierDivisor)
+    ambient_scheme(C::EffectiveCartierDivisor)
     
 Return the ambient scheme containing `C`. 
 """
-scheme(C::EffectiveCartierDivisor) = C.X
+ambient_scheme(C::EffectiveCartierDivisor) = C.X
 
 @doc raw"""
     trivializing_covering(C::EffectiveCartierDivisor)
@@ -40,11 +42,11 @@ function EffectiveCartierDivisor(I::AbsIdealSheaf;
 end
 
 @doc raw"""
-    scheme(C::CartierDivisor)
+    ambient_scheme(C::CartierDivisor)
     
 Return the ambient scheme containing `C`. 
 """
-scheme(C::CartierDivisor) = C.X
+ambient_scheme(C::CartierDivisor) = C.X
 
 @doc raw"""
     coefficient_ring(C::CartierDivisor)
@@ -64,7 +66,7 @@ Return a list of effective Cartier divisors $C_i$ such that $C$ is a linear comb
 components(C::CartierDivisor) = collect(keys(coefficient_dict(C)))
 
 function +(C::CartierDivisor, D::CartierDivisor) 
-  scheme(C) === scheme(D) || error("divisors must be defined over the same scheme")
+  ambient_scheme(C) === ambient_scheme(D) || error("divisors must be defined over the same scheme")
   coefficient_ring(C) === coefficient_ring(D) || error("divisors must have the same coefficient rings")
   R = coefficient_ring(C)
   coeff_dict = IdDict{EffectiveCartierDivisor, elem_type(R)}()
@@ -83,7 +85,7 @@ function +(C::CartierDivisor, D::CartierDivisor)
       coeff_dict[k] = D[k]
     end
   end
-  return CartierDivisor(scheme(C), coefficient_ring(C), coeff_dict)
+  return CartierDivisor(ambient_scheme(C), coefficient_ring(C), coeff_dict)
 end
 
 function +(C::CartierDivisor, D::EffectiveCartierDivisor) 
@@ -98,6 +100,8 @@ function +(C::EffectiveCartierDivisor, D::CartierDivisor)
   return CartierDivisor(C) + D
 end
 
+zero(D::CartierDivisor) = CartierDivisor(ambient_scheme(D),coefficient_ring(D))
+
 function *(a::RingElem, C::CartierDivisor)
   parent(a) === coefficient_ring(C) || return coefficient_ring(C)(a)*C
   coeff_dict = IdDict{EffectiveCartierDivisor, typeof(a)}()
@@ -109,7 +113,7 @@ function *(a::RingElem, C::CartierDivisor)
       coeff_dict[k] = c
     end
   end
-  return CartierDivisor(scheme(C), coefficient_ring(C), coeff_dict)
+  return CartierDivisor(ambient_scheme(C), coefficient_ring(C), coeff_dict)
 end
 
 function *(a::Integer, C::CartierDivisor)
@@ -121,7 +125,21 @@ function -(C::CartierDivisor, D::CartierDivisor)
 end
 
 function iszero(C::CartierDivisor)
-  return iszero(length(keys(coefficient_dict(C)))) || all(k->iszero(C[k]), components(C))
+  iszero(length(keys(coefficient_dict(C)))) && return true
+  all(iszero, values(coefficient_dict(C))) && return true
+  all(iszero, keys(coefficient_dict(C))) && return true
+  
+  # write C = P - M with P and M effective.
+  # TODO: Multiplying everything together is quick and dirty
+  P = sum(ai*Ci for (Ci,ai) in coefficient_dict(C) if ai>0; init=zero(C))
+  M = sum(-ai*Ci for (Ci,ai) in coefficient_dict(C) if ai<0; init=zero(C))
+  if iszero(length(coefficient_dict(P))) || iszero(length(coefficient_dict(M)))
+      # we know that there is at least one non-zero summand and no cancellation
+      return false
+  end
+  IP = prod(ideal_sheaf(Ci)^ai for (Ci,ai) in coefficient_dict(P))
+  IM = prod(ideal_sheaf(Ci)^ai for (Ci,ai) in coefficient_dict(M))
+  return IP==IM
 end
 
 @doc raw"""
@@ -167,7 +185,7 @@ defined by the formal sum of
 cartier_divisor(E::EffectiveCartierDivisor) = CartierDivisor(E)
 
 function CartierDivisor(C::EffectiveCartierDivisor)
-  return CartierDivisor(scheme(C), ZZ, IdDict([C => one(ZZ)]))
+  return CartierDivisor(ambient_scheme(C), ZZ, IdDict([C => one(ZZ)]))
 end
 
 function CartierDivisor(X::AbsCoveredScheme, kk::Ring)
@@ -175,21 +193,14 @@ function CartierDivisor(X::AbsCoveredScheme, kk::Ring)
 end
 
 function *(a::RingElem, C::EffectiveCartierDivisor)
-  return CartierDivisor(scheme(C), parent(a), IdDict{EffectiveCartierDivisor, typeof(a)}([C => a]))
+  return CartierDivisor(ambient_scheme(C), parent(a), IdDict{EffectiveCartierDivisor, typeof(a)}([C => a]))
 end
 function *(a::Integer, C::EffectiveCartierDivisor)
-  return CartierDivisor(scheme(C), ZZ, IdDict{EffectiveCartierDivisor, elem_type(ZZ)}([C => ZZ(a)]))
+  return CartierDivisor(ambient_scheme(C), ZZ, IdDict{EffectiveCartierDivisor, elem_type(ZZ)}([C => ZZ(a)]))
 end
 
 function ==(C::CartierDivisor, D::CartierDivisor)
-  C === D && return true
-  for k in components(C)
-    iszero(C[k]) || (haskey(coefficient_dict(D), k) && D[k] == C[k]) || error("equality check not implemented in this complicated case")
-  end
-  for k in components(D) 
-    iszero(D[k]) || (haskey(coefficient_dict(C), k) && D[k] == C[k]) || error("equality check not implemented in this complicated case")
-  end
-  return true
+  return iszero(C-D)
 end
 
 @doc raw"""
@@ -262,7 +273,7 @@ end
 Return a `Vector` of pairs ``(I,k)`` corresponding to the irreducible components of ``C``. More precisely,  each ``I`` is a prime  `AbsIdealSheaf` corresponding to an irreducible component of ``C`` and ``k``is the multiplicity of this component in ``C``.
 """
 function irreducible_decomposition(C::EffectiveCartierDivisor)
-  X = scheme(C)
+  X = ambient_scheme(C)
   cov = default_covering(X)
   OOX = OO(X)
 
@@ -301,7 +312,7 @@ function weil_divisor(C::EffectiveCartierDivisor;
 end
 
 function weil_divisor(C::CartierDivisor)
-  X = scheme(C)
+  X = ambient_scheme(C)
   kk = coefficient_ring(C)
   result = WeilDivisor(X, kk)
   for c in components(C)
@@ -310,11 +321,16 @@ function weil_divisor(C::CartierDivisor)
   return result
 end
 
+@doc raw"""
+    intersect(W::WeilDivisor, C::EffectiveCartierDivisor; check::Bool=true)
+    
+Computes the intersection of ``W`` and ``C`` as in [Ful98](@cite) and 
+returns an `AbsAlgebraicCycle` of codimension ``2``.
+"""
 function intersect(W::WeilDivisor, C::EffectiveCartierDivisor; check::Bool=true)
-  X = scheme(W)
+  X = ambient_scheme(W)
   result = zero(W)
-  for I in components(W)
-    @check is_prime(I) "all components of the first argument must be sheaves of prime ideals"
+  for I in components(irreducible_decomposition(W))
     inc_Y = CoveredClosedEmbedding(X, I, check=false)
     #inc_Y = CoveredClosedEmbedding(X, I, covering=trivializing_covering(C), check=false)
     Y = domain(inc_Y)
@@ -333,9 +349,8 @@ returns an `AbsAlgebraicCycle` of codimension ``2``.
 """
 function intersect(W::AbsWeilDivisor, C::CartierDivisor; check::Bool=true)
   result = zero(W)
-  iC = irreducible_decomposition(C)
-  for c in components(iC)
-    result = result + iC[c] * intersect(W, c, check=check)
+  for c in components(C)
+    result = result + C[c] * intersect(W, c, check=check)
   end
   return result
 end
@@ -356,22 +371,8 @@ function intersect(D::CartierDivisor, C::CartierDivisor)
   return intersect(irreducible_decomposition(weil_divisor(D)), C)
 end
 
-
-function pushforward(inc::CoveredClosedEmbedding, W::WeilDivisor)
-  X = domain(inc)
-  Y = codomain(inc)
-  X === scheme(W) || error("divisor not defined on the domain")
-  kk = coefficient_ring(W)
-  ideal_dict = IdDict{AbsIdealSheaf, elem_type(kk)}()
-  for I in components(W)
-    pfI = pushforward(inc)(I)
-    ideal_dict[pfI] = W[I]
-  end
-  return WeilDivisor(Y, kk, ideal_dict, check=false)
-end
-
-dim(C::EffectiveCartierDivisor) = dim(scheme(C))-1
-dim(C::CartierDivisor) = dim(scheme(C))-1
+dim(C::EffectiveCartierDivisor) = dim(ambient_scheme(C))-1
+dim(C::CartierDivisor) = dim(ambient_scheme(C))-1
 
 ###########################################################################
 ## show functions for Cartier divisors
@@ -379,7 +380,7 @@ dim(C::CartierDivisor) = dim(scheme(C))-1
 function Base.show(io::IO, C::EffectiveCartierDivisor)
   io = pretty(io)
   if get(io, :show_semi_compact, false)
-    cov = Oscar._covering_for_printing(io, scheme(C))
+    cov = Oscar._covering_for_printing(io, ambient_scheme(C))
     n = get(io, :label, "")
     _show_semi_compact(io, C, cov, n)
   elseif is_terse(io)
@@ -388,7 +389,7 @@ function Base.show(io::IO, C::EffectiveCartierDivisor)
     print(io, get_attribute(C, :name))
   else
     print(io, "Effective cartier divisor on ", Lowercase())
-    show(io, scheme(C))
+    show(io, ambient_scheme(C))
   end
 end
 
@@ -397,7 +398,7 @@ end
 function Base.show(io::IO, ::MIME"text/plain", C::EffectiveCartierDivisor)
   io = pretty(io)
   I = ideal_sheaf(C)
-  X = scheme(C)
+  X = ambient_scheme(C)
   cov = Oscar._covering_for_printing(io, X)
 
   print(io, "Effective cartier divisor")
@@ -424,7 +425,7 @@ end
 # We usually use "a" for the domain and "b" for the codomain
 function _show_semi_compact(io::IO, C::EffectiveCartierDivisor, cov::Covering, n::String)
   io = pretty(io)
-  X = scheme(C)
+  X = ambient_scheme(C)
   print(io, "Effective cartier divisor")
   if has_attribute(C, :name)
     print(io, " ", get_attribute(C, :name))
@@ -438,7 +439,7 @@ end
 function Base.show(io::IO, C::CartierDivisor)
   io = pretty(io)
   if get(io, :show_semi_compact, false)
-    cov = Oscar._covering_for_printing(io, scheme(C))
+    cov = Oscar._covering_for_printing(io, ambient_scheme(C))
     n = get(io, :label, "")
     _show_semi_compact(io, C, cov, n)
   elseif is_terse(io)
@@ -447,7 +448,7 @@ function Base.show(io::IO, C::CartierDivisor)
     print(io, get_attribute(C, :name))
   else
     print(io, "Cartier divisor on ", Lowercase())
-    show(io, scheme(C))
+    show(io, ambient_scheme(C))
   end
 end
 
@@ -458,13 +459,13 @@ end
 # right.
 function Base.show(io::IO, ::MIME"text/plain", C::CartierDivisor)
   io = pretty(io)
-  X = scheme(C)
+  X = ambient_scheme(C)
   cov = Oscar._covering_for_printing(io, X)
   cc = components(C)
   if length(cc) == 0
     print(io, "Zero cartier divisor ")
     print(io, Indent(), "on ", Lowercase())
-    show(IOContext(io, :covering => cov), scheme(C))
+    show(IOContext(io, :covering => cov), ambient_scheme(C))
     print(io, Dedent())
   else
     print(io, "Cartier divisor")
@@ -473,7 +474,7 @@ function Base.show(io::IO, ::MIME"text/plain", C::CartierDivisor)
     end
     println(io)
     print(io, Indent(), "on ", Lowercase())
-    show(IOContext(io, :covering => cov), scheme(C))
+    show(IOContext(io, :covering => cov), ambient_scheme(C))
     println(io)
     println(io, Dedent(), "with coefficients in ", Lowercase(), coefficient_ring(C))
     print(io, "defined by the formal sum of")
@@ -503,7 +504,7 @@ end
 # We usually use "a" for the domain and "b" for the codomain
 function _show_semi_compact(io::IO, C::CartierDivisor, cov::Covering, n::String)
   io = pretty(io)
-  X = scheme(C)
+  X = ambient_scheme(C)
   cc = components(C)
   if length(cc) == 0
     print(io, "Zero cartier divisor")
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/Types.jl b/src/AlgebraicGeometry/Schemes/Divisors/Types.jl
index 3fe37a19da3..ff82ed4b12c 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/Types.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/Types.jl
@@ -73,7 +73,7 @@ ideal sheaves on ``X``.
   end
 
   function WeilDivisor(C::AlgebraicCycle; check::Bool=true)
-    X = scheme(C)
+    X = ambient_scheme(C)
     @check begin
       for D in keys(coefficient_dict(C))
         is_equidimensional(D) || error("components of a divisor must be sheaves of equidimensional ideals")
@@ -99,7 +99,7 @@ generated by rational functions ``f₁,…,fᵣ ∈ K(X)``.
     length(f) == 0 && return new{typeof(D)}(D, Vector{VarietyFunctionFieldElem}())
     KK = parent(f[1])
     all(g -> (parent(g) === KK), f[2:end]) || error("elements must have the same parent")
-    X = scheme(D)
+    X = ambient_scheme(D)
     X === variety(KK) || error("input not compatible")
 
     @check begin
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl b/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
index 80a3a6d2fea..f9b0194f360 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
@@ -24,7 +24,7 @@ function WeilDivisor(X::AbsCoveredScheme, R::Ring)
 end
 
 function zero(W::WeilDivisor)
-  return WeilDivisor(scheme(W), coefficient_ring(W))
+  return WeilDivisor(ambient_scheme(W), coefficient_ring(W))
 end
 
 # provide non-camelcase methods
@@ -109,7 +109,7 @@ function copy(D::AbsWeilDivisor)
   for I in keys(coefficient_dict(D))
     new_dict[I] = D[I]
   end
-  return WeilDivisor(scheme(D), coefficient_ring(D), new_dict, check=false)
+  return WeilDivisor(ambient_scheme(D), coefficient_ring(D), new_dict, check=false)
 end
 
 function irreducible_decomposition(D::AbsWeilDivisor)
@@ -121,7 +121,7 @@ end
 # relevant information for free
 function Base.show(io::IO, D::AbsWeilDivisor)
   io = pretty(io)
-  X = scheme(D)
+  X = ambient_scheme(D)
   if get(io, :show_semi_compact, false)
     cov = Oscar._covering_for_printing(io, X)
     _show_semi_compact(io, D, cov)
@@ -157,7 +157,7 @@ end
 # everything consistently.
 function _show_semi_compact(io::IO, D::AbsWeilDivisor, cov::Covering)
   io = pretty(io)
-  X = scheme(D)
+  X = ambient_scheme(D)
   C = underlying_cycle(D)
   eff = all(i >= 0 for i in collect(values(coefficient_dict(C))))
   prim = eff && get_attribute(D, :is_prime, false)
@@ -181,7 +181,7 @@ end
 # on the right.
 function Base.show(io::IO, ::MIME"text/plain", D::AbsWeilDivisor)
   io = pretty(io)
-  X = scheme(D)
+  X = ambient_scheme(D)
   cov = Oscar._covering_for_printing(io, X)
   C = underlying_cycle(D)
   eff = all(i >= 0 for i in collect(values(coefficient_dict(C))))
@@ -294,7 +294,7 @@ function ==(D::WeilDivisor, E::WeilDivisor)
 end
 
 @doc raw"""
-    intersect(D::AbsWeilDivisor, E::AbsWeilDivisor; covering::Covering=default_covering(scheme(D)))
+    intersect(D::AbsWeilDivisor, E::AbsWeilDivisor; covering::Covering=default_covering(ambient_scheme(D)))
 
 Return the intersection number of the the Weil divisors `D` and `E`
 on a complete smooth surface as defined in [Har77](@cite).
@@ -303,11 +303,11 @@ on a complete smooth surface as defined in [Har77](@cite).
 The optional keyword argument `covering` specifies the covering to be used for the computation.
 """
 function intersect(D::AbsWeilDivisor, E::AbsWeilDivisor;
-    covering::Covering=default_covering(scheme(D))
+    covering::Covering=default_covering(ambient_scheme(D))
   )
-  X = scheme(D)
+  X = ambient_scheme(D)
   @req dim(X) == 2 "intersection of Weil divisors is only implemented for surfaces."
-  X === scheme(E) || error("divisors do not live on the same scheme")
+  X === ambient_scheme(E) || error("divisors do not live on the same scheme")
   R = coefficient_ring(D)
   R === coefficient_ring(E) || error("divisors do not have the same coefficient ring")
   result = zero(R)
@@ -334,115 +334,21 @@ function intersect(D::AbsWeilDivisor, E::AbsWeilDivisor;
   return result
 end
 
-@attr Bool function has_dimension_leq_zero(I::Ideal)
-  is_one(I) && return true
-  return dim(I) <= 0
-end
-
-@attr Bool function has_dimension_leq_zero(I::MPolyLocalizedIdeal)
-  R = base_ring(I)
-  P = base_ring(R)::MPolyRing
-  J = ideal(P, numerator.(gens(I)))
-  has_dimension_leq_zero(J) && return true
-  is_one(I) && return true
-  return dim(I) <= 0
-end
-
-@attr Bool function has_dimension_leq_zero(I::MPolyQuoLocalizedIdeal)
-  R = base_ring(I)
-  P = base_ring(R)::MPolyRing
-  J = ideal(P, lifted_numerator.(gens(I)))
-  has_dimension_leq_zero(J) && return true
-  is_one(I) && return true
-  return dim(I) <= 0
-end
-
 
-function has_dimension_leq_zero(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
-  for U in keys(object_cache(I))
-    has_dimension_leq_zero(I(U)) || return false
+function pushforward(inc::CoveredClosedEmbedding, W::WeilDivisor)
+  X = domain(inc)
+  Y = codomain(inc)
+  X === ambient_scheme(W) || error("divisor not defined on the domain")
+  kk = coefficient_ring(W)
+  ideal_dict = IdDict{AbsIdealSheaf, elem_type(kk)}()
+  for I in components(W)
+    pfI = pushforward(inc)(I)
+    ideal_dict[pfI] = W[I]
   end
-
-  all_patches = patches(covering)
-  for U in all_patches
-    if !has_dimension_leq_zero(cheap_sub_ideal(I, U))
-      has_dimension_leq_zero(I(U)) || return false
-    end
-  end
-  return true
+  return WeilDivisor(Y, kk, ideal_dict, check=false)
 end
 
-function has_dimension_leq_zero(I::SumIdealSheaf; 
-    covering::Covering=default_covering(scheme(I)),
-    use_decomposition_info::Bool=true
-  )
-  J = summands(I)
-
-  common_patches = keys(object_cache(first(J)))
-  for JJ in J[2:end]
-    common_patches = [U for U in common_patches if U in keys(object_cache(JJ))]
-  end
-
-  for U in common_patches
-    has_dimension_leq_zero(I(U)) || return false
-  end
 
-  all_patches = patches(covering)
-  for U in all_patches
-    patch_ok = false
-    # go through the object cache of the summands
-    if U in keys(object_cache(I))
-      i = I(U)
-      if has_decomposition_info(covering) && use_decomposition_info
-        K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
-        has_dimension_leq_zero(K + i) && (patch_ok = true)
-      else
-        has_dimension_leq_zero(i) && (patch_ok = true)
-      end
-      patch_ok && continue
-    end
-
-    # Do the same for the summands
-    for j in summands(I)
-      if U in keys(object_cache(j))
-        j_loc = j(U)
-        if has_decomposition_info(covering) && use_decomposition_info
-          K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
-          has_dimension_leq_zero(K + j_loc) && (patch_ok = true)
-        else
-          has_dimension_leq_zero(j_loc) && (patch_ok = true)
-        end
-        patch_ok && break
-      end
-    end
-    patch_ok && continue
-
-    # repeat with cheap sub-ideals
-    for j in summands(I)
-      j_cheap = cheap_sub_ideal(j, U)
-      if has_decomposition_info(covering) && use_decomposition_info
-        K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
-        has_dimension_leq_zero(K + j_cheap) && (patch_ok = true)
-      else
-        has_dimension_leq_zero(j_cheap) && (patch_ok = true)
-      end
-      patch_ok && break
-    end
-    patch_ok && continue
-
-    if has_decomposition_info(covering) && use_decomposition_info
-      K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
-      if !has_dimension_leq_zero(cheap_sub_ideal(I, U) + K)
-        has_dimension_leq_zero(I(U) + K) || return false
-      end
-    else
-      if !has_dimension_leq_zero(cheap_sub_ideal(I, U))
-        has_dimension_leq_zero(I(U)) || return false
-      end
-    end
-  end
-  return true
-end
 
 """
     _self_intersection(I::AbsIdealSheaf) -> Integer
@@ -495,91 +401,6 @@ function _is_known_to_be_one(I::SumIdealSheaf, U::AbsAffineScheme;
   return false
 end
 
-@doc raw"""
-    colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
-    
-Return the colength of `I`.
-
-
-"""
-function colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
-  X = scheme(I)
-  all_patches = copy(patches(covering))
-  patches_done = AbsAffineScheme[]
-  patches_todo = AbsAffineScheme[]
-  @vprintln :Divisors 2 "checking colength of $(I)"
-  for U in all_patches
-    dec_inf = (has_decomposition_info(covering) ? elem_type(OO(U))[OO(U)(g) for g in decomposition_info(covering)[U]] : elem_type(OO(U))[])
-    if _is_known_to_be_one(I, U; dec_inf)
-      push!(patches_done, U)
-    else
-      push!(patches_todo, U)
-    end
-  end
-
-  result = 0
-  while length(patches_todo) != 0
-    U = pop!(patches_todo)
-
-    # First do a cheaper test whether this chart needs to be looked at
-    J_cheap = cheap_sub_ideal(I, U)
-    if has_decomposition_info(covering)
-      h = decomposition_info(covering)[U]
-      if isone(J_cheap + ideal(OO(U), elem_type(OO(U))[OO(U)(a) for a in h])) # R(a) is not type stable for R::MPolyQuoRing
-        # push!(patches_done, U)
-        continue
-      end
-    end
-
-    J = I(U)
-    if has_decomposition_info(covering)
-      h = decomposition_info(covering)[U]
-      # The elements in h indicate where components must 
-      # be located so that they can not be spotted in other charts.
-      # We iteratively single out these components by adding a sufficiently high 
-      # power of the equation to the ideal.
-      for f in h
-        g = f
-        while !(g in ideal(OO(U), g*f) + J)
-          g = g * g
-        end
-        J = J + ideal(OO(U), g)
-        if isone(J)
-          push!(patches_done, U)
-          continue
-        end
-      end
-    else
-      # To avoid overcounting, throw away all components that 
-      # were already visible in other charts.
-      for V in patches_done
-        if !haskey(gluings(covering), (U, V))
-          continue
-        end
-        G = covering[U, V]
-        (UV, VU) = gluing_domains(G)
-        UV isa PrincipalOpenSubset || error("method is only implemented for simple gluings")
-        f = complement_equation(UV)
-        # Find a sufficiently high power of f such that it throws
-        # away all components away from the horizon, but does not affect
-        # those on the horizon itself.
-        g = f
-        while !(g in ideal(OO(U), g*f) + J)
-          g = g * g
-        end
-        J = J + ideal(OO(U), g)
-        isone(J) && break
-      end
-    end
-    if !isone(J)
-      JJ = leading_ideal(saturated_ideal(J))
-      A, _ = quo(base_ring(JJ), JJ)
-      result = result + ngens(vector_space(coefficient_ring(base_ring(A)), A)[1])
-    end
-    push!(patches_done, U)
-  end
-  return result
-end
 
 @doc raw"""
     is_in_linear_system(f::VarietyFunctionFieldElem, D::WeilDivisor; regular_on_complement::Bool=true, check::Bool=true) -> Bool
@@ -590,7 +411,7 @@ Return whether the rational function `f` is in the linear system ``|D|``, i.e. i
 - `regular_on_complement` -- set to `true` if `f` is regular on the complement of the support of `D`. 
 """
 function is_in_linear_system(f::VarietyFunctionFieldElem, D::AbsWeilDivisor; regular_on_complement::Bool=false, check::Bool=true)
-  X = scheme(D) 
+  X = ambient_scheme(D) 
   X === variety(parent(f)) || error("schemes not compatible")
   C = simplified_covering(X)
   for I in components(D)
@@ -671,7 +492,7 @@ gen(L::LinearSystem, i::Int) = L.f[i]
 
 Return the variety on which `L` is defined.
 """
-variety(L::LinearSystem) = scheme(weil_divisor(L))
+variety(L::LinearSystem) = ambient_scheme(weil_divisor(L))
 # an alias for the user's convenience 
 scheme(L::LinearSystem) = variety(L)
 
diff --git a/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl b/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
index b23cd710743..81f808122a1 100644
--- a/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
+++ b/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
@@ -344,6 +344,116 @@ end
   return maximum(dims)
 end
   
+@attr Bool function has_dimension_leq_zero(I::Ideal)
+  is_one(I) && return true
+  return dim(I) <= 0
+end
+
+@attr Bool function has_dimension_leq_zero(I::MPolyLocalizedIdeal)
+  R = base_ring(I)
+  P = base_ring(R)::MPolyRing
+  J = ideal(P, numerator.(gens(I)))
+  has_dimension_leq_zero(J) && return true
+  is_one(I) && return true
+  return dim(I) <= 0
+end
+
+@attr Bool function has_dimension_leq_zero(I::MPolyQuoLocalizedIdeal)
+  R = base_ring(I)
+  P = base_ring(R)::MPolyRing
+  J = ideal(P, lifted_numerator.(gens(I)))
+  has_dimension_leq_zero(J) && return true
+  is_one(I) && return true
+  return dim(I) <= 0
+end
+
+
+function has_dimension_leq_zero(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
+  for U in keys(object_cache(I))
+    has_dimension_leq_zero(I(U)) || return false
+  end
+
+  all_patches = patches(covering)
+  for U in all_patches
+    if !has_dimension_leq_zero(cheap_sub_ideal(I, U))
+      has_dimension_leq_zero(I(U)) || return false
+    end
+  end
+  return true
+end
+
+function has_dimension_leq_zero(I::SumIdealSheaf; 
+    covering::Covering=default_covering(scheme(I)),
+    use_decomposition_info::Bool=true
+  )
+  J = summands(I)
+
+  common_patches = keys(object_cache(first(J)))
+  for JJ in J[2:end]
+    common_patches = [U for U in common_patches if U in keys(object_cache(JJ))]
+  end
+
+  for U in common_patches
+    has_dimension_leq_zero(I(U)) || return false
+  end
+
+  all_patches = patches(covering)
+  for U in all_patches
+    patch_ok = false
+    # go through the object cache of the summands
+    if U in keys(object_cache(I))
+      i = I(U)
+      if has_decomposition_info(covering) && use_decomposition_info
+        K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
+        has_dimension_leq_zero(K + i) && (patch_ok = true)
+      else
+        has_dimension_leq_zero(i) && (patch_ok = true)
+      end
+      patch_ok && continue
+    end
+
+    # Do the same for the summands
+    for j in summands(I)
+      if U in keys(object_cache(j))
+        j_loc = j(U)
+        if has_decomposition_info(covering) && use_decomposition_info
+          K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
+          has_dimension_leq_zero(K + j_loc) && (patch_ok = true)
+        else
+          has_dimension_leq_zero(j_loc) && (patch_ok = true)
+        end
+        patch_ok && break
+      end
+    end
+    patch_ok && continue
+
+    # repeat with cheap sub-ideals
+    for j in summands(I)
+      j_cheap = cheap_sub_ideal(j, U)
+      if has_decomposition_info(covering) && use_decomposition_info
+        K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
+        has_dimension_leq_zero(K + j_cheap) && (patch_ok = true)
+      else
+        has_dimension_leq_zero(j_cheap) && (patch_ok = true)
+      end
+      patch_ok && break
+    end
+    patch_ok && continue
+
+    if has_decomposition_info(covering) && use_decomposition_info
+      K = ideal(OO(U), OO(U).(decomposition_info(covering)[U]))
+      if !has_dimension_leq_zero(cheap_sub_ideal(I, U) + K)
+        has_dimension_leq_zero(I(U) + K) || return false
+      end
+    else
+      if !has_dimension_leq_zero(cheap_sub_ideal(I, U))
+        has_dimension_leq_zero(I(U)) || return false
+      end
+    end
+  end
+  return true
+end
+  
 @doc raw"""
     extend!(C::Covering, D::Dict{AffineSchemeType, IdealType}) where {AffineSchemeType<:AffineScheme, IdealType<:Ideal}
 
@@ -708,7 +818,7 @@ Return the order of the rational function `f` on the prime divisor `D`.
 function order_of_vanishing(f::VarietyFunctionFieldElem, D::WeilDivisor; check::Bool=true)
   @check is_prime(D) "divisor must be prime for the order of vanishing to be defined"
   P = components(D)[1]
-  return order_of_vanishing(f, P)
+  return order_of_vanishing(f, P; check=false)
 end 
 
 @doc raw"""
@@ -726,14 +836,11 @@ function order_of_vanishing(
   X = space(I)::AbsCoveredScheme
   X == variety(parent(f)) || error("schemes not compatible")
 
-  #order_dict = Dict{AbsAffineScheme, Int}()
-
   # Since X is integral and I is a sheaf of prime ideals,
   # it suffices to find one chart in which I is non-trivial.
 
   # We look for the chart with the least complexity
   V = first(affine_charts(X))
-  #complexity = Vector{Tuple{AbsAffineScheme, Int}}()
   complexity = inf
   for U in keys(Oscar.object_cache(underlying_presheaf(I))) # Those charts on which I is known.
     U in default_covering(X) || continue
@@ -766,27 +873,6 @@ function order_of_vanishing(
   num_mult = _minimal_power_such_that(J, x->(issubset(quotient(x+K, aR), J)))[1]-1
   den_mult = _minimal_power_such_that(J, x->(issubset(quotient(x+K, bR), J)))[1]-1
   return num_mult - den_mult
-#    # Deprecated code computing symbolic powers explicitly:
-#    L, map = localization(OO(U),
-#                          MPolyComplementOfPrimeIdeal(saturated_ideal(I(U)))
-#                         )
-#    L isa Union{MPolyLocRing{<:Any, <:Any, <:Any, <:Any,
-#                                        <:MPolyComplementOfPrimeIdeal},
-#                     MPolyQuoLocRing{<:Any, <:Any, <:Any, <:Any,
-#                                           <:MPolyComplementOfPrimeIdeal}
-#                    } || error("localization was not successful")
-#
-#    floc = f[U]
-#    a = numerator(floc)
-#    b = denominator(floc)
-#    # TODO: cache groebner bases in a reasonable way.
-#    P = L(prime_ideal(inverted_set(L)))
-#    if one(L) in P
-#      continue # the multiplicity is -∞ in this case and does not count
-#    end
-#    upper = _minimal_power_such_that(P, x->!(L(a) in x))[1]-1
-#    lower = _minimal_power_such_that(P, x->!(L(b) in x))[1]-1
-#    order_dict[U] = upper-lower
 end
 
 @doc raw"""
@@ -2145,3 +2231,88 @@ is_one(II::SingularLocusIdealSheaf) = all(is_one(produce_non_radical_ideal_of_si
 
 in_radical(J::AbsIdealSheaf, II::SingularLocusIdealSheaf) = all(all(radical_membership(g, produce_non_radical_ideal_of_singular_locus(II, U)) for g in gens(J(U))) for U in affine_charts(scheme(J)))
 
+@doc raw"""
+    colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
+    
+Return the colength of `I`.
+
+
+"""
+function colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I)))
+  X = scheme(I)
+  all_patches = copy(patches(covering))
+  patches_done = AbsAffineScheme[]
+  patches_todo = AbsAffineScheme[]
+  @vprintln :Divisors 2 "checking colength of $(I)"
+  for U in all_patches
+    dec_inf = (has_decomposition_info(covering) ? elem_type(OO(U))[OO(U)(g) for g in decomposition_info(covering)[U]] : elem_type(OO(U))[])
+    if _is_known_to_be_one(I, U; dec_inf)
+      push!(patches_done, U)
+    else
+      push!(patches_todo, U)
+    end
+  end
+
+  result = 0
+  while length(patches_todo) != 0
+    U = pop!(patches_todo)
+
+    # First do a cheaper test whether this chart needs to be looked at
+    J_cheap = cheap_sub_ideal(I, U)
+    if has_decomposition_info(covering)
+      h = decomposition_info(covering)[U]
+      if isone(J_cheap + ideal(OO(U), elem_type(OO(U))[OO(U)(a) for a in h])) # R(a) is not type stable for R::MPolyQuoRing
+        # push!(patches_done, U)
+        continue
+      end
+    end
+
+    J = I(U)
+    if has_decomposition_info(covering)
+      h = decomposition_info(covering)[U]
+      # The elements in h indicate where components must 
+      # be located so that they can not be spotted in other charts.
+      # We iteratively single out these components by adding a sufficiently high 
+      # power of the equation to the ideal.
+      for f in h
+        g = f
+        while !(g in ideal(OO(U), g*f) + J)
+          g = g * g
+        end
+        J = J + ideal(OO(U), g)
+        if isone(J)
+          push!(patches_done, U)
+          continue
+        end
+      end
+    else
+      # To avoid overcounting, throw away all components that 
+      # were already visible in other charts.
+      for V in patches_done
+        if !haskey(gluings(covering), (U, V))
+          continue
+        end
+        G = covering[U, V]
+        (UV, VU) = gluing_domains(G)
+        UV isa PrincipalOpenSubset || error("method is only implemented for simple gluings")
+        f = complement_equation(UV)
+        # Find a sufficiently high power of f such that it throws
+        # away all components away from the horizon, but does not affect
+        # those on the horizon itself.
+        g = f
+        while !(g in ideal(OO(U), g*f) + J)
+          g = g * g
+        end
+        J = J + ideal(OO(U), g)
+        isone(J) && break
+      end
+    end
+    if !isone(J)
+      JJ = leading_ideal(saturated_ideal(J))
+      A, _ = quo(base_ring(JJ), JJ)
+      result = result + ngens(vector_space(coefficient_ring(base_ring(A)), A)[1])
+    end
+    push!(patches_done, U)
+  end
+  return result
+end
diff --git a/src/deprecations.jl b/src/deprecations.jl
index 943252dd9b8..e48a6cc42ec 100644
--- a/src/deprecations.jl
+++ b/src/deprecations.jl
@@ -140,3 +140,6 @@ Base.@deprecate_binding QQAbElem QQAbFieldElem
 Base.@deprecate_binding FreeAssAlgIdeal FreeAssociativeAlgebraIdeal
 
 Base.@deprecate_binding in_linear_system is_in_linear_system
+@deprecate scheme(W::AbsAlgebraicCycle) ambient_scheme(W)
+@deprecate scheme(W::CartierDivisor) ambient_scheme(W)
+@deprecate scheme(W::EffectiveCartierDivisor) ambient_scheme(W)
diff --git a/test/AlgebraicGeometry/Schemes/CartierDivisor.jl b/test/AlgebraicGeometry/Schemes/CartierDivisor.jl
index ff96a5f79d8..ce3ea584db5 100644
--- a/test/AlgebraicGeometry/Schemes/CartierDivisor.jl
+++ b/test/AlgebraicGeometry/Schemes/CartierDivisor.jl
@@ -51,11 +51,15 @@ end
   X = covered_scheme(IPX)
   h = (x+y+z+w)^3
   u = (x-y+4*w)
-  C = Oscar.cartier_divisor(IPX, h)
-  D = Oscar.cartier_divisor(IPX, u)
+  C = cartier_divisor(IPX, h)
+  D = cartier_divisor(IPX, u)
   @test 2*C == C+C
   @test iszero(D-D)
   @test 3*(C + D) == 3*C + 3*D
+  Z = Oscar.cartier_divisor(IPX, S(1))
+  @test iszero(Z)
+  C2 = cartier_divisor(IPX, h^2)
+  @test iszero(2*C - C2)
 end
 
 @testset "conversion of Cartier to Weil divisors" begin
diff --git a/test/AlgebraicGeometry/Schemes/elliptic_surface_disabled.jl b/test/AlgebraicGeometry/Schemes/elliptic_surface_disabled.jl
new file mode 100644
index 00000000000..31635576438
--- /dev/null
+++ b/test/AlgebraicGeometry/Schemes/elliptic_surface_disabled.jl
@@ -0,0 +1,276 @@
+# The tests for elliptic surfaces have lead to random failure of the CI-tests, 
+# see issue no. 3676. 
+#
+# They are now disabled, also because in general they tend to take quite long.
+# We keep them, however, to allow for running them locally.
+@testset "elliptic surfaces" begin
+  @testset "trivial lattice" begin
+    k = GF(29)
+    # The generic fiber of the elliptic fibration
+    # as an elliptic curve over k(t)
+    kt, t = polynomial_ring(k, :t)
+    kP1 = fraction_field(kt)
+    t = gen(kP1)
+    E = elliptic_curve(kP1, [3*t^8 + 10*t^7 + 6*t^6 + 17*t^5 + 25*t^4 + 4*t^3 + 23*t^2 + 9*t + 14, 5*t^12 + 25*t^11 + 2*t^10 + 28*t^9 + 28*t^8 + 19*t^7 + 3*t^6 + 17*t^5 + 19*t^4 + 12*t^3 + 25*t^2 + 12*t + 6])
+    # A basis for the Mordell-Weil group of E
+    mwl_basis = [E(collect(i)) for i in [
+    (12*t^4 + 21*t^3 + 5*t^2 + 12*t + 18, 23*t^5 + 7*t^4 + 22*t^3 + 13*t^2),
+    (12*t^4 + 20*t^3 + 22*t^2 + 27*t + 18, 15*t^4 + 12*t^3 + 12*t),
+    (12*t^4 + 20*t^3 + 27*t^2 + 11*t + 18, -(16*t^4 + 24*t^3 + 3*t^2 + 24*t)),
+    (4*t^4 + 5*t^3 + 5*t^2 + 6*t + 13,  9*t^6 + 21*t^5 + 17*t^4 + 12*t^2 + 3*t + 6)]]
+    S = elliptic_surface(E, 2, mwl_basis)
+    weierstrass_model(S)
+    weierstrass_contraction(S)
+    trivial_lattice(S)
+
+    E = elliptic_curve(kP1, [0,0,0,1,t^10])
+    X = elliptic_surface(E, 2)
+    triv = trivial_lattice(X)
+    @test det(triv[2])==-3
+  end  
+  
+
+  # This test takes about 1 minute
+  @testset "mordel weil lattices" begin
+    k = GF(29,2)
+    # The generic fiber of the elliptic fibration
+    # as an elliptic curve over k(t)
+    kt, t = polynomial_ring(k, :t)
+    kP1 = fraction_field(kt)
+    t = gen(kP1)
+
+    E = elliptic_curve(kP1, [0, 2*t^4 + 28*t^2, 0, 27*t^6 + 19, 10*t^2])
+    P = E([0,sqrt(k(10))*t,1])
+    X = elliptic_surface(E, 2, [P])
+    triv = trivial_lattice(X)
+    @test det(triv[2]) == 256
+    @test length(Oscar.mordell_weil_torsion(X)) == 1
+    alg = algebraic_lattice(X)
+    @test det(alg[3]) == -192
+    @test det(mordell_weil_lattice(X)) == 3
+    
+    X1 = elliptic_surface(short_weierstrass_model(E)[1],2)
+    Oscar.isomorphism_from_generic_fibers(X,X1)
+  end
+  #=
+  # this test is quite expensive
+  # probably because it is over QQ
+  # ... and because there are loads of complicated singular fibers
+  Qt, t = polynomial_ring(QQ, :t)
+  Qtf = fraction_field(Qt)
+  E = elliptic_curve(Qtf, [0,0,0,0,t^5*(t-1)^2])
+  X3 = elliptic_surface(E, 2)
+  weierstrass_contraction(X3)
+  trivial_lattice(X3)
+  =#
+
+  @testset "elliptic parameter" begin
+    k = GF(29)
+    kt, _ = polynomial_ring(k, :t)
+    kP1 = fraction_field(kt); t = gen(kP1)
+    E = elliptic_curve(kP1, [(21*t^7+6*t^6+11*t^4),(21*t^10+15*t^9+17*t^7+18*t^6+t^5)])
+    mwl_basis = E.([
+    [t^3 + 24*t^2 + 22*t + 5, 10*t^5 + 6*t^4 + 6*t^3 + 8*t^2 + 14*t + 3],
+    [7*t^3 + 24*t^2 + 9, 9*t^5 + 18*t^4 + 12*t^3 + 8*t^2 + 2],
+    [13*t^3 + 9*t + 24, 2*t^5 + 7*t^4 + 6*t^3 + 16*t^2 + 16*t + 22],
+    [(17*t^5 + 14*t^4 + 28*t^2 + 2*t + 1)//(t^2 + 11*t + 23), (28*t^8 + 10*t^7 + 27*t^6 + 28*t^5 + 3*t^4 + 27*t^3 + 3*t + 1)//(t^3 + 2*t^2 + 11*t + 25)],
+    [(11*t^5 + 22*t^4 + 23*t^3 + 14*t^2 + 17*t + 16)//(t^2 + 24*t + 28), (22*t^8 + 6*t^7 + 21*t^6 + 24*t^5 + 4*t^4 + 23*t^3 + 25*t^2 + 15*t + 6)//(t^3 + 7*t^2 + 26*t + 17)]])
+
+    X = elliptic_surface(E, 2, mwl_basis)
+    ff = QQFieldElem[9, 4, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, -1]
+    @test det(algebraic_lattice(X)[3])==-1183
+    @test length(Oscar.mordell_weil_torsion(X)) == 0 # no torsion points
+    # u = elliptic_parameter(X, ff)
+    g, phi = two_neighbor_step(X, ff)
+
+
+    mwl_gens_new = vcat([mwl_basis[1] + mwl_basis[2], mwl_basis[1] - mwl_basis[2]])
+    set_mordell_weil_basis!(X, mwl_gens_new)
+    @test det(algebraic_lattice(X)[3])==96
+    Oscar.algebraic_lattice_primitive_closure!(X)
+    @test det(algebraic_lattice(X)[3])==24
+
+  end
+end
+
+@testset "normalize_quartic and transform_to_weierstrass" begin
+  R, (x, y) = polynomial_ring(QQ, [:x, :y])
+  P, (u, v) = polynomial_ring(QQ, [:u, :v])
+  f = (3*x^2 - 5)^2*(y - 5*x^3 + 30*x - 5)^2 - (7*x^3 + x^2 -5*x + 2)
+  g, trans = Oscar._normalize_hyperelliptic_curve(f, parent=P)
+  @test trans(f) == g
+
+  R, (x, y) = polynomial_ring(QQ, [:x, :y])
+  f = y^2 - 4*x^4 + 5*x^3 - 3*x^2 + 7*x - 4
+  f_trans, trans = Oscar.transform_to_weierstrass(f, x, y, QQ.([0, 2]))
+  @test trans(f//1) == (16*x^6 + 1//8*x^5 + 14*x^4*y - 16383//4096*x^4 - 16*x^3*y^2 + 647//128*x^3*y)//y^4
+end
+
+@testset "normalize_quartic and transform_to_weierstrass over function fields" begin
+  pt, t = QQ[:t]
+  kt = fraction_field(pt)
+  R, (x, y) = polynomial_ring(kt, [:x, :y])
+  P, (u, v) = polynomial_ring(kt, [:u, :v])
+  f = (3*x^2 - 5*t^2)^2*(y - 5*t*x^3 + 30*x - 5)^2 - (x^3 + t*x^2 -5*x + 2*t^2)
+  g, trans = Oscar._normalize_hyperelliptic_curve(f, parent=P)
+  @test trans(f) == g
+
+  R, (x, y) = polynomial_ring(kt, [:x, :y])
+  f = y^2 - 4*x^4 + 5*x^3 - 3*x^2 + 7*x - 4*t^8
+  f_trans, trans = Oscar.transform_to_weierstrass(f, x, y, kt.([0, 2*t^4]))
+  @test f_trans == x^3 + (-3//8*t^8 + 49//128)//t^16*x^2 + 7//8//t^8*x*y + (-1//4*t^24 + 9//256*t^16 - 147//2048*t^8 + 2401//65536)//t^32*x - y^2 + (5//16*t^16 - 21//128*t^8 + 343//2048)//t^24*y
+  
+  R, (x, y) = polynomial_ring(kt, [:x, :y])
+  f = y^2 - 4*x^4 + 5*x^3 - 3*x^2 + 7*x - 4*t^8
+  Y,trafo = elliptic_surface(f, kt.([0, 2*t^4]))
+end
+
+
+
+# These tests are disabled, because they take too long, about 5 minutes. But one can run them if in doubt.
+# most of the time is spent in decomposing the fibers
+@testset "two neighbor steps" begin
+  K = GF(7)
+  Kt, t = polynomial_ring(K, :t)
+  Ktf = fraction_field(Kt)
+  E = elliptic_curve(Ktf, [0, -t^3, 0, t^3, 0])
+  P = E([t^3, t^3])
+  X2 = elliptic_surface(E, 2, [P]);
+  KX2 = function_field(X2; check=false)
+  U = weierstrass_chart(X2)
+  (xx, yy, tt) = ambient_coordinates(U)
+  u4 = KX2(yy, xx*tt)
+  u5 = KX2(yy + tt^3, xx*tt - tt^4)
+  u6 = KX2(yy + tt^3, xx - tt^3)
+  fibers_in_X2 = [QQ.(vec(collect(i))) for i in [
+                                                 [4   2   0   0   0   0   0   0   0   -4   -4   -8   -7   -6   -5   -4   -3   -2   -1   0],
+                                                 [1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0],
+                                                 [5   2   -2   -3   -4   -3   -2   -1   -2   -5   -4   -8   -7   -6   -5   -4   -3   -2   -1   0],
+                                                 [4   2   -2   -4   -6   -9//2   -3   -3//2   -7//2   -3   -5//2   -5   -9//2   -4   -7//2   -3   -5//2   -2   -3//2   0],
+                                                 [2   1   -1   -2   -3   -2   -1   0   -2   -1   -1   -2   -2   -2   -2   -2   -2   -1   0   1],
+                                                 [2   1   0   0   0   0   0   0   0   -2   -2   -4   -4   -4   -4   -3   -2   -1   0   1]
+                                                ]]
+  g4,_ = two_neighbor_step(X2, fibers_in_X2[4])  # this should be the 2-torsion case
+  g5,_ = two_neighbor_step(X2, fibers_in_X2[5])  # the non-torsion case
+  g6,_ = two_neighbor_step(X2, fibers_in_X2[6])  # the non-torsion case
+end
+
+
+
+# The following tests take roughly 10 minutes which is too much for the CI testsuite.
+# We keep them for long running tests to be checked every now and then.
+@testset "translation on elliptic surfaces" begin
+  P, t = GF(113)[:t]
+  kt = fraction_field(P)
+
+  R, (x, y) = kt[:x, :y]
+
+  f = y^2 - (x^3 + (37*t^8 + 56*t^7 + 102*t^6 + 97*t^5 + 21*t^4 + 14*t^3 + 24*t^2 + 23*t + 59)*x + 18*t^12 + 111*t^11 + 17*t^10 + 25*t^9 + 108*t^8 + 91*t^7 + 90*t^6 + 21*t^5 + 68*t^4 + 8*t^3 + 92*t^2 + 66*t + 31)
+
+  E = elliptic_curve(f, x, y)
+
+  P = E([(10*t^6 + 23*t^5 + 94*t^4 + 32*t^3 + t^2 + 40*t + 52)//(t^2 + 77*t + 98), (22*t^9 + 40*t^8 + 63*t^7 + 42*t^6 + 34*t^5 + 48*t^4 + 72*t^3 + 92*t^2 + 85*t + 91)//(t^3 + 59*t^2 + 68*t + 44)])
+  
+  pts = E.([
+  [111*t^4 + 94*t^3 + 57*t^2 + 55*t + 31, 106*t^6 + 2*t^5 + 61*t^4 + 92*t^3 + 105*t^2 + 42*t + 24],   
+  [26*t^3 + 32*t^2 + 41*t + 95, 73*t^6 + 101*t^5 + 83*t^4 + 45*t^3 + 53*t^2 + 97], 
+  [23*t^4 + 35*t^3 + 27*t^2 + 106*t + 40, 43*t^6 + 103*t^5 + 15*t^4 + 44*t^3 + 34*t^2 + 30*t + 88],
+  [109*t^4 + 49*t^3 + 82*t^2 + 69*t + 80, 22*t^6 + 51*t^5 + 27*t^4 + 5*t^3 + 27*t^2 + 80*t + 62],
+  [103*t^4 + 27*t^3 + 44*t^2 + 111*t + 1, 2*t^6 + 72*t^5 + 37*t^4 + 59*t^3 + 55*t^2 + 106*t + 59],
+  [6*t^4 + 48*t^3 + 70*t^2 + 112*t + 61, 111*t^6 + 41*t^5 + 76*t^4 + 54*t^3 + 58*t^2 + 7*t + 54],
+  [42*t^4 + 7*t^3 + 106*t^2 + 112*t + 69, 96*t^6 + 100*t^5 + 48*t^4 + 34*t^3 + 112*t^2 + 83*t + 74],
+  [59*t^4 + 55*t^3 + 50*t^2 + 36*t + 8, 15*t^6 + 23*t^5 + 92*t^4 + 64*t^3 + 103*t^2 + 17*t + 87]
+  ])
+
+  X = elliptic_surface(E, 2, pts)
+
+  #set_verbose_level(:EllipticSurface, 5)
+
+  # The following should not take more than at most two minutes.
+  # But it broke the tests at some point leading to timeout, 
+  # so we put it here to indicate regression.
+  for (i, g) in enumerate(values(gluings(default_covering(X))))
+    gluing_domains(g) # Trigger the computation once
+  end
+
+  D_P = section(X, P)
+ 
+  II = first(components(D_P));
+ 
+  trans = Oscar.translation_morphism(X, P; divisor=D_P)
+ 
+  JJ = Oscar._pushforward_section(trans, P; divisor=D_P)
+ 
+  IIX = first(components(section(X, 2*P)));
+  # We have little chance to get through with the computations on all charts. 
+  # But it suffices to compare the result on the weierstrass charts.
+  weier = weierstrass_chart_on_minimal_model(X)
+  @test IIX(weier) == JJ(weier)
+
+  
+  # pushforward of the whole algebraic lattice with verification of the result.
+  lat, _, A = algebraic_lattice(X)
+  A = gram_matrix(ambient_space(A))
+ 
+  ll = Oscar._pushforward_lattice_along_isomorphism(trans)
+  res_mat = [Oscar.basis_representation(X, d) for d in ll]
+  
+  res_mat = matrix(QQ, res_mat)
+  # Check that the base change matrix is indeed orthogonal for the given lattice
+  @test res_mat*A*transpose(res_mat) == A
+end
+
+
+@testset "moebius transformations" begin
+  k = GF(113)
+  kt,t = polynomial_ring(k,:t)
+  Ft = fraction_field(kt)
+
+  E = elliptic_curve(Ft, Ft.([0,0,0,(56*t^8 + 56),0]));
+  basis_mwl = [[112*t^4 + 112*t^3 + 56*t^2 + 15*t + 1, 100*t^6 + 24*t^5 + 56*t^4 + 13*t^3 + 64*t^2 + 89*t + 31],
+               [31*t^4 + 15*t^2 + 82, 44*t^5 + 14*t^3 + 69*t],
+               [82*t^4 + 13, 37*t^4 + 10],
+               [91*t^4 + 16*t^3 + 25*t^2 + 14*t + 22, 18*t^6 + 55*t^5 + 45*t^4 + 44*t^3 + 110*t^2 + 58*t + 69],
+               [32*t^4 + 72*t^2 + 81, 78*t^6 + 79*t^4 + 58*t^2 + 40],
+               [56*t^4 + 49*t^3 + 85*t^2 + 57*t + 57, 72*t^6 + 25*t^5 + 22*t^4 + 101*t^3 + 104*t^2 + 88*t + 50],
+               [22*t^4 + 87*t^3 + 77*t^2 + 26*t + 22, 44*t^6 + 27*t^5 + 68*t^4 + 98*t^3 + 45*t^2 + 27*t + 69],
+               [112*t^4 + 44*t^3 + 49*t^2 + 44*t + 112, 100*t^6 + 74*t^5 + 49*t^4 + 8*t^3 + 49*t^2 + 74*t + 100]]
+  basis_mwl = [E(i) for i in basis_mwl];  
+
+  X = elliptic_surface(E, 2, basis_mwl[1:0])  # speed up test by hiding the sections
+  
+  moeb = Oscar.admissible_moebius_transformations(X, X)
+
+  @test length(moeb) == 16 # This must really be the number for this particular surface.
+
+  for phi in moeb
+    @test Oscar.check_isomorphism_on_generic_fibers(phi)
+  end
+
+  phi = moeb[5] # Just pick one which is not the identity
+
+  W = weierstrass_chart_on_minimal_model(X)
+
+  for (i, U) in enumerate(affine_charts(X))
+    phi_loc = Oscar.cheap_realization(phi, W, U)
+    @test all(iszero(pullback(phi_loc)(lifted_numerator(g))) for g in gens(modulus(OO(U))))
+  end
+  
+  phi_star = Oscar.pushforward_on_algebraic_lattices(phi)
+  # Check that the pushforward preserves the numerical lattice of X
+  Triv = algebraic_lattice(X)[3]
+  @test phi_star(Triv)==Triv
+  
+  # test translation
+  P = Oscar.mordell_weil_torsion(X)[1]
+  torsion_translation = Oscar.translation_morphism(X, P)
+  tors_of_P = pushforward(torsion_translation,zero_section(X)) 
+  Psect = section(X, P)
+  @test tors_of_P == Psect
+  
+  # add a section so that there is something to compute (although rather trivial)
+  set_mordell_weil_basis!(X, basis_mwl[1:1])
+  P = Oscar.extract_mordell_weil_basis(torsion_translation)
+  @test length(P) == 2
+end
+=#