diff --git a/cython/core/basic.pyx b/cython/core/basic.pyx
index b910d27d..db3e3680 100644
--- a/cython/core/basic.pyx
+++ b/cython/core/basic.pyx
@@ -52,6 +52,7 @@ from . import twister
 from . import snap
 from . import verify
 from .verify import complex_volume as verify_complex_volume
+from .cusps import compute_cusp_shapes as cusps_compute_cusp_shapes
 from . import decorated_isosig
 from .ptolemy import manifoldMethods as ptolemyManifoldMethods
 from .export_stl import stl
diff --git a/cython/core/manifold.pyx b/cython/core/manifold.pyx
index 26a281ca..1425db06 100644
--- a/cython/core/manifold.pyx
+++ b/cython/core/manifold.pyx
@@ -1176,7 +1176,7 @@ cdef class Manifold(Triangulation):
                                 "of cusp_info only "
                                 "implemented for cusp shapes. Pass 'shape' "
                                 "as first argument to cusp_info().")
-            return verify.compute_cusp_shapes(self, verified = verified,
+            return cusps_compute_cusp_shapes(self, verified = verified,
                                               bits_prec = bits_prec)
 
         if data_spec is None:
diff --git a/python/cusps/__init__.py b/python/cusps/__init__.py
index 1a701da0..20046539 100644
--- a/python/cusps/__init__.py
+++ b/python/cusps/__init__.py
@@ -1,4 +1,38 @@
 """
-Picking cusp neighborhoods and exceptional slopes.
-
+Computing data about cusps such as cusp matrix, shape, translations
+and exceptional slopes.
 """
+
+from ..geometric_structure.cusp_neighborhood.complex_cusp_cross_section import ComplexCuspCrossSection
+from ..verify.shapes import compute_hyperbolic_shapes
+from ..exceptions import NonorientableManifoldError
+
+def compute_cusp_shapes(manifold, verified, bits_prec=None):
+    """
+    Compute verified cusp shapes (following the SnapPea kernel convention,
+    it returns the conjugate of the quotient of the translations
+    corresponding to the longitude and meridian for each cusp.
+
+    >>> M = Manifold('s843')
+    >>> M.cusp_info('shape', bits_prec = 100) # doctest: +NUMERIC21
+    [0.46738227586341496791816972792 + 1.1903600506742881207098973751*I, 0.084187324414612694374797271558 + 1.0506945576790020048456757228*I]
+
+    sage: M = Manifold('s843')
+    sage: M.cusp_info('shape', verified = True) # doctest: +NUMERIC12
+    [0.46738227587? + 1.19036005068?*I, 0.0841873244146? + 1.0506945576790?*I]
+
+    sage: M.cusp_info('shape', verified = True, bits_prec = 100) # doctest: +NUMERIC21
+    [0.4673822758634149679181698? + 1.1903600506742881207098974?*I, 0.084187324414612694374797272? + 1.050694557679002004845675723?*I]
+    """
+
+    if not manifold.is_orientable():
+        raise NonorientableManifoldError(manifold)
+
+    shapes = compute_hyperbolic_shapes(
+        manifold, verified=verified, bits_prec=bits_prec)
+
+    # Compute cusp cross section
+    c = ComplexCuspCrossSection.fromManifoldAndShapes(manifold, shapes)
+
+    # Compute shape
+    return c.cusp_shapes()
diff --git a/python/exceptions.py b/python/exceptions.py
index d46eb5bf..0d2db4f9 100644
--- a/python/exceptions.py
+++ b/python/exceptions.py
@@ -1,13 +1,24 @@
 # Exceptions from the SnapPea kernel
 class SnapPeaFatalError(Exception):
     """
-    This exception is raised by SnapPy when the SnapPea kernel
-    encounters a fatal error.
+    Exception raised by SnapPy when the SnapPea kernel encounters a fatal
+    error.
     """
 
-
 class InsufficientPrecisionError(Exception):
     """
-    This exception is raised when a computation fails and is likely
-    to succeed if higher precision is used.
+    Exception raised when a computation fails and is likely to succeed if
+    higher precision is used.
+    """
+
+class NonorientableManifoldError(ValueError):
     """
+    Exception raised when a non-orientable manifold is given to a method
+    only supporting orientable manifolds.
+    """
+    def __init__(self, manifold):
+        self.manifold = manifold
+
+    def __str__(self):
+        return ('Computation only supports orientable manifolds but %s is '
+                'non-orientable') % self.manifold
diff --git a/python/verify/__init__.py b/python/verify/__init__.py
index 94445240..1ca0a83c 100644
--- a/python/verify/__init__.py
+++ b/python/verify/__init__.py
@@ -4,7 +4,6 @@
 from .hyperbolicity import *
 from .canonical import *
 from .cusp_translations import *
-from .cusp_shapes import *
 from .volume import *
 from .complex_volume import *
 from .interval_tree import *
diff --git a/python/verify/cusp_shapes.py b/python/verify/cusp_shapes.py
deleted file mode 100644
index 162a2b39..00000000
--- a/python/verify/cusp_shapes.py
+++ /dev/null
@@ -1,48 +0,0 @@
-from ..geometric_structure.cusp_neighborhood.complex_cusp_cross_section import ComplexCuspCrossSection
-from .shapes import compute_hyperbolic_shapes
-
-__all__ = ['NonorientableManifoldError', 'compute_cusp_shapes']
-
-
-class NonorientableManifoldError(RuntimeError):
-    """
-    Exception raised when trying to compute cusp shapes for a non-orientable
-    manifold.
-    """
-    def __init__(self, manifold):
-        self.manifold = manifold
-
-    def __str__(self):
-        return (('Cannot compute cusp shapes for non-orientable '
-                 'manifold %s') % self.manifold)
-
-
-def compute_cusp_shapes(manifold, verified, bits_prec=None):
-    """
-    Compute verified cusp shapes (following the SnapPea kernel convention,
-    it returns the conjugate of the quotient of the translations
-    corresponding to the longitude and meridian for each cusp.
-
-    >>> M = Manifold('s843')
-    >>> M.cusp_info('shape', bits_prec = 100) # doctest: +NUMERIC21
-    [0.46738227586341496791816972792 + 1.1903600506742881207098973751*I, 0.084187324414612694374797271558 + 1.0506945576790020048456757228*I]
-
-    sage: M = Manifold('s843')
-    sage: M.cusp_info('shape', verified = True) # doctest: +NUMERIC12
-    [0.46738227587? + 1.19036005068?*I, 0.0841873244146? + 1.0506945576790?*I]
-
-    sage: M.cusp_info('shape', verified = True, bits_prec = 100) # doctest: +NUMERIC21
-    [0.4673822758634149679181698? + 1.1903600506742881207098974?*I, 0.084187324414612694374797272? + 1.050694557679002004845675723?*I]
-    """
-
-    if not manifold.is_orientable():
-        raise NonorientableManifoldError(manifold)
-
-    shapes = compute_hyperbolic_shapes(
-        manifold, verified=verified, bits_prec=bits_prec)
-
-    # Compute cusp cross section
-    c = ComplexCuspCrossSection.fromManifoldAndShapes(manifold, shapes)
-
-    # Compute shape
-    return c.cusp_shapes()