diff --git a/src/NumFieldOrd/NfOrd/Clgp/Saturate.jl b/src/NumFieldOrd/NfOrd/Clgp/Saturate.jl
index e4b9c33ead..2ede5227c4 100644
--- a/src/NumFieldOrd/NfOrd/Clgp/Saturate.jl
+++ b/src/NumFieldOrd/NfOrd/Clgp/Saturate.jl
@@ -12,10 +12,11 @@ function mod_p(R::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}, Q::
   pp, e = Hecke.ppio(oF, p)
   #We now want the discrete logarithm of the images of the elements of R in F
   #We don't need the full group, just the quotient by p
-  #We compute a generator and cache the powers in dl
-  #Then we will compute the discrete logarithms by checking the values in dl
+  #We compute a generator and
+  # - cache the powers in dl if pp is small
+  # -  we will compute the discrete logarithms by checking the values in dl
+  #    or via BSGS
   dl = Dict{Hecke.Nemo.fpFieldElem, Int}()
-  dl[F(1)] = 0
   exp_to_test = divexact(pp, p)
   x = rand(F)
   while true
@@ -28,12 +29,21 @@ function mod_p(R::Vector{FacElem{AbsSimpleNumFieldElem, AbsSimpleNumField}}, Q::
     end
     x = rand(F)
   end
-  y = x
-  for i = 1:pp-1
-    dl[y] = i
-    y *= x
+  if nbits(pp) < 19 # seems to be a reasonable cutoff
+    y = one(F)
+    for i = 0:pp-1
+      dl[y] = i
+      y *= x
+    end
+  end
+  ma = Vector{Int}(undef, length(R))
+  for i in 1:length(R)
+    imgd = image(mF1, R[i], D[i], cached, pp)^e
+    ma[i] = get!(dl, imgd) do
+      Hecke.baby_step_giant_step(x, pp, imgd, dl) % p
+    end
   end
-  return matrix(T, 1, length(R), Int[dl[image(mF1, R[i], D[i], cached, pp)^e] % p for i in 1:length(R)])
+  return matrix(T, 1, length(R), ma)
 end
 
 #= idea
diff --git a/src/NumFieldOrd/NfOrd/ResidueRingMultGrp.jl b/src/NumFieldOrd/NfOrd/ResidueRingMultGrp.jl
index a4e5435ac2..b097317cf8 100644
--- a/src/NumFieldOrd/NfOrd/ResidueRingMultGrp.jl
+++ b/src/NumFieldOrd/NfOrd/ResidueRingMultGrp.jl
@@ -714,6 +714,8 @@ of the first step. This allows to speed up subsequent calls with
 the same $g$ and $n$.
 """
 function baby_step_giant_step(g, n, h, cache::Dict)
+  # if cache is given, it must be empty or contain
+  # *all* big steps
   @assert typeof(g) == typeof(h)
   m = ZZRingElem(isqrt(n) + 1)
   if isempty(cache)
diff --git a/test/NfOrd/Clgp.jl b/test/NfOrd/Clgp.jl
index c2fddd16d7..de7e57b614 100644
--- a/test/NfOrd/Clgp.jl
+++ b/test/NfOrd/Clgp.jl
@@ -285,4 +285,11 @@ end
       @test order(c) == 892
     end
   end
+
+  # saturation at large primes
+  
+  # the cyclotomic units are saturated at any prime l, since the class number h_149^+ is one.
+  K, = cyclotomic_real_subfield(149, "a")
+  u = Hecke._cyclotomic_units_totally_real_prime_conductor(K, 149)
+  @test nrows(Hecke.RelSaturate.compute_candidates_for_saturate(FacElem.(u[2:end]), next_prime(2^25))) == 0
 end