docs: def eig(C):
          """eigendecomposition of a symmetric matrix, much slower than
          `numpy.linalg.eigh`, return ``(EVals, Basis)``, the eigenvalues
          and an orthonormal basis of the corresponding eigenvectors, where
  
              ``Basis[i]``
                  the i-th row of ``Basis``
              columns of ``Basis``, ``[Basis[j][i] for j in range(len(Basis))]``
                  the i-th eigenvector with eigenvalue ``EVals[i]``
  
          """
  
      # class eig(object):
      #     def __call__(self, C):
  
      # Householder transformation of a symmetric matrix V into tridiagonal form.
          # -> n             : dimension
          # -> V             : symmetric nxn-matrix
          # <- V             : orthogonal transformation matrix:
          #                    tridiag matrix == V * V_in * V^t
          # <- d             : diagonal
          # <- e[0..n-1]     : off diagonal (elements 1..n-1)
  
          # Symmetric tridiagonal QL algorithm, iterative
          # Computes the eigensystem from a tridiagonal matrix in roughtly 3N^3 operations
          # -> n     : Dimension.
          # -> d     : Diagonale of tridiagonal matrix.
          # -> e[1..n-1] : off-diagonal, output from Householder
          # -> V     : matrix output von Householder
          # <- d     : eigenvalues
          # <- e     : garbage?
          # <- V     : basis of eigenvectors, according to d
  
  
          #  tred2(N, B, diagD, offdiag); B=C on input
          #  tql2(N, diagD, offdiag, B);
  
          #  private void tred2 (int n, double V[][], double d[], double e[]) {
          def tred2 (n, V, d, e):
              #  This is derived from the Algol procedures tred2 by
              #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
              #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
              #  Fortran subroutine in EISPACK.
  
              num_opt = False  # factor 1.5 in 30-D
  
              for j in xrange(n):
                  d[j] = V[n - 1][j]  # d is output argument
  
              # Householder reduction to tridiagonal form.
  
              for i in xrange(n - 1, 0, -1):
                  # Scale to avoid under/overflow.
                  h = 0.0
                  if not num_opt:
                      scale = 0.0
                      for k in xrange(i):
                          scale = scale + abs(d[k])
                  else:
                      scale = sum(abs(d[0:i]))
  
                  if scale == 0.0:
                      e[i] = d[i - 1]
                      for j in xrange(i):
                          d[j] = V[i - 1][j]
                          V[i][j] = 0.0
                          V[j][i] = 0.0
                  else:
  
                      # Generate Householder vector.
                      if not num_opt:
                          for k in xrange(i):
                              d[k] /= scale
                              h += d[k] * d[k]
                      else:
                          d[:i] /= scale
                          h = np.dot(d[:i], d[:i])
  
                      f = d[i - 1]
                      g = h**0.5
  
                      if f > 0:
                          g = -g
  
                      e[i] = scale * g
                      h = h - f * g
                      d[i - 1] = f - g
                      if not num_opt:
                          for j in xrange(i):
                              e[j] = 0.0
                      else:
                          e[:i] = 0.0
  
                      # Apply similarity transformation to remaining columns.
  
                      for j in xrange(i):
                          f = d[j]
                          V[j][i] = f
                          g = e[j] + V[j][j] * f
                          if not num_opt:
                              for k in xrange(j + 1, i):
                                  g += V[k][j] * d[k]
                                  e[k] += V[k][j] * f
                              e[j] = g
                          else:
                              e[j + 1:i] += V.T[j][j + 1:i] * f
                              e[j] = g + np.dot(V.T[j][j + 1:i], d[j + 1:i])
  
                      f = 0.0
                      if not num_opt:
                          for j in xrange(i):
                              e[j] /= h
                              f += e[j] * d[j]
                      else:
                          e[:i] /= h
                          f += np.dot(e[:i], d[:i])
  
                      hh = f / (h + h)
                      if not num_opt:
                          for j in xrange(i):
                              e[j] -= hh * d[j]
                      else:
                          e[:i] -= hh * d[:i]
  
                      for j in xrange(i):
                          f = d[j]
                          g = e[j]
                          if not num_opt:
                              for k in xrange(j, i):
                                  V[k][j] -= (f * e[k] + g * d[k])
                          else:
                              V.T[j][j:i] -= (f * e[j:i] + g * d[j:i])
  
                          d[j] = V[i - 1][j]
                          V[i][j] = 0.0
  
                  d[i] = h
              # end for i--
  
              # Accumulate transformations.
  
              for i in xrange(n - 1):
                  V[n - 1][i] = V[i][i]
                  V[i][i] = 1.0
                  h = d[i + 1]
                  if h != 0.0:
                      if not num_opt:
                          for k in xrange(i + 1):
                              d[k] = V[k][i + 1] / h
                      else:
                          d[:i + 1] = V.T[i + 1][:i + 1] / h
  
                      for j in xrange(i + 1):
                          if not num_opt:
                              g = 0.0
                              for k in xrange(i + 1):
                                  g += V[k][i + 1] * V[k][j]
                              for k in xrange(i + 1):
                                  V[k][j] -= g * d[k]
                          else:
                              g = np.dot(V.T[i + 1][0:i + 1], V.T[j][0:i + 1])
                              V.T[j][:i + 1] -= g * d[:i + 1]
  
                  if not num_opt:
                      for k in xrange(i + 1):
                          V[k][i + 1] = 0.0
                  else:
                      V.T[i + 1][:i + 1] = 0.0
  
  
              if not num_opt:
                  for j in xrange(n):
                      d[j] = V[n - 1][j]
                      V[n - 1][j] = 0.0
              else:
                  d[:n] = V[n - 1][:n]
                  V[n - 1][:n] = 0.0
  
              V[n - 1][n - 1] = 1.0
              e[0] = 0.0
  
  
          # Symmetric tridiagonal QL algorithm, taken from JAMA package.
          # private void tql2 (int n, double d[], double e[], double V[][]) {
          # needs roughly 3N^3 operations
          def tql2 (n, d, e, V):
  
              #  This is derived from the Algol procedures tql2, by
              #  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
              #  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
              #  Fortran subroutine in EISPACK.
  
              num_opt = False  # using vectors from numpy makes it faster
  
              if not num_opt:
                  for i in xrange(1, n):  # (int i = 1; i < n; i++):
                      e[i - 1] = e[i]
              else:
                  e[0:n - 1] = e[1:n]
              e[n - 1] = 0.0
  
              f = 0.0
              tst1 = 0.0
              eps = 2.0**-52.0
              for l in xrange(n):  # (int l = 0; l < n; l++) {
  
                  # Find small subdiagonal element
  
                  tst1 = max(tst1, abs(d[l]) + abs(e[l]))
                  m = l
                  while m < n:
                      if abs(e[m]) <= eps * tst1:
                          break
                      m += 1
  
                  # If m == l, d[l] is an eigenvalue,
                  # otherwise, iterate.
  
                  if m > l:
                      iiter = 0
                      while 1:  # do {
                          iiter += 1  # (Could check iteration count here.)
  
                          # Compute implicit shift
  
                          g = d[l]
                          p = (d[l + 1] - g) / (2.0 * e[l])
                          r = (p**2 + 1)**0.5  # hypot(p,1.0)
                          if p < 0:
                              r = -r
  
                          d[l] = e[l] / (p + r)
                          d[l + 1] = e[l] * (p + r)
                          dl1 = d[l + 1]
                          h = g - d[l]
                          if not num_opt:
                              for i in xrange(l + 2, n):
                                  d[i] -= h
                          else:
                              d[l + 2:n] -= h
  
                          f = f + h
  
                          # Implicit QL transformation.
  
                          p = d[m]
                          c = 1.0
                          c2 = c
                          c3 = c
                          el1 = e[l + 1]
                          s = 0.0
                          s2 = 0.0
  
                          # hh = V.T[0].copy()  # only with num_opt
                          for i in xrange(m - 1, l - 1, -1):  # (int i = m-1; i >= l; i--) {
                              c3 = c2
                              c2 = c
                              s2 = s
                              g = c * e[i]
                              h = c * p
                              r = (p**2 + e[i]**2)**0.5  # hypot(p,e[i])
                              e[i + 1] = s * r
                              s = e[i] / r
                              c = p / r
                              p = c * d[i] - s * g
                              d[i + 1] = h + s * (c * g + s * d[i])
  
                              # Accumulate transformation.
  
                              if not num_opt:  # overall factor 3 in 30-D
                                  for k in xrange(n):  # (int k = 0; k < n; k++) {
                                      h = V[k][i + 1]
                                      V[k][i + 1] = s * V[k][i] + c * h
                                      V[k][i] = c * V[k][i] - s * h
                              else:  # about 20% faster in 10-D
                                  hh = V.T[i + 1].copy()
                                  # hh[:] = V.T[i+1][:]
                                  V.T[i + 1] = s * V.T[i] + c * hh
                                  V.T[i] = c * V.T[i] - s * hh
                                  # V.T[i] *= c
                                  # V.T[i] -= s * hh
  
                          p = -s * s2 * c3 * el1 * e[l] / dl1
                          e[l] = s * p
                          d[l] = c * p
  
                          # Check for convergence.
                          if abs(e[l]) <= eps * tst1:
                              break
                      # } while (Math.abs(e[l]) > eps*tst1);
  
                  d[l] = d[l] + f
                  e[l] = 0.0
  
  
              # Sort eigenvalues and corresponding vectors.
          # tql2
  
          N = len(C[0])
          if 1 < 3:
              V = [[x[i] for i in xrange(N)] for x in C]  # copy each "row"
              d = N * [0.]
              e = N * [0.]
  
          tred2(N, V, d, e)
          tql2(N, d, e, V)
          return (array(d), array(V))