'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]``'

@staticmethod
def eig(C):

    def tred2(n, V, d, e):
        num_opt = False
        for j in range(n):
            d[j] = V[(n - 1)][j]
        for i in range((n - 1), 0, (-1)):
            h = 0.0
            if (not num_opt):
                scale = 0.0
                for k in range(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 range(i):
                    d[j] = V[(i - 1)][j]
                    V[i][j] = 0.0
                    V[j][i] = 0.0
            else:
                if (not num_opt):
                    for k in range(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 range(i):
                        e[j] = 0.0
                else:
                    e[:i] = 0.0
                for j in range(i):
                    f = d[j]
                    V[j][i] = f
                    g = (e[j] + (V[j][j] * f))
                    if (not num_opt):
                        for k in range((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 range(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 range(i):
                        e[j] -= (hh * d[j])
                else:
                    e[:i] -= (hh * d[:i])
                for j in range(i):
                    f = d[j]
                    g = e[j]
                    if (not num_opt):
                        for k in range(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
        for i in range((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 range((i + 1)):
                        d[k] = (V[k][(i + 1)] / h)
                else:
                    d[:(i + 1)] = (V.T[(i + 1)][:(i + 1)] / h)
                for j in range((i + 1)):
                    if (not num_opt):
                        g = 0.0
                        for k in range((i + 1)):
                            g += (V[k][(i + 1)] * V[k][j])
                        for k in range((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 range((i + 1)):
                    V[k][(i + 1)] = 0.0
            else:
                V.T[(i + 1)][:(i + 1)] = 0.0
        if (not num_opt):
            for j in range(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
    def tql2(n, d, e, V):
        num_opt = False
        if (not num_opt):
            for i in range(1, n):
                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 range(n):
            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):
                iiter = 0
                while 1:
                    iiter += 1
                    g = d[l]
                    p = ((d[(l + 1)] - g) / (2.0 * e[l]))
                    r = (((p ** 2) + 1) ** 0.5)
                    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 range((l + 2), n):
                            d[i] -= h
                    else:
                        d[(l + 2):n] -= h
                    f = (f + h)
                    p = d[m]
                    c = 1.0
                    c2 = c
                    c3 = c
                    el1 = e[(l + 1)]
                    s = 0.0
                    s2 = 0.0
                    for i in range((m - 1), (l - 1), (-1)):
                        c3 = c2
                        c2 = c
                        s2 = s
                        g = (c * e[i])
                        h = (c * p)
                        r = (((p ** 2) + (e[i] ** 2)) ** 0.5)
                        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]))))
                        if (not num_opt):
                            for k in range(n):
                                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:
                            hh = V.T[(i + 1)].copy()
                            V.T[(i + 1)] = ((s * V.T[i]) + (c * 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)
                    if (abs(e[l]) <= (eps * tst1)):
                        break
            d[l] = (d[l] + f)
            e[l] = 0.0
        if (11 < 3):
            for i in range((n - 1)):
                k = i
                p = d[i]
                for j in range((i + 1), n):
                    if (d[j] < p):
                        k = j
                        p = d[j]
                if (k != i):
                    d[k] = d[i]
                    d[i] = p
                    for j in range(n):
                        p = V[j][i]
                        V[j][i] = V[j][k]
                        V[j][k] = p
    N = len(C[0])
    if (11 < 3):
        V = np.array([x[:] for x in C])
        N = V[0].size
        d = np.zeros(N)
        e = np.zeros(N)
    else:
        V = [[x[i] for i in xrange(N)] for x in C]
        d = (N * [0.0])
        e = (N * [0.0])
    tred2(N, V, d, e)
    tql2(N, d, e, V)
    return (array(d), array(V))
