Adding speedup for compute_length.

In the process, had to add a subset of QUADPACK and `d1mach.f`
in order to do adaptive quadrature to compute the length.

F77 source was retrieved directly from:

  http://www.netlib.org/quadpack/dqpsrt.f
  http://www.netlib.org/quadpack/dqk21.f
  http://www.netlib.org/quadpack/dqelg.f
  http://www.netlib.org/quadpack/dqagse.f
  http://www.netlib.org/blas/d1mach.f

It's worth noting that all of QUADPACK (and I believe `d1mach`) are
"Public Domain".

Used `dqagse` because that is what is called by `scipy.integrate.quad`
when there are no `points` (where the function behaves badly), there
is no special weight function and the bounds are finite.

I also used the default (from SciPy) `limit` of 50 subintervals in
the adaptive quadrature and both an absolute (`epsabs`) and relative
(`epsrel`) tolerance of `2^{-26}` (i.e. `sqrt(machine precision`).

It's worth noting that `setup_helpers.py` is a little unwieldy now
because there is so much copy-pasting in the dependencies.

For now, the error code `ier` is ignored in the Fortran implementation.
As comparision, it'd be at the very least printed as a warning from
the SciPy implementation.

setup_helpers.py

@@ -48,17 +48,51 @@
 useful, for example, to track changes across different systems or simple
 to make sure the build is occurring as expected.
 """
+QUADPACK_DIR = 'quadpack'
+QUADPACK_SOURCE_FILENAME = os.path.join('src', 'bezier', QUADPACK_DIR, '{}.f')
+# NOTE: QUADPACK module dependencies: order is important.
+QUADPACK_MODULES = (
+    'd1mach',
+    'dqelg',
+    'dqpsrt',
+    'dqk21',
+    'dqagse',
+)
 # NOTE: This represents the Fortran module dependency graph. Order is
 #       important both of the keys and of the dependencies that are in
 #       each value.
 FORTRAN_MODULES = collections.OrderedDict()
 FORTRAN_MODULES['types'] = ('types',)
 FORTRAN_MODULES['helpers'] = ('types', 'helpers')
-FORTRAN_MODULES['curve'] = ('types', 'helpers', 'curve')
-FORTRAN_MODULES['surface'] = ('types', 'helpers', 'curve', 'surface')
+FORTRAN_MODULES['curve'] = (
+    'types',
+    'helpers',
+    os.path.join(QUADPACK_DIR, 'd1mach'),
+    os.path.join(QUADPACK_DIR, 'dqelg'),
+    os.path.join(QUADPACK_DIR, 'dqpsrt'),
+    os.path.join(QUADPACK_DIR, 'dqk21'),
+    os.path.join(QUADPACK_DIR, 'dqagse'),
+    'curve',
+)
+FORTRAN_MODULES['surface'] = (
+    'types',
+    'helpers',
+    os.path.join(QUADPACK_DIR, 'd1mach'),
+    os.path.join(QUADPACK_DIR, 'dqelg'),
+    os.path.join(QUADPACK_DIR, 'dqpsrt'),
+    os.path.join(QUADPACK_DIR, 'dqk21'),
+    os.path.join(QUADPACK_DIR, 'dqagse'),
+    'curve',
+    'surface',
+)
 FORTRAN_MODULES['curve_intersection'] = (
     'types',
     'helpers',
+    os.path.join(QUADPACK_DIR, 'd1mach'),
+    os.path.join(QUADPACK_DIR, 'dqelg'),
+    os.path.join(QUADPACK_DIR, 'dqpsrt'),
+    os.path.join(QUADPACK_DIR, 'dqk21'),
+    os.path.join(QUADPACK_DIR, 'dqagse'),
     'curve',
     'curve_intersection',
 )
@@ -312,10 +346,15 @@ def save_journal(self):
             print(msg, file=sys.stderr)
 
     def compile_fortran_obj_files(self):
-        source_files = [
+        source_files_quadpack = [
+            QUADPACK_SOURCE_FILENAME.format(mod_name)
+            for mod_name in QUADPACK_MODULES
+        ]
+        source_files_bezier = [
             FORTRAN_SOURCE_FILENAME.format(mod_name)
             for mod_name in FORTRAN_MODULES
         ]
+        source_files = source_files_quadpack + source_files_bezier
         obj_files = self.F90_COMPILER.compile(
             source_files,
             output_dir=None,

src/bezier/_curve.pxd

@@ -61,3 +61,5 @@ cdef extern from "bezier/curve.h":
     void full_reduce(
         int *num_nodes, int *dimension, double *nodes,
         int *num_reduced_nodes, double *reduced, bool_t *not_implemented)
+    void compute_length(
+        int *num_nodes, int *dimension, double *nodes, double *length)

src/bezier/_curve_helpers.py

@@ -303,7 +303,7 @@ def vec_size(nodes, s_val):
     return np.linalg.norm(result_vec[0, :], ord=2)
 
 
-def compute_length(nodes, degree):
+def _compute_length(nodes, degree):
     r"""Approximately compute the length of a curve.
 
     .. _QUADPACK: https://en.wikipedia.org/wiki/QUADPACK
@@ -939,6 +939,7 @@ def _full_reduce(nodes):
     subdivide_nodes = _subdivide_nodes
     evaluate_multi = _evaluate_multi
     evaluate_multi_barycentric = _evaluate_multi_barycentric
+    compute_length = _compute_length
     elevate_nodes = _elevate_nodes
     specialize_curve = _specialize_curve
     evaluate_hodograph = _evaluate_hodograph
@@ -953,6 +954,7 @@ def _full_reduce(nodes):
     subdivide_nodes = _curve_speedup.subdivide_nodes
     evaluate_multi = _curve_speedup.evaluate_multi
     evaluate_multi_barycentric = _curve_speedup.evaluate_multi_barycentric
+    compute_length = _curve_speedup.compute_length
     elevate_nodes = _curve_speedup.elevate_nodes
     specialize_curve = _curve_speedup.specialize_curve
     evaluate_hodograph = _curve_speedup.evaluate_hodograph