Track incremental perf improvements / refactorings in `RegularGridInterpolator(..., method="linear")`

https://github.com/scipy/scipy/pull/17291 implemented a fast path for the 2D linear mode of the RegularGridInterpolator. That is an initial implementation, and this issue is to keep track of possible incremental improvements.

First, several possible experiments to see if the current implementation (see interpolate/_rgi.py and interpolate/_rgi_cython.pyx) can be further sped up or streamlined:

• experiment with the tuple grid argument in cython: try to either pad the individual elements (so that the grid becomes an N-D array) and pass around the real sizes, or bite the bullet and switch over to double** (maybe not this though). This would allow to remove last bits of python interaction from cython code, the latter would then be possible to be declare nogil. Will need to see if this helps with perf or not.
• see if we can require memoryviews to be C contiguous. This would require I think to loop over the trailing dimensions of values explicitly (not a big deal if somewhat messy) and maybe get the transposes right to make sure that cython loops respect memory locality. Also this would require some np.ascontiguousarrays, which might involve copying, so there are memory considerations. Again, the impact will need to be measured.
• see if we can simplify the cython logic: ATM both size-1 axes and nan inputs make find_indices return -1. Maybe we could simplify that. One complication is that the current fast path cannot fully replace python-level evaluate_linear because of "duck-typed arrays" (looking at the original PR from 2014 which introduced RGI, these are meant to take care of out-of-memory arrays like NetCDF files etc).

Additional points brought up in gh-17291 review:

• Consider adding a 3D fast path.
• Consider adding a generic N-D cythonized evaluation. This will hopefully obviate the need for fast paths, but so far it seems to be difficult to loop over runtime-specified number of dimensions in Cython. In a nutshell, this needs a performant no-python itertools.product replacement.

An alternative if quite bigger project would be to implement an NdBSpline class to mirror NdPPoly. This has a potential to significantly improve both linear and spline RGI modes: the former could then be replaced by an "slinear" mode , and the latter would no longer need recursive make_interp_spline calls.

Update 02/2024: https://github.com/scipy/scipy/pull/19633 has been merged in the 1.13 release cycle, so two last items (an N-D evaluation and NdBSplin) are done, and are available as mode="slinear"