diff --git a/experimental/Schemes/src/Auxiliary.jl b/experimental/Schemes/src/Auxiliary.jl
index b8cf7d827365..1e245b36460e 100644
--- a/experimental/Schemes/src/Auxiliary.jl
+++ b/experimental/Schemes/src/Auxiliary.jl
@@ -128,12 +128,11 @@ end
 
 function pullback(f::AbsCoveredSchemeMorphism, C::CartierDivisor)
   R = coefficient_ring(C)
-  C = CartierDivisor(domain(f), R)
-  pb = pullback(f)
-  for (c,D) in coefficient_dict(C)
-    C += c*pb(C)
+  result = CartierDivisor(domain(f), R)
+  for (D, c) in coefficient_dict(C)
+    result += c*pullback(f, D)
   end
-  return C
+  return result
 end
 
 function pullback(f::AbsCoveredSchemeMorphism, CC::Covering)
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl b/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
index ea7b74e85290..84cf9e4b7131 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/AlgebraicCycles.jl
@@ -395,7 +395,7 @@ end
 # Note that we need one minimal concrete type for the return values, 
 # so the implementation can not be truly generic. 
 
-function +(D::T, E::T) where {T<:AbsAlgebraicCycle}
+function +(D::AbsAlgebraicCycle, E::AbsAlgebraicCycle)
   X = ambient_scheme(D)
   X === ambient_scheme(E) || error("divisors do not live on the same scheme")
   R = coefficient_ring(D)
@@ -454,9 +454,13 @@ Return a cycle ``E`` equal to ``D`` but as a formal sum ``E = ∑ₖ aₖ ⋅ I
 where the `components` ``Iₖ`` of ``E`` are all sheaves of prime ideals.
 """
 function irreducible_decomposition(D::AbsAlgebraicCycle)
-  all(is_prime, keys(coefficient_dict(D))) && return D
+  @vprint :Divisors 4 "computing irreducible decomposition for $D"
   result = zero(D)
   for (I, a) in coefficient_dict(D)
+    if is_prime(I)
+      result[I] = a
+      continue
+    end
     next_dict = IdDict{AbsIdealSheaf, elem_type(coefficient_ring(D))}()
     decomp = maximal_associated_points(I)
     for P in decomp
@@ -465,7 +469,34 @@ function irreducible_decomposition(D::AbsAlgebraicCycle)
     end
     result = result + a * AlgebraicCycle(ambient_scheme(D), coefficient_ring(D), next_dict, check=false)
   end
-  return result
+  return _unique_prime_components(result)
+end
+
+# Given a cycle `D` with only prime components, compare the components
+# and gather the coefficients of equal ones so that the result has 
+# pairwise distinct components.
+function _unique_prime_components(D::AbsAlgebraicCycle)
+  @hassert :Divisors 2 all(is_prime, keys(coefficient_dict(D))) 
+  buckets = Vector{Vector{AbsIdealSheaf}}()
+  for P in keys(coefficient_dict(D))
+    found = false
+    for bucket in buckets
+      if P == first(bucket)
+        push!(bucket, P)
+        found = true
+        break
+      end
+    end
+    !found && push!(buckets, [P])
+  end
+  R = coefficient_ring(D)
+  coeff_dict = IdDict{AbsIdealSheaf, elem_type(R)}()
+  for bucket in buckets
+    c = sum(D[P] for P in bucket; init=zero(R))
+    is_zero(c) && continue
+    coeff_dict[first(bucket)] = c
+  end
+  return AlgebraicCycle(ambient_scheme(D), coefficient_ring(D), coeff_dict; check=false)
 end
 
 function _colength_in_localization(Q::AbsIdealSheaf, P::AbsIdealSheaf; covering=simplified_covering(scheme(P)))
@@ -534,6 +565,25 @@ function integral(W::AbsAlgebraicCycle; check::Bool=true)
   return result
 end
 
+# Getters for the components as honest algebraic cycles, not ideal sheaves.
+function components(::Type{T}, D::AbsAlgebraicCycle) where {T <: AbsAlgebraicCycle}
+  X = scheme(D)
+  R = coefficient_ring(D)
+  return [AlgebraicCycle(X, R, IdDict{AbsIdealSheaf, elem_type(R)}([I=>one(R)]); check=false)::T for I in components(D)]
+end
+
+function components(::Type{T}, D::AbsWeilDivisor) where {T <: AbsWeilDivisor}
+  X = scheme(D)
+  R = coefficient_ring(D)
+  return [WeilDivisor(X, R, IdDict{AbsIdealSheaf, elem_type(R)}([I=>one(R)]); check=false)::T for I in components(D)]
+end
+
+function getindex(D::AbsAlgebraicCycle, C::AbsAlgebraicCycle)
+  comps = components(C)
+  @assert isone(length(comps))
+  return D[first(comps)]
+end
+
 function Base.hash(X::AbsAlgebraicCycle, u::UInt)
   return u
 end
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl b/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
index 36f7e74ae448..bde1c268e288 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/CartierDivisor.jl
@@ -88,6 +88,10 @@ function +(C::CartierDivisor, D::CartierDivisor)
   return CartierDivisor(ambient_scheme(C), coefficient_ring(C), coeff_dict)
 end
 
+function -(C::CartierDivisor, E::EffectiveCartierDivisor)
+  return C - 1*E
+end
+
 function +(C::CartierDivisor, D::EffectiveCartierDivisor) 
   return C + CartierDivisor(D)
 end
@@ -322,17 +326,16 @@ function weil_divisor(C::CartierDivisor)
 end
 
 @doc raw"""
-    intersect(W::WeilDivisor, C::EffectiveCartierDivisor; check::Bool=true)
+    intersect(W::AbsWeilDivisor, 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)
+function intersect(W::AbsWeilDivisor, C::EffectiveCartierDivisor; check::Bool=true)
   X = ambient_scheme(W)
   result = zero(W)
   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)
     pbC = pullback(inc_Y)(C) # Will complain if the defining equation of C is vanishing identically on Y
     W_sub = weil_divisor(pbC)
@@ -342,7 +345,7 @@ function intersect(W::WeilDivisor, C::EffectiveCartierDivisor; check::Bool=true)
 end
 
 @doc raw"""
-    intersect(W::WeilDivisor, C::CartierDivisor; check::Bool=true)
+    intersect(W::AbsWeilDivisor, C::CartierDivisor; check::Bool=true)
 
 Computes the intersection of ``W`` and ``C`` as in [Ful98](@cite) and 
 returns an `AbsAlgebraicCycle` of codimension ``2``.
diff --git a/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl b/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
index 33a5ce98411e..b8b384f42767 100644
--- a/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
+++ b/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl
@@ -630,3 +630,215 @@ function order_of_vanishing(
   I = components(D)[1]
   return order_of_vanishing(f, I, check=false)
 end
+
+# Compute the colength of I_P in the localization R_P.
+# This assumes that R itself is a domain and that I is of height 1. 
+# Then there exists a regular 
+# point on Spec(R) and hence R_P is also regular and a UFD. 
+# Being of height 1, I_P must then be principal.
+function _colength_in_localization(I::Ideal, P::Ideal)
+  R = base_ring(I)
+  @assert R === base_ring(P)
+  U = MPolyComplementOfPrimeIdeal(saturated_ideal(P); check=false)
+  L, loc = localization(R, U)
+  I_loc = loc(I)
+  @assert base_ring(I_loc) === L
+  P_loc = loc(P)
+  x = _find_principal_generator(P_loc)
+  y = one(L)
+  k = 0
+  while true
+    (y in I_loc) && return k
+    y = y*x
+    k += 1
+  end
+end
+
+# same assumptions as above apply.
+function _find_principal_generator(I::Union{<:MPolyLocalizedIdeal, <:MPolyQuoLocalizedIdeal})
+  L = base_ring(I)
+  g = gens(I)
+  g = sort!(g, by=x->total_degree(lifted_numerator(x)))
+  for x in g
+    is_zero(x) && continue
+    ideal(L, x) == I && return x
+  end
+  error("no principal generator found")
+end
+
+# produce the principal divisor associated to a rational function
+principal_divisor(::Type{WeilDivisor}, f::VarietyFunctionFieldElem; 
+    ring::Ring=ZZ, covering::Covering=default_covering(scheme(f))
+  ) = weil_divisor(f; ring, covering)
+  
+function weil_divisor(
+    f::VarietyFunctionFieldElem; 
+    ring::Ring=ZZ, covering::Covering=default_covering(scheme(f))
+  )
+  @vprint :Divisors 4 "calculating principal divisor for $f\n"
+  X = scheme(f)
+  ideal_dict = IdDict{AbsIdealSheaf, elem_type(ring)}()
+  for U in patches(covering)
+    # TODO: We compute the dimensions again and again. 
+    # Instead we should store these things in some matured version of `decomposition_info`!
+    has_decomposition_info(covering) && (dim(ideal(OO(U), decomposition_info(covering)[U])) <= dim(X) - 2) && continue
+
+    # covering to take a shortcut in the
+    @vprint :Divisors 4 "doing patch $U\n"
+    inc_dict = IdDict{AbsIdealSheaf, elem_type(ring)}()
+    f_loc = f[U]
+    num = numerator(f_loc)
+    den = denominator(f_loc)
+    num_ideal = ideal(OO(U), num)
+    den_ideal = ideal(OO(U), den)
+    num_dec = primary_decomposition(num_ideal)
+    den_dec = primary_decomposition(den_ideal)
+    @vprint :Divisors 4 "  numerator:\n"
+    for (_, P) in num_dec
+      @vprint :Divisors 4 "    $P\n"
+      # If this component was already seen in another patch, skip it.
+      new_comp = PrimeIdealSheafFromChart(X, U, P)
+      @vprint :Divisors 4 "    $(any(new_comp == PP for PP in keys(ideal_dict)) ? "already found" : "new component")\n"
+      any(new_comp == PP for PP in keys(ideal_dict)) && continue 
+      c = _colength_in_localization(num_ideal, P)
+      @vprint :Divisors 4 "    multiplicity $c\n"
+      inc_dict[new_comp] = c
+    end
+    @vprint :Divisors 4 "  denominator:\n"
+    for (_, P) in den_dec
+      # If this component was already seen in another patch, skip it.
+      new_comp = PrimeIdealSheafFromChart(X, U, P)
+      @vprint :Divisors 4 "    $(any(new_comp == PP for PP in keys(ideal_dict)) ? "already found" : "new component")\n"
+      any(new_comp == PP for PP in keys(ideal_dict)) && continue 
+      c = _colength_in_localization(den_ideal, P)
+      @vprint :Divisors 4 "    multiplicity $c\n"
+      key_list = collect(keys(inc_dict))
+      k = findfirst(==(new_comp), key_list)
+      if k === nothing
+        is_zero(c) && continue
+        inc_dict[new_comp] = -c
+      else
+        d = inc_dict[key_list[k]]
+        if c == d
+          delete!(inc_dict, key_list[k])
+          continue
+        end
+        inc_dict[key_list[k]] = d - c
+      end
+    end
+    for (pp, c) in inc_dict
+      ideal_dict[pp] = c
+    end
+  end
+  return WeilDivisor(X, ring, ideal_dict; check=false)
+end
+
+@doc raw"""
+    move_divisor(D::AbsWeilDivisor; check::Bool=false)
+
+Given an `AbsWeilDivisor` `D` on a scheme `X`, create a principal divisor 
+`div(f)` for some rational function, so that `D - div(f)` is supported 
+in (non-closed) scheme-theoretic points which are different from those of `D`.
+
+Keyword arguments:
+  * `randomization`: By default, the choices made to create `f` keep it as simple as possible. However, one might encounter constellations where this will just swap two components of `D`. In order to avoid this, one can then switch on randomization here. 
+  * `is_prime`: Set this to `true` if you know your divisor `D` to already be prime and avoid expensive internal checks.
+"""
+function move_divisor(
+    D::AbsWeilDivisor; 
+    check::Bool=false,
+    randomization::Bool=false,
+    is_prime::Bool=false
+  )
+  X = scheme(D)
+  @check is_irreducible(X) && is_reduced(X) "scheme must be irreducible and reduced"
+  is_zero(D) && return D
+
+  if !is_prime && !Oscar.is_prime(D)
+    R = coefficient_ring(D)
+    return sum(a*move_divisor(WeilDivisor(D, R; check=false)) for (D, a) in coefficient_dict(irreducible_decomposition(D)); init=WeilDivisor(X, R))
+  end
+
+  # We may assume that `D` is prime.
+  P = first(components(D))
+  # find a chart where the support is visible
+  i = findfirst(U->!isone(P(U)), affine_charts(X))
+  i === nothing && error("divisor is trivial")
+  U = affine_charts(X)[i]
+  I = P(U)
+  L, loc = localization(OO(U), complement_of_prime_ideal(I))
+  LP = loc(I)
+  # Find a principal generator for the local ideal.
+  # This works, because we assume that `X` is irreducible and reduced.
+  # Then there is at least one smooth point on `U` and therefore `LP` 
+  # is regular. Regular domains are UFD and an ideal of height 1 is 
+  # principal there. 
+  g = gens(saturated_ideal(I))
+  g = sort!(g, by=total_degree)
+  i = findfirst(f->(ideal(L, f) == LP), g)
+  x = g[i]
+  f = function_field(X; check)(x)
+  if randomization
+    kk = base_ring(X)
+    R = ambient_coordinate_ring(U)
+    y = rand(kk, 1:10) + sum(rand(kk, 1:10)*a for a in gens(R); init=zero(R))
+    f = f*inv(parent(f)(y))
+  end
+  result = irreducible_decomposition(D - weil_divisor(f))
+  # Check whether the supports are really different
+  if any(any(P == Q for Q in components(D)) for P in components(result))
+    return move_divisor(D; randomization=true, check, is_prime)
+  end
+  return result
+end
+
+function is_zero(D::AbsAlgebraicCycle)
+  all(is_zero(c) || is_one(I) for (I, c) in coefficient_dict(D)) && return true
+  return all(is_zero(c) || is_one(I) for (I, c) in coefficient_dict(irreducible_decomposition(D))) && return true
+end
+
+# Internal method which performs an automated move of the second argument
+# if things are not in sufficiently general position. 
+# The problem is: We have no way to check whether a variety is proper over 
+# its base field. But the output is only valid if that is true. 
+# Therefore, we have no choice, but to hide it from the user for now.
+function _intersect(C::CartierDivisor, D::AbsWeilDivisor)
+  X = scheme(C)
+  @assert X === scheme(D)
+  R = coefficient_ring(C)
+  @assert R === coefficient_ring(D)
+  result = AlgebraicCycle(X, R)
+  for (E, c) in coefficient_dict(C)
+    result = result + c*_intersect(E, D)
+  end
+  return result
+end
+
+function _intersect(E::EffectiveCartierDivisor, D::AbsWeilDivisor; check::Bool=true)
+  X = scheme(E)
+  @assert X === scheme(D)
+  R = coefficient_ring(D)
+  result = AlgebraicCycle(X, R)
+  cpcd = copy(coefficient_dict(D))
+  DD = irreducible_decomposition(D)
+  cpcd = copy(coefficient_dict(DD))
+  for (P, a) in coefficient_dict(DD)
+    _, inc_P = sub(P)
+    # WARNING: The heuristic implemented below might still run into an infinite loop!
+    # We have to think about how this can effectively be avoided. See the test file 
+    # for an example.
+    if is_zero(pullback(inc_P, ideal_sheaf(E)))
+      P_moved = move_divisor(WeilDivisor(X, R, 
+        IdDict{AbsIdealSheaf, elem_type(R)}(
+          [P=>one(R)]); check=false); check=false, randomization=true, is_prime=true)
+      result = result + a*_intersect(E, P_moved)
+    else
+      result = result + a*AlgebraicCycle(X, R, 
+        IdDict{AbsIdealSheaf, elem_type(R)}(
+          [P + ideal_sheaf(E)=>one(R)]); check=false)
+    end
+  end
+  return result
+end
+
+
diff --git a/src/AlgebraicGeometry/Schemes/FunctionField/FunctionFields.jl b/src/AlgebraicGeometry/Schemes/FunctionField/FunctionFields.jl
index 91c018e0b6e0..710e252f4d97 100644
--- a/src/AlgebraicGeometry/Schemes/FunctionField/FunctionFields.jl
+++ b/src/AlgebraicGeometry/Schemes/FunctionField/FunctionFields.jl
@@ -517,3 +517,6 @@ function pullback(f::AbsCoveredSchemeMorphism, a::VarietyFunctionFieldElem)
   end
 end
 
+scheme(f::VarietyFunctionFieldElem) = scheme(parent(f))
+variety(f::VarietyFunctionFieldElem) = variety(parent(f))
+
diff --git a/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl b/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
index bba517f1af49..3c59da570bf5 100644
--- a/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
+++ b/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl
@@ -757,7 +757,7 @@ is_locally_prime(I::PrimeIdealSheafFromChart) = true
 
 function is_equidimensional(I::AbsIdealSheaf; covering=default_covering(scheme(I)))
   has_attribute(I, :is_equidimensional) && return get_attribute(I, :is_equidimensional)::Bool
-  has_attribute(I, :is_prime) && return get_attribute(I, :is_prime)::Bool
+  has_attribute(I, :is_prime) && get_attribute(I, :is_prime)::Bool && return true
   local_dims = [dim(I(U)) for U in patches(covering) if !isone(I(U))]
   length(local_dims) == 0 && return true # This only happens if I == OO(X)
   d = first(local_dims)
@@ -2329,3 +2329,12 @@ function colength(I::AbsIdealSheaf; covering::Covering=default_covering(scheme(I
   end
   return result
 end
+
+function is_zero(II::AbsIdealSheaf)
+  return all(iszero(II(U)) for U in affine_charts(scheme(II)))
+end
+
+function is_zero(II::PrimeIdealSheafFromChart)
+  return is_zero(II(original_chart(II)))
+end
+
diff --git a/test/AlgebraicGeometry/Schemes/self_intersection.jl b/test/AlgebraicGeometry/Schemes/self_intersection.jl
new file mode 100644
index 000000000000..af50152a8972
--- /dev/null
+++ b/test/AlgebraicGeometry/Schemes/self_intersection.jl
@@ -0,0 +1,27 @@
+@testset "self intersection" begin
+  IP3 = projective_space(QQ, [:x, :y, :z, :w])
+  S = homogeneous_coordinate_ring(IP3)
+  (x, y, z, w) = gens(S)
+  # The following line produces an example where the current heuristic 
+  # for moving around divisors to general position still runs into an 
+  # infinite loop. We would like to solve this problem in the long run.
+  #I = ideal(S, x^3 + y*z*w)
+  # This example is fine, though.
+  I = ideal(S, x^2 + y*z)
+  IPX, inc_IPX = sub(IP3, I)
+  X = covered_scheme(IPX)
+  J = ideal(homogeneous_coordinate_ring(IPX), [x, y, z])
+  JJ = ideal_sheaf(IPX, J)
+  bl = blow_up(JJ)
+  Y = domain(bl)
+  @test is_smooth(Y)
+  E = exceptional_divisor(bl)
+  E_weil = irreducible_decomposition(weil_divisor(E))
+  @test integral(Oscar._intersect(E, E_weil)) == -2
+  mov = Oscar.move_divisor(E_weil)
+  mov = Oscar.move_divisor(mov)
+  # The following line will also produce an infinite loop. TODO: Fix this!
+  #mov = Oscar.move_divisor(Oscar.move_divisor(Oscar.move_divisor(E_weil)))
+  @test integral(Oscar._intersect(E, mov)) == -2
+end
+