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

Code before:
{
    auto t = ConstRealArray1D(tptr, len_t);
    auto c = ConstRealArray2D(cptr, n, m);
    auto xp = ConstRealArray1D(xp_ptr, s);
    auto out = RealArray2D(out_ptr, s, m);

    int64_t interval = k;
    for(int64_t ip=0; ip < s; ip++) {

        double xval = xp(ip);
        interval = _find_interval(t.data, len_t, k, xval, interval, extrapolate);
        if (interval < 0) {
            // xval was nan etc
            for (int64_t jp=0; jp < m; jp++) {
                out(ip, jp) = std::numeric_limits<double>::quiet_NaN();
            }
        }
        else {
            // Evaluate (k+1) b-splines which are non-zero on the interval.
            // on return, first k+1 elements of work are B_{m-k},..., B_{m}
            if (nu > k + 1) {
                throw std::runtime_error(
                    "Requested derivative order nu exceeds the maximum allowed (k+1). "
                    "Cannot compute derivative of order nu for spline of order k."
                );
            }
            _deBoor_D(t.data, xval, k, interval, nu, wrk);

            // Form linear combinations
            for (int64_t jp=0; jp < m; jp++) {
                out(ip, jp) = 0.0;
                for (int64_t a=0; a < k+1 ; a++){
                    out(ip, jp) += c(interval + a -k, jp) * wrk[a];
                }
            }
        }
    }
}


/*
 * Spline colocation matrix in the LAPACK banded storage
 */
void
_coloc_matrix(const double *xptr, int64_t m,       // x, shape(m,)
              const double *tptr, int64_t len_t,   // t, shape(len_t,)
              int k,
              double *abT_ptr, int64_t nbands,     // ab(nbands, len_t - k - 1) in F order!
              int offset,
              double *wrk                          // scratch, shape (2k+2)
)

Review comment (suggestion):
This is not desired on two counts:

1. If you evaluate a spline on an array of 100500 points, this runs the check 100500 times, instead of just once.

2. Evaluating a too-high derivative of a cubic spline should return zeros, not give a RuntimeError:

```
In [1]: from scipy.interpolate import make_interp_spline

In [2]: import numpy as np

In [3]: spl = make_interp_spline(np.arange(8), np.arange(8))

In [4]: spl.k
Out[4]: 3

In [5]: spl(1)
Out[5]: array(1.)

In [6]: spl(1, nu=8)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[6], line 1
----> 1 spl(1, nu=8)

File ~/repos/scipy/scipy/build-install/usr/lib/python3.12/site-packages/scipy/interpolate/_bsplines.py:564, in BSpline.__call__(self, x, nu, extrapolate)
    561 cc = cc.reshape(cc.shape[0], -1)
    563 # heavy lifting: actually perform the evaluations
--> 564 out = _dierckx.evaluate_spline(self._t, cc, self.k, x, nu, extrapolate)
    566 if is_complex:
    567     out = out.view(complex)

RuntimeError: Requested derivative order nu exceeds the maximum allowed (k+1). Cannot compute derivative of order nu for spline of order k.

```

Code after:
{
    auto t = ConstRealArray1D(tptr, len_t);
    auto c = ConstRealArray2D(cptr, n, m);
    auto xp = ConstRealArray1D(xp_ptr, s);
    auto out = RealArray2D(out_ptr, s, m);

    /*
     * If nu > k+1, the B-spline derivative is identically zero.
     * Avoid the de Boor recursion (which assumes nu <= k+1) and
     * return a zero-filled result directly.
     */
    if( nu > k + 1 ) {
        for( int64_t i = 0; i < s; i++ ) {
            for( int64_t j = 0; j < m; j++ ) {
                out(i, j) = 0.0;
            }
        }
        return ;
    }

    int64_t interval = k;
    for(int64_t ip=0; ip < s; ip++) {

        double xval = xp(ip);
        interval = _find_interval(t.data, len_t, k, xval, interval, extrapolate);
        if (interval < 0) {
            // xval was nan etc
            for (int64_t jp=0; jp < m; jp++) {
                out(ip, jp) = std::numeric_limits<double>::quiet_NaN();
            }
        }
        else {
            // Evaluate (k+1) b-splines which are non-zero on the interval.
            // on return, first k+1 elements of work are B_{m-k},..., B_{m}
            _deBoor_D(t.data, xval, k, interval, nu, wrk);

            // Form linear combinations
            for (int64_t jp=0; jp < m; jp++) {
                out(ip, jp) = 0.0;
                for (int64_t a=0; a < k+1 ; a++){
                    out(ip, jp) += c(interval + a -k, jp) * wrk[a];
                }
            }
        }
    }
}


/*
 * Spline colocation matrix in the LAPACK banded storage
 */