docs: 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)))