Repository: scipy/scipy

Issue:
BUG: numerical instability in interp1d linear interpolation
### Describe your issue.

When interp1d._call_linear(self, x_new) is invoked (which happens when scipy.interpolate.interp1d(x,y) is called, where y has dimension greater than 1), the current implementation uses a numerically unstable formulation of the interpolation calculation, which can incorrectly fail to return exact values when evaluated at existing data points due to floating point error (as shown below). This can be easily resolved by using the analytical expression for de Boor's algorithm instead (used under the hood for linear interpolation with the now-recommended make_interp_spline), which should always return exact values at data points.

This issue can cause interp1d to return a negative number when interpolating between positive numbers (as shown below).

I know that interp1d is legacy code, but this issue can cause unexpected negative numbers and therefore correctness issues (I found this bug because I was interpolating variances). Since there aren't currently plans to remove interp1d and it is still widely used, I think making this fix is reasonable (will submit a pull request).

### Reproducing Code Example

```python
import numpy as np
import scipy
from scipy import interpolate

x = np.array([0.0007499999999999, 0.002])
y = np.array([[0.0, 0.0], [0.0004164930555555557, 0.0]])

x_new = 0.002

# Issue present with interp1d with y having more than 1 dimension (delegates to _call_linear)
test_interp = interpolate.interp1d(x, y, axis = -1)
print('Incorrect with interp1d: Result should be [0,0], but interpolated value is actually:', test_interp(x_new), "\n")

# Issue is not present with interpolate.make_interp_spline:
bspl = interpolate.make_interp_spline(x, y, k=1, axis=-1)
print("Correct with make_interp_spline: Result should be [0,0], interpolated value is:", bspl(x_new), "\n")

# Issue is not present if y is 1D (delegates to numpy.interp):
y = np.array([0.0004164930555555557, 0.0])
test_interp = interpolate.interp1d(x, y, axis = -1)
print("Correct with interp1d if y is 1D: Result should be 0, interpolated value is:", test_interp(x_new), "\n")

# Issue if not present if using kind = "slinear" option (delegates to make_interp_spline)
test_interp = interpolate.interp1d(x, y, kind = "slinear", axis = -1) #create interpolation function for pdf variances
print('Correct with interp1d if "slinear" instead of "linear": Result should be [0,0], interpolated value is:', test_interp(x_new), "\n")
```

### Error message

```shell
Result should be [0,0], but interpolated value is actually: [ 0.00000000e+00 -5.42101086e-20] 

Correct with make_interp_spline: Result should be [0,0], interpolated value is: [0. 0.] 

Correct with interp1d if y is 1D: Result should be 0, interpolated value is: 0.0 

Correct with interp1d if "slinear" instead of "linear": Result should be [0,0], interpolated value is: 0.0
```

### SciPy/NumPy/Python version and system information

```shell
1.16.3 2.4.0 sys.version_info(major=3, minor=12, micro=12, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  lapack:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 3.0.1
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 15.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 15.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.1.6
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 13.4.0
  pythran:
    include directory: ../../../../../../private/var/folders/q1/lmdr88vx1k90l53fxl9p8lwc0000gn/T/build-env-w5rl_jw7/lib/python3.12/site-packages/pythran
    version: 0.18.0
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /private/var/folders/q1/lmdr88vx1k90l53fxl9p8lwc0000gn/T/build-env-w5rl_jw7/bin/python
  version: '3.12'
```

Discussion hints:

Solution patch:
diff --git a/scipy/interpolate/_interpolate.py b/scipy/interpolate/_interpolate.py
index a138e705aa7f..37e54578f7c6 100644
--- a/scipy/interpolate/_interpolate.py
+++ b/scipy/interpolate/_interpolate.py
@@ -487,7 +487,7 @@ def _call_linear(self, x_new):
         #    of x_new[n] = x[0]
         x_new_indices = x_new_indices.clip(1, len(self.x)-1).astype(int)
 
-        # 4. Calculate the slope of regions that each x_new value falls in.
+        # 4. Calculate y_new using de Boor's algorithm for linear interpolation.
         lo = x_new_indices - 1
         hi = x_new_indices
 
@@ -496,13 +496,13 @@ def _call_linear(self, x_new):
         y_lo = self._y[lo]
         y_hi = self._y[hi]
 
-        # Note that the following two expressions rely on the specifics of the
+        # Note that the following expression relies on the specifics of the
         # broadcasting semantics.
-        slope = (y_hi - y_lo) / (x_hi - x_lo)[:, None]
-
-        # 5. Calculate the actual value for each entry in x_new.
-        y_new = slope*(x_new - x_lo)[:, None] + y_lo
-
+        y_new = (
+            ((x_new - x_lo)/(x_hi - x_lo))[:, None] * y_hi
+            + ((x_hi - x_new)/(x_hi - x_lo))[:, None] * y_lo
+            )
+        
         return y_new
 
     def _call_nearest(self, x_new):


Test patch:
diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py
index 8078ac01109a..f40375a39ac5 100644
--- a/scipy/interpolate/tests/test_interpolate.py
+++ b/scipy/interpolate/tests/test_interpolate.py
@@ -235,6 +235,16 @@ def test_linear_dtypes(self):
         yp = interp1d(x, y)(x)
         xp_assert_close(yp, y, atol=1e-15)
 
+    def test_linear_numerical_stability(self):
+        # regression test for gh-24281: Using de Boor's algorithm, there 
+        # should be no floating point error for query points contained 
+        # exactly in the x input array
+        x = np.array([0.0007499999999999, 0.002])
+        y = np.array([[0.0, 0.0], [0.0004164930555555557, 0.0]])
+        yp = interp1d(x, y)(x)
+        xp_assert_close(yp, y, atol=1e-20)
+        assert np.all(yp >= 0) # less stable computations give yp[1,1] = -5e-20
+
     def test_slinear_dtypes(self):
         # regression test for gh-7273: 1D slinear interpolation fails with
         # float32 inputs