Cacher is the code snippet organizer for pro developers

%%cython
import numpy as np
cdef extern from "<math.h>" nogil:
    float sqrtf(float)

def sum_tidal(int[:] ix, float[:, :] rs, int length=0):
    cdef:
        int h, i, j, k
        int n = ix.size, m = max(np.max(ix) + 1, length)
        float r1, r2, r3, r5
        float[:, :, :] res = np.zeros((m, 3, 3), dtype='f4')

    for i in range(n):
        h = ix[i]
        r2 = 0
        for j in range(3):
            r2 += rs[i, j] * rs[i, j]
        r1 = sqrtf(r2)
        r3 = r1 * r2
        r5 = r2 * r3
        for j in range(3):
            for k in range(3):
                if j == k:
                    res[h, j, k] += 3 * rs[i, j] * rs[i, k] / r5 - 1 / r3
                else:
                    res[h, j, k] += 3 * rs[i, j] * rs[i, k] / r5
    return res.base