如何用Numpy快速获取单调数组的近线性间距元素索引？

我有一个约10^7个元素的单调递增Numpy数组：

example_array = np.array([10, 10.5, 13, 15, 20, 35, 37.1, 40, 50])

我希望找到一组索引，使得原数组按该索引取值后，元素间距最接近指定的线性间隔（无需显式插值）。例如，当指定间隔为10.5时，期望输出索引为：

output_index = np.array([1, 4, 5, 7, 8])
print(example_array[output_index]) # [10.5, 20. , 35. , 40. , 50. ]

如何用Numpy高效实现该需求？我曾尝试取模运算：

In [1]: import numpy as np

In [2]: example_array = np.array([10, 10.5, 13, 15, 20, 35, 37.1, 40, 50])

In [3]: spacing = 10.5

In [4]: example_array/spacing

Out[4]: array([0.95238095, 1.        , 1.23809524, 1.42857143, 1.9047619 ,
       3.33333333, 3.53333333, 3.80952381, 4.76190476])

In [5]: example_array%spacing

Out[5]: array([10. ,  0. ,  2.5,  4.5,  9.5,  3.5,  5.6,  8.5,  8. ])

但无法区分目标索引，实际我需要找到每个整数倍间隔对应的最接近元素的索引（每个整数仅对应一个）。

::编辑说明::
已补充遗漏的索引，该问题与《在两个Numpy数组中查找最接近值》有相似性。

以下是测试代码（感谢AJ Biffl提供的实现）：

In [1]: import numpy as np

In [2]: def nearest_integers2D(arr, spacing):
   ...:     nmax = int(arr[-1]/spacing)
   ...:     range_arr = np.arange(1, nmax+1)
   ...:     difference = np.abs(arr[:, np.newaxis]/spacing - range_arr)
   ...:     return np.argmin(difference, axis=0)
   ...: 

In [3]: def nearest_integers_sort_fast(arr, spacing):
   ...:     nmin = int(arr[0]//spacing) + 1
   ...:     nmax = int(arr[-1]//spacing) + 1
   ...:     idx = np.zeros(nmax - nmin + 1, dtype = int)
   ...:     n = nmin
   ...:     dprev = spacing
   ...:     xprev = -nmin*spacing
   ...:     for i, x in enumerate(arr):
   ...:         d = abs(n*spacing - x)
   ...:         if d > dprev:
   ...:             idx[n-nmin] = i-1
   ...:             if n < len(arr):
   ...:                 n += 1
   ...:                 d = abs(n*spacing - x)
   ...:                 dprev += abs(n*spacing - xprev)
   ...:         else:
   ...:             dprev = d
   ...:         xprev = x
   ...:     idx[-1] = len(arr)-1
   ...:     return idx
   ...: 

In [4]: example_array = np.array([10, 10.5, 13, 15, 20, 35, 37.1, 40, 50])

In [5]: rand_array = np.cumsum(np.random.random(1_000_000))

In [6]: %timeit nearest_integers2D(example_array, 10.5)
9.21 µs ± 9.22 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [7]: %timeit nearest_integers_sort_fast(example_array, 10.5)
7.88 µs ± 63.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [8]: %timeit nearest_integers_sort_fast(rand_array, 10.5)
373 ms ± 3.22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

内容的提问来源于stack exchange，提问作者Diego