ENH: improve interp() speed in some cases.

The interp function was computing slopes for all intervals, even when there
were only a few points to be interpolated. Now it only does so when the
number of interpolation points exceeds the number of sample points.

2 files changed, 49 insertions, 19 deletions

diff --git a/numpy/lib/src/_compiled_base.c b/numpy/lib/src/_compiled_base.c
index a66e5da4c..2e290f322 100644
--- a/numpy/lib/src/_compiled_base.c
+++ b/numpy/lib/src/_compiled_base.c
@@ -537,7 +537,7 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
     PyObject *fp, *xp, *x;
     PyObject *left = NULL, *right = NULL;
     PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
-    npy_intp i, lenx, lenxp, indx;
+    npy_intp i, lenx, lenxp;
     double lval, rval;
     double *dy, *dx, *dz, *dres, *slopes;
 
@@ -604,29 +604,54 @@ arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
         }
     }
 
-    slopes = (double *) PyDataMem_NEW((lenxp - 1)*sizeof(double));
-    NPY_BEGIN_ALLOW_THREADS;
-    for (i = 0; i < lenxp - 1; i++) {
-        slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]);
-    }
-    for (i = 0; i < lenx; i++) {
-        indx = binary_search(dz[i], dx, lenxp);
-        if (indx == -1) {
-            dres[i] = lval;
-        }
-        else if (indx == lenxp - 1) {
-            dres[i] = dy[indx];
+    /* only pre-calculate slopes if there are relatively few of them. */
+    if (lenxp <= lenx) {
+        slopes = (double *) PyDataMem_NEW((lenxp - 1)*sizeof(double));
+        NPY_BEGIN_ALLOW_THREADS;
+        for (i = 0; i < lenxp - 1; i++) {
+            slopes[i] = (dy[i + 1] - dy[i])/(dx[i + 1] - dx[i]);
         }
-        else if (indx == lenxp) {
-            dres[i] = rval;
+        for (i = 0; i < lenx; i++) {
+            npy_intp j = binary_search(dz[i], dx, lenxp);
+
+            if (j == -1) {
+                dres[i] = lval;
+            }
+            else if (j == lenxp - 1) {
+                dres[i] = dy[j];
+            }
+            else if (j == lenxp) {
+                dres[i] = rval;
+            }
+            else {
+                dres[i] = slopes[j]*(dz[i] - dx[j]) + dy[j];
+            }
         }
-        else {
-            dres[i] = slopes[indx]*(dz[i] - dx[indx]) + dy[indx];
+        NPY_END_ALLOW_THREADS;
+        PyDataMem_FREE(slopes);
+    }
+    else {
+        NPY_BEGIN_ALLOW_THREADS;
+        for (i = 0; i < lenx; i++) {
+            npy_intp j = binary_search(dz[i], dx, lenxp);
+
+            if (j == -1) {
+                dres[i] = lval;
+            }
+            else if (j == lenxp - 1) {
+                dres[i] = dy[j];
+            }
+            else if (j == lenxp) {
+                dres[i] = rval;
+            }
+            else {
+                double slope = (dy[j + 1] - dy[j])/(dx[j + 1] - dx[j]);
+                dres[i] = slope*(dz[i] - dx[j]) + dy[j];
+            }
         }
+        NPY_END_ALLOW_THREADS;
     }
 
-    NPY_END_ALLOW_THREADS;
-    PyDataMem_FREE(slopes);
     Py_DECREF(afp);
     Py_DECREF(axp);
     Py_DECREF(ax);
diff --git a/numpy/lib/tests/test_function_base.py b/numpy/lib/tests/test_function_base.py
index 145c83b77..67a5204a1 100644
--- a/numpy/lib/tests/test_function_base.py
+++ b/numpy/lib/tests/test_function_base.py
@@ -1179,6 +1179,11 @@ class TestInterp(TestCase):
         x0 = np.array(.3, dtype=object)
         assert_almost_equal(np.interp(x0, x, y), .3)
 
+    def test_if_len_x_is_small(self):
+        xp = np.arange(0, 1000, 0.0001)
+        fp = np.sin(xp)
+        assert_almost_equal(np.interp(np.pi, xp, fp), 0.0)
+
 
 def compare_results(res, desired):
     for i in range(len(desired)):