diff --git a/docs/algorithms/curve-curve-intersection.rst b/docs/algorithms/curve-curve-intersection.rst
index d3902031..48883ad4 100644
--- a/docs/algorithms/curve-curve-intersection.rst
+++ b/docs/algorithms/curve-curve-intersection.rst
@@ -11,7 +11,7 @@ Curve-Curve Intersection
            return -np.inf
        _, result = np.frexp(value)
        # Shift [1/2, 1) --> [1, 2) borrows one from exponent
-       return result - 1
+       return int(result - 1)
 
 The problem of intersecting two curves is a difficult one
 in computational geometry. The :meth:`.Curve.intersect` method (when using
diff --git a/src/python/bezier/curve.py b/src/python/bezier/curve.py
index 94c4e9b9..080024dc 100644
--- a/src/python/bezier/curve.py
+++ b/src/python/bezier/curve.py
@@ -28,7 +28,7 @@ def binary_exponent(value):
            return -np.inf
        _, result = np.frexp(value)
        # Shift [1/2, 1) --> [1, 2) borrows one from exponent
-       return result - 1
+       return int(result - 1)
 """
 
 import numpy as np
@@ -764,10 +764,10 @@ def specialize(self, start, end):
 
            >>> left, right = curve.subdivide()
            >>> also_left = curve.specialize(0.0, 0.5)
-           >>> np.all(also_left.nodes == left.nodes)
+           >>> bool(np.all(also_left.nodes == left.nodes))
            True
            >>> also_right = curve.specialize(0.5, 1.0)
-           >>> np.all(also_right.nodes == right.nodes)
+           >>> bool(np.all(also_right.nodes == right.nodes))
            True
 
         Args:
diff --git a/src/python/bezier/hazmat/algebraic_intersection.py b/src/python/bezier/hazmat/algebraic_intersection.py
index 9c7225a6..70364fbb 100644
--- a/src/python/bezier/hazmat/algebraic_intersection.py
+++ b/src/python/bezier/hazmat/algebraic_intersection.py
@@ -976,7 +976,7 @@ def lu_companion(top_row, value):
        array([[ 1. , -0.5,  0. ],
               [ 0. ,  1. , -0.5],
               [-3.5,  1. , -2. ]])
-       >>> np.linalg.norm(a_mat, ord=1)
+       >>> float(np.linalg.norm(a_mat, ord=1))
        4.5
 
     Args:
@@ -999,7 +999,7 @@ def lu_companion(top_row, value):
     lu_mat = np.zeros((degree, degree), order="F")
     if degree == 1:
         lu_mat[0, 0] = top_row[0] - value
-        return lu_mat, abs(lu_mat[0, 0])
+        return lu_mat, float(abs(lu_mat[0, 0]))
 
     # Column 0: Special case since it doesn't have ``-t`` above the diagonal.
     horner_curr = top_row[0] - value
@@ -1022,7 +1022,7 @@ def lu_companion(top_row, value):
     one_norm = max(one_norm, abs(value) + abs(curr_coeff))
     lu_mat[last_row - 1, last_row] = -value
     lu_mat[last_row, last_row] = horner_curr
-    return lu_mat, one_norm
+    return lu_mat, float(one_norm)
 
 
 def _reciprocal_condition_number(lu_mat, one_norm):
diff --git a/src/python/bezier/hazmat/clipping.py b/src/python/bezier/hazmat/clipping.py
index 89ceda5e..d8867c01 100644
--- a/src/python/bezier/hazmat/clipping.py
+++ b/src/python/bezier/hazmat/clipping.py
@@ -102,7 +102,7 @@ def compute_implicit_line(nodes):
     coeff_b = delta[0]
     # c = - ax - by = delta[1] x - delta[0] y
     coeff_c = delta[1] * nodes[0, 0] - delta[0] * nodes[1, 0]
-    return coeff_a, coeff_b, coeff_c
+    return float(coeff_a), float(coeff_b), float(coeff_c)
 
 
 def compute_fat_line(nodes):
@@ -172,7 +172,7 @@ def compute_fat_line(nodes):
             d_min = curr_dist
         elif curr_dist > d_max:
             d_max = curr_dist
-    return coeff_a, coeff_b, coeff_c, d_min, d_max
+    return coeff_a, coeff_b, coeff_c, float(d_min), float(d_max)
 
 
 def _update_parameters(s_min, s_max, start0, end0, start1, end1):
diff --git a/src/python/bezier/hazmat/curve_helpers.py b/src/python/bezier/hazmat/curve_helpers.py
index 7c900f81..c62c0d1e 100644
--- a/src/python/bezier/hazmat/curve_helpers.py
+++ b/src/python/bezier/hazmat/curve_helpers.py
@@ -668,7 +668,7 @@ def get_curvature(nodes, tangent_vec, s):
     )
     # NOTE: We convert to 1D to make sure NumPy uses vector norm.
     curvature /= np.linalg.norm(tangent_vec[:, 0], ord=2) ** 3
-    return curvature
+    return float(curvature)
 
 
 def newton_refine(nodes, point, s):
@@ -765,22 +765,22 @@ def newton_refine(nodes, point, s):
        array([[0.],
               [0.]])
        >>> s_vals = [0.625, None, None, None, None, None]
-       >>> np.log2(abs(expected - s_vals[0]))
+       >>> float(np.log2(abs(expected - s_vals[0])))
        -3.0
        >>> s_vals[1] = newton_refine(nodes, point, s_vals[0])
-       >>> np.log2(abs(expected - s_vals[1]))
+       >>> float(np.log2(abs(expected - s_vals[1])))
        -3.983...
        >>> s_vals[2] = newton_refine(nodes, point, s_vals[1])
-       >>> np.log2(abs(expected - s_vals[2]))
+       >>> float(np.log2(abs(expected - s_vals[2])))
        -4.979...
        >>> s_vals[3] = newton_refine(nodes, point, s_vals[2])
-       >>> np.log2(abs(expected - s_vals[3]))
+       >>> float(np.log2(abs(expected - s_vals[3])))
        -5.978...
        >>> s_vals[4] = newton_refine(nodes, point, s_vals[3])
-       >>> np.log2(abs(expected - s_vals[4]))
+       >>> float(np.log2(abs(expected - s_vals[4])))
        -6.978...
        >>> s_vals[5] = newton_refine(nodes, point, s_vals[4])
-       >>> np.log2(abs(expected - s_vals[5]))
+       >>> float(np.log2(abs(expected - s_vals[5])))
        -7.978...
 
     .. testcleanup:: newton-refine-curve-cusp
@@ -842,7 +842,7 @@ def newton_refine(nodes, point, s):
     delta_s = np.vdot(pt_delta[:, 0], derivative[:, 0]) / np.vdot(
         derivative[:, 0], derivative[:, 0]
     )
-    return s + delta_s
+    return float(s + delta_s)
 
 
 def locate_point(nodes, point):
diff --git a/src/python/bezier/hazmat/geometric_intersection.py b/src/python/bezier/hazmat/geometric_intersection.py
index 2e6a2090..628bde96 100644
--- a/src/python/bezier/hazmat/geometric_intersection.py
+++ b/src/python/bezier/hazmat/geometric_intersection.py
@@ -240,7 +240,7 @@ def linearization_error(nodes):
     # max_{0 <= s <= 1} s(1 - s)/2 = 1/8 = 0.125
     multiplier = 0.125 * degree * (degree - 1)
     # NOTE: worst_case is 1D due to np.max(), so this is the vector norm.
-    return multiplier * np.linalg.norm(worst_case, ord=2)
+    return float(multiplier * np.linalg.norm(worst_case, ord=2))
 
 
 def segment_intersection(start0, end0, start1, end1):
@@ -408,7 +408,7 @@ def segment_intersection(start0, end0, start1, end1):
         start_delta = start1 - start0
         s = _py_helpers.cross_product(start_delta, delta1) / cross_d0_d1
         t = _py_helpers.cross_product(start_delta, delta0) / cross_d0_d1
-        return s, t, True
+        return float(s), float(t), True
 
 
 def parallel_lines_parameters(start0, end0, start1, end1):
diff --git a/src/python/bezier/hazmat/intersection_helpers.py b/src/python/bezier/hazmat/intersection_helpers.py
index 646dcf65..e133d18b 100644
--- a/src/python/bezier/hazmat/intersection_helpers.py
+++ b/src/python/bezier/hazmat/intersection_helpers.py
@@ -187,16 +187,16 @@ def realroots(*coeffs):
        >>> expected, = realroots(28, -30, 9, -1)
        >>> s_vals = [0.625, None, None, None, None]
        >>> t = 0.625
-       >>> np.log2(abs(expected - s_vals[0]))
+       >>> float(np.log2(abs(expected - s_vals[0])))
        -4.399...
        >>> s_vals[1], t = newton_refine(s_vals[0], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[1]))
+       >>> float(np.log2(abs(expected - s_vals[1])))
        -7.901...
        >>> s_vals[2], t = newton_refine(s_vals[1], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[2]))
+       >>> float(np.log2(abs(expected - s_vals[2])))
        -16.010...
        >>> s_vals[3], t = newton_refine(s_vals[2], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[3]))
+       >>> float(np.log2(abs(expected - s_vals[3])))
        -32.110...
        >>> s_vals[4], t = newton_refine(s_vals[3], nodes1, t, nodes2)
        >>> np.allclose(s_vals[4], expected, rtol=6 * machine_eps, atol=0.0)
@@ -229,22 +229,22 @@ def realroots(*coeffs):
        >>> expected = 0.5
        >>> s_vals = [0.375, None, None, None, None, None]
        >>> t = 0.375
-       >>> np.log2(abs(expected - s_vals[0]))
+       >>> float(np.log2(abs(expected - s_vals[0])))
        -3.0
        >>> s_vals[1], t = newton_refine(s_vals[0], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[1]))
+       >>> float(np.log2(abs(expected - s_vals[1])))
        -4.0
        >>> s_vals[2], t = newton_refine(s_vals[1], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[2]))
+       >>> float(np.log2(abs(expected - s_vals[2])))
        -5.0
        >>> s_vals[3], t = newton_refine(s_vals[2], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[3]))
+       >>> float(np.log2(abs(expected - s_vals[3])))
        -6.0
        >>> s_vals[4], t = newton_refine(s_vals[3], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[4]))
+       >>> float(np.log2(abs(expected - s_vals[4])))
        -7.0
        >>> s_vals[5], t = newton_refine(s_vals[4], nodes1, t, nodes2)
-       >>> np.log2(abs(expected - s_vals[5]))
+       >>> float(np.log2(abs(expected - s_vals[5])))
        -8.0
 
     .. testcleanup:: newton-refine3
@@ -276,12 +276,12 @@ def realroots(*coeffs):
     .. doctest:: newton-refine3-continued
 
        >>> s1 = t1 = 0.5 - 0.5**27
-       >>> np.log2(0.5 - s1)
+       >>> float(np.log2(0.5 - s1))
        -27.0
        >>> s2, t2 = newton_refine(s1, nodes1, t1, nodes2)
        >>> s2 == t2
        True
-       >>> np.log2(0.5 - s2)
+       >>> float(np.log2(0.5 - s2))
        -28.0
        >>> s3, t3 = newton_refine(s2, nodes1, t2, nodes2)
        >>> s3 == t3 == s2
@@ -364,24 +364,24 @@ def modified_update(s, t):
            RHS = helpers.matrix_product(DG_t, minus_G)
            delta_params = np.linalg.solve(LHS, RHS)
            delta_s, delta_t = delta_params.flatten()
-           return s + delta_s, t + delta_t
+           return float(s + delta_s), float(t + delta_t)
 
     .. doctest:: newton-refine4
 
        >>> s0, t0 = 0.375, 0.375
-       >>> np.log2(0.5 - s0)
+       >>> float(np.log2(0.5 - s0))
        -3.0
        >>> s1, t1 = modified_update(s0, t0)
        >>> s1 == t1
        True
        >>> 1040.0 * s1
        519.0
-       >>> np.log2(0.5 - s1)
+       >>> float(np.log2(0.5 - s1))
        -10.022...
        >>> s2, t2 = modified_update(s1, t1)
        >>> s2 == t2
        True
-       >>> np.log2(0.5 - s2)
+       >>> float(np.log2(0.5 - s2))
        -31.067...
        >>> s3, t3 = modified_update(s2, t2)
        >>> s3 == t3 == 0.5
@@ -417,7 +417,7 @@ def modified_update(s, t):
     if singular:
         raise ValueError("Jacobian is singular.")
 
-    return s + delta_s, t + delta_t
+    return float(s + delta_s), float(t + delta_t)
 
 
 class NewtonSimpleRoot:  # pylint: disable=too-few-public-methods
diff --git a/src/python/bezier/hazmat/triangle_intersection.py b/src/python/bezier/hazmat/triangle_intersection.py
index ef24225f..1fb6680b 100644
--- a/src/python/bezier/hazmat/triangle_intersection.py
+++ b/src/python/bezier/hazmat/triangle_intersection.py
@@ -168,7 +168,7 @@ def newton_refine(nodes, degree, x_val, y_val, s, t):
        >>> triangle.is_valid
        True
        >>> (x_val,), (y_val,) = triangle.evaluate_cartesian(0.25, 0.5)
-       >>> x_val, y_val
+       >>> float(x_val), float(y_val)
        (1.25, 1.25)
        >>> s, t = 0.5, 0.25
        >>> new_s, new_t = newton_refine(nodes, 2, x_val, y_val, s, t)
@@ -215,7 +215,7 @@ def newton_refine(nodes, degree, x_val, y_val, s, t):
     delta_s, delta_t = newton_refine_solve(
         jac_both, x_val, triangle_x, y_val, triangle_y
     )
-    return s + delta_s, t + delta_t
+    return float(s + delta_s), float(t + delta_t)
 
 
 def update_locate_candidates(candidate, next_candidates, x_val, y_val, degree):
diff --git a/src/python/bezier/triangle.py b/src/python/bezier/triangle.py
index 5e7da5a2..8d6d9034 100644
--- a/src/python/bezier/triangle.py
+++ b/src/python/bezier/triangle.py
@@ -838,7 +838,7 @@ def _compute_valid(self):
                 self._degree, supported=(1, 2, 3)
             )
 
-        return poly_sign == 1
+        return bool(poly_sign == 1)
 
     @property
     def is_valid(self):