def pir3(a, d):
    sidx = a.argsort()
p1 = a.searchsorted(a + d, 'right', sorter=sidx)
p2 = a.searchsorted(a - d, sorter=sidx)
return p1 - p2

def pir4(a, d):
    s = a.argsort()
y = np.empty(s.size, dtype=np.int64)
y[s] = np.arange(s.size)
a_ = a[s]
return (a_.searchsorted(a_ + d, 'right') - a_.searchsorted(a_ - d))[y]
count how many elements in a numpy array are within delta of every other element