Repository: scipy/scipy
File: scipy/interpolate/src/__fitpack.cc (C++)

Code before:
        memcpy(hh, h, j*sizeof(double));
        h[0] = 0.0;
        for (n = 1; n <= j; n++) {
            ind = ell + n;
            xb = t[ind];
            xa = t[ind - j];
            if (xb == xa) {
                h[n] = 0.0;
                continue;
            }
            w = hh[n - 1]/(xb - xa);
            h[n - 1] += w*(xb - x);
            h[n] = w*(x - xa);
        }
    }

    /*
     * Now do m "derivative" recursions
     * to convert the values of beta into the mth derivative
     */
    if (m > k + 1) {
        throw std::runtime_error(
            "Requested derivative order m exceeds the maximum allowed (k+1). "
            "Cannot compute derivative of order m for spline of order k."
        );
    }
    for (j = k - m + 1; j <= k; j++) {
        memcpy(hh, h, j*sizeof(double));
        h[0] = 0.0;
        for (n = 1; n <= j; n++) {
            ind = ell + n;
            xb = t[ind];
            xa = t[ind - j];
            if (xb == xa) {
                h[m] = 0.0;
                continue;
            }
            w = j*hh[n - 1]/(xb - xa);
            h[n - 1] -= w;
            h[n] = w;
        }
    }
}


/*
 *  Find an interval such that t[interval] <= xval < t[interval+1].
 */
int64_t
_find_interval(const double* tptr, int64_t len_t,
               int k,

Review comment (suggestion):
While raising here is correct, this effectively puts the check inside a hot loop over the data points: typically we call evaluations for arrays of data points.
Thus, could you please remove the `if` check from `_deBoor_D`, and add it to call sites?
Grepping shows three relevant places:

```
$ git grep -n deBoor_D
scipy/interpolate/src/__fitpack.cc:16:_deBoor_D(const double *t, double x, int k, int ell, int m, double *result) {
scipy/interpolate/src/__fitpack.cc:166:        _deBoor_D(t.data, xval, k, ind, 0, wrk);
scipy/interpolate/src/__fitpack.cc:253:        _deBoor_D(t.data, xval, k, ind, 0, wrk);
scipy/interpolate/src/__fitpack.cc:858:        _deBoor_D(t.data, xval, k, ind, 0, wrk.data());
scipy/interpolate/src/__fitpack.cc:1208:            _deBoor_D(t.data, xval, k, interval, nu, wrk);    # evaluate_spline
scipy/interpolate/src/__fitpack.cc:1246:        _deBoor_D(t.data, xval, k, left, 0, wrk);
scipy/interpolate/src/__fitpack.cc:1291:        _deBoor_D(t.data, xval, k, left, 0, wrk);
scipy/interpolate/src/__fitpack.cc:1378:            _deBoor_D(td, xd, kd, i_d, nu(d), wrk.data());   # evaluate_ndbspline
scipy/interpolate/src/__fitpack.cc:1498:            _deBoor_D(td, xd, kd, i_d, 0, wrk.data());
scipy/interpolate/src/__fitpack.h:116:_deBoor_D(const double *t, double x, int k, int ell, int m, double *result);
scipy/interpolate/src/_dierckxmodule.cc:1049:    fitpack::_deBoor_D(        # eval_all_bspl
```

The last one, `eval_all_bspline` is special in that it's directly user facing, so it should just return the correct answer straight from C instead of raising an error.

Code after:
        memcpy(hh, h, j*sizeof(double));
        h[0] = 0.0;
        for (n = 1; n <= j; n++) {
            ind = ell + n;
            xb = t[ind];
            xa = t[ind - j];
            if (xb == xa) {
                h[n] = 0.0;
                continue;
            }
            w = hh[n - 1]/(xb - xa);
            h[n - 1] += w*(xb - x);
            h[n] = w*(x - xa);
        }
    }

    /*
     * Now do m "derivative" recursions
     * to convert the values of beta into the mth derivative
     */
    for (j = k - m + 1; j <= k; j++) {
        memcpy(hh, h, j*sizeof(double));
        h[0] = 0.0;
        for (n = 1; n <= j; n++) {
            ind = ell + n;
            xb = t[ind];
            xa = t[ind - j];
            if (xb == xa) {
                h[m] = 0.0;
                continue;
            }
            w = j*hh[n - 1]/(xb - xa);
            h[n - 1] -= w;
            h[n] = w;
        }
    }
}


/*
 *  Find an interval such that t[interval] <= xval < t[interval+1].
 */
int64_t
_find_interval(const double* tptr, int64_t len_t,
               int k,
               double xval,
               int64_t prev_l,
               int extrapolate)
{
    ConstRealArray1D t = ConstRealArray1D(tptr, len_t);