def _separable_series3(h, N=1, verbose=False):
    """ finds separable approximations to the 3d kernel h
    returns res = (hx,hy,hz)[N]
    s.t. h \approx sum_i einsum("i,j,k",res[i,0],res[i,1],res[i,2])

    FIXME: This is just a naive and slow first try!
    """

    hx, hy, hz = [], [], []
    res = h.copy()
    for i in range(N):
        _hx, _hy, _hz, P = _splitrank3(res, verbose=verbose)
        res -= P
        hx.append(_hx)
        hy.append(_hy)
        hz.append(_hz)
    return np.array(list(zip(hx, hy, hz)))
finds separable approximations to the 3d kernel h
    returns res = (hx,hy,hz)[N]
    s.t. h \approx sum_i einsum("i,j,k",res[i,0],res[i,1],res[i,2])

    FIXME: This is just a naive and slow first try!