Using Fortran speedup for surface locate_point().

src/bezier/_curve_speedup.pyx

@@ -171,9 +171,9 @@ def locate_point(double[::1, :] nodes, double[::1, :] point):
         &s_approx,
     )
 
-    if s_approx == -1.0:
+    if s_approx == -1.0:  # LOCATE_MISS
         return None
-    elif s_approx == -2.0:
+    elif s_approx == -2.0:  # LOCATE_INVALID
         raise ValueError(
             'Parameters not close enough to one another')
     else:

src/bezier/_surface_intersection.pxd

@@ -18,3 +18,6 @@ cdef extern from "bezier/surface_intersection.h":
         int *num_nodes, double *nodes, int *degree,
         double *x_val, double *y_val, double *s, double *t,
         double *updated_s, double *updated_t)
+    void locate_point_surface(
+        int *num_nodes, double *nodes, int *degree,
+        double *x_val, double *y_val, double *s_val, double *t_val)

src/bezier/_surface_intersection.py

@@ -305,9 +305,14 @@ def mean_centroid(candidates):
     return sum_x / denom, sum_y / denom
 
 
-def locate_point(nodes, degree, x_val, y_val):
+def _locate_point(nodes, degree, x_val, y_val):
     r"""Locate a point on a surface.
 
+    .. note::
+
+       There is also a Fortran implementation of this function, which
+       will be used if it can be built.
+
     Does so by recursively subdividing the surface and rejecting
     sub-surfaces with bounding boxes that don't contain the point.
     After the sub-surfaces are sufficiently small, uses Newton's
@@ -360,6 +365,8 @@ def locate_point(nodes, degree, x_val, y_val):
 # pylint: disable=invalid-name
 if _surface_intersection_speedup is None:  # pragma: NO COVER
     newton_refine = _newton_refine
+    locate_point = _locate_point
 else:
     newton_refine = _surface_intersection_speedup.newton_refine
+    locate_point = _surface_intersection_speedup.locate_point
 # pylint: enable=invalid-name