diff --git a/doc_src/verify.rst b/doc_src/verify.rst
index 04fc2c46..241a62d7 100644
--- a/doc_src/verify.rst
+++ b/doc_src/verify.rst
@@ -164,7 +164,7 @@ Some examples of verified computations are:
 
 * An example of finding all geodesics up to length 1::
 
-    sage: from sage.all import RIF
+    sage: from snappy.sage_helper import RIF
     sage: L = RIF(1)
     sage: M = Manifold("m003")
     sage: spec = M.length_spectrum_alt_gen(verified=True)
diff --git a/python/cusps/__init__.py b/python/cusps/__init__.py
index 76b2c246..a13f6d10 100644
--- a/python/cusps/__init__.py
+++ b/python/cusps/__init__.py
@@ -210,7 +210,7 @@ def short_slopes(manifold,
     The given :attr:`length` is cast to a SageMath ``RealIntervalField`` of the
     given precision if :attr:`verified = True`::
 
-        sage: from sage.all import pi
+        sage: from sage.symbolic.constants import pi
         sage: M.short_slopes(length = 2 * pi, verified = True, bits_prec = 100)
         [[(1, 0), (-5, 1), (-4, 1), (-3, 1), (-2, 1), (-1, 1), (0, 1), (1, 1), (2, 1), (3, 1), (4, 1), (5, 1)]]
 
diff --git a/python/cusps/cusp_areas_from_matrix.py b/python/cusps/cusp_areas_from_matrix.py
index 45b9fb21..468e4fef 100644
--- a/python/cusps/cusp_areas_from_matrix.py
+++ b/python/cusps/cusp_areas_from_matrix.py
@@ -8,7 +8,7 @@ def unbiased_cusp_areas_from_cusp_area_matrix(cusp_area_matrix):
 
     Examples::
 
-        sage: from sage.all import matrix, RIF
+        sage: from snappy.sage_helper import matrix, RIF
         sage: unbiased_cusp_areas_from_cusp_area_matrix(
         ...             matrix([[RIF(9.0,9.0005),RIF(6.0, 6.001)],
         ...                     [RIF(6.0,6.001 ),RIF(4.0, 4.001)]]))
@@ -64,7 +64,7 @@ def unbiased_cusp_areas_from_cusp_area_matrix(cusp_area_matrix):
 def greedy_cusp_areas_from_cusp_area_matrix(cusp_area_matrix, first_cusps=[]):
     """
 
-        sage: from sage.all import matrix, RIF
+        sage: from snappy.sage_helper import matrix, RIF
         sage: greedy_cusp_areas_from_cusp_area_matrix(
         ...             matrix([[RIF(9.0,9.0005),RIF(6.0, 6.001)],
         ...                     [RIF(6.0,6.001 ),RIF(10.0, 10.001)]]))
diff --git a/python/cusps/short_slopes_for_cusp.py b/python/cusps/short_slopes_for_cusp.py
index de4787a4..65d190d1 100644
--- a/python/cusps/short_slopes_for_cusp.py
+++ b/python/cusps/short_slopes_for_cusp.py
@@ -20,7 +20,7 @@ def short_slopes_from_cusp_shape_and_area(
     cusp_shape is longitude over meridian (conjugate).
     l/m
 
-    sage: from sage.all import RIF, CIF
+    sage: from snappy.sage_helper import RIF, CIF
     sage: short_slopes_from_cusp_shape_and_area(CIF(RIF(1.0),RIF(1.3333,1.3334)), RIF(12.0))
     [(1, 0), (-2, 1), (-1, 1), (0, 1)]
 
diff --git a/python/margulis/__init__.py b/python/margulis/__init__.py
index 89d832b1..208b4c5a 100644
--- a/python/margulis/__init__.py
+++ b/python/margulis/__init__.py
@@ -263,7 +263,7 @@ def margulis(
     number. Here is an example to prove that :math:`\\mu(M)>0.8`::
 
         sage: M=ManifoldHP("s479(-3,1)")
-        sage: from sage.all import RealIntervalField
+        sage: from snappy.sage_helper import RealIntervalField
         sage: M.margulis(verified=True, stopper=0.9)>RealIntervalField()('0.8')
         True
 
diff --git a/python/math_basics.py b/python/math_basics.py
index 3cad4c64..94e9942a 100644
--- a/python/math_basics.py
+++ b/python/math_basics.py
@@ -104,7 +104,7 @@ def correct_min(l):
     This is needed because python's min returns the wrong result
     for intervals, for example:
 
-        sage: from sage.all import RIF
+        sage: from snappy.sage_helper import RIF
         sage: min(RIF(4,5), RIF(3,6)).endpoints()
         (4.00000000000000, 5.00000000000000)
 
diff --git a/python/matrix.py b/python/matrix.py
index 8ef883a2..0bc16bcd 100644
--- a/python/matrix.py
+++ b/python/matrix.py
@@ -373,7 +373,7 @@ def mat_solve(m, v, epsilon=0):
 
     Setup a complex interval for example::
 
-        sage: from sage.all import RealIntervalField, ComplexIntervalField
+        sage: from snappy.sage_helper import RealIntervalField, ComplexIntervalField
         sage: RIF = RealIntervalField(80)
         sage: CIF = ComplexIntervalField(80)
         sage: fuzzy_four = CIF(RIF(3.9999,4.0001),RIF(-0.0001,0.0001))
diff --git a/python/tiling/floor.py b/python/tiling/floor.py
index 05928b23..67ef6fe0 100644
--- a/python/tiling/floor.py
+++ b/python/tiling/floor.py
@@ -8,7 +8,7 @@ def floor_as_integers(x) -> Sequence[int]:
     Computes floor of a number or interval, returning a list of integers
     if evaluating floor is ambiguous.
 
-        sage: from sage.all import RIF
+        sage: from snappy.sage_helper import RIF
         sage: floor_as_integers(RIF(1.8, 1.9))
         [1]
         sage: floor_as_integers(RIF(1.9, 2.1))
diff --git a/python/tiling/real_hash_dict.py b/python/tiling/real_hash_dict.py
index ecb7da9c..ae46edb9 100644
--- a/python/tiling/real_hash_dict.py
+++ b/python/tiling/real_hash_dict.py
@@ -58,7 +58,7 @@ class RealHashDict:
     >>> d.get(to_number(10.0000001))
     ['A', 'B']
 
-    sage: from sage.all import RIF
+    sage: from snappy.sage_helper import RIF
     sage: def equality_predicate(x, y):
     ...       d = abs(x - y)
     ...       if d < RIF(0.1):
diff --git a/python/verify/interval_newton_shapes_engine.py b/python/verify/interval_newton_shapes_engine.py
index f635f4f2..f06cca7c 100644
--- a/python/verify/interval_newton_shapes_engine.py
+++ b/python/verify/interval_newton_shapes_engine.py
@@ -350,7 +350,7 @@ def certified_newton_iteration(equations, shape_intervals,
 
         Intervals containing the true solution::
 
-            sage: from sage.all import vector
+            sage: from snappy.sage_helper import vector
             sage: good_shapes = vector([
             ...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
             ...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
@@ -368,7 +368,7 @@ def certified_newton_iteration(equations, shape_intervals,
 
         Intervals not containing a true solution::
 
-            sage: from sage.all import vector
+            sage: from snappy.sage_helper import vector
             sage: bad_shapes = vector([
             ...       C.CIF(C.RIF(0.78054, 0.78055), C.RIF(0.91447, 0.91448)),
             ...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
diff --git a/python/verify/interval_tree.py b/python/verify/interval_tree.py
index 08ae51dd..b57dbda8 100644
--- a/python/verify/interval_tree.py
+++ b/python/verify/interval_tree.py
@@ -12,7 +12,7 @@ class IntervalTree():
 
     Create an interval tree and add pairs to it::
 
-        sage: from sage.all import RIF
+        sage: from snappy.sage_helper import RIF
         sage: t = IntervalTree()
         sage: t.insert(RIF(1.01,1.02),'1')
         sage: t.insert(RIF(3.01,3.02),'3')
diff --git a/python/verify/krawczyk_shapes_engine.py b/python/verify/krawczyk_shapes_engine.py
index 1313a3db..5843a401 100644
--- a/python/verify/krawczyk_shapes_engine.py
+++ b/python/verify/krawczyk_shapes_engine.py
@@ -288,7 +288,7 @@ def krawczyk_interval(self, shape_intervals):
 
         Do several Krawczyk operations to get a better solution::
 
-            sage: from sage.all import vector
+            sage: from snappy.sage_helper import vector
             sage: M = Manifold("m019")
             sage: shapes = vector(ComplexIntervalField(53), [ 0.5+0.8j, 0.5+0.8j, 0.5+0.8j])
             sage: for i in range(15):