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/cusp_areas_from_matrix.py b/python/cusps/cusp_areas_from_matrix.py
index bf53d58e..21627dce 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)]]))
@@ -35,7 +35,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/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/sage_helper.py b/python/sage_helper.py
index 8f9f382d..c75e7768 100644
--- a/python/sage_helper.py
+++ b/python/sage_helper.py
@@ -41,6 +41,10 @@ def sage_method(function):
         from sage.all import gap
     except ImportError:
         # Modularized Sage library
+        try:
+            import sage.all__sagemath_modules
+        except ImportError:
+            pass
         from sage.algebras.group_algebra import GroupAlgebra
         from sage.arith.misc import (gcd, xgcd, is_prime, is_prime_power,
                                      is_power_of_two, prime_range, prime_powers)
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):
diff --git a/setup.py b/setup.py
index 71f18a85..f1a8fe21 100644
--- a/setup.py
+++ b/setup.py
@@ -562,12 +562,8 @@ def split_version(s : str):
                     'pypng', # Used to save OpenGL images.
                     'PyX', # Used to save PDF images of links.
                     'pickleshare', # To avoid https://github.com/ipython/ipython/issues/14416
+                    'ipython>=5.0',
                     ]
-try:
-    import sage
-except ImportError:
-    install_requires.append('cypari>=2.3')
-    install_requires.append('ipython>=5.0')
 
 # Get version number:
 exec(open('python/version.py').read())
@@ -583,6 +579,22 @@ def split_version(s : str):
        force = True,
        python_requires = '>=3.8',
        install_requires = install_requires,
+       extras_require = {
+           "cypari": [
+               "cypari>=2.3",
+           ],
+           "passagemath": [
+               "fpylll",
+               "passagemath-pari",       # provides cypari2
+               "passagemath-combinat",   # sage.algebras.free_algebra
+               "passagemath-flint",      # for sage.rings.imaginary_unit
+               "passagemath-repl",       # --> sage_eval
+               "passagemath-groups",
+               "passagemath-symbolics",  # for sage.symbolic.constants
+               "passagemath-plot",
+               "passagemath-graphs",
+            ],
+       },
        packages = ['snappy',
                    'snappy/manifolds',
                    'snappy/twister',
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 00000000..6e1d75c6
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,8 @@
+[tox]
+envlist = passagemath
+
+[testenv:passagemath]
+extras = passagemath
+
+commands =
+     python -m snappy.test