def tridiag_eigen(ctx, d, e, z = False):

    n = len(d)
    e[n-1] = 0
    iterlim = 2 * ctx.dps

    for l in xrange(n):
        j = 0
        while 1:
            m = l
            while 1:
                # look for a small subdiagonal element
                if m + 1 == n:
                    break
                if abs(e[m]) <= ctx.eps * (abs(d[m]) + abs(d[m + 1])):
                    break
                m = m + 1
            if m == l:
                break

            if j >= iterlim:
                raise RuntimeError("tridiag_eigen: no convergence to an eigenvalue after %d iterations" % iterlim)

            j += 1

            # form shift

            p = d[l]
            g = (d[l + 1] - p) / (2 * e[l])
            r = ctx.hypot(g, 1)

            if g < 0:
                s = g - r
            else:
                s = g + r

            g = d[m] - p + e[l] / s

            s, c, p = 1, 1, 0

            for i in xrange(m - 1, l - 1, -1):
                f = s * e[i]
                b = c * e[i]
                if abs(f) > abs(g):             # this here is a slight improvement also used in gaussq.f or acm algorithm 726.
                    c = g / f
                    r = ctx.hypot(c, 1)
                    e[i + 1] = f * r
                    s = 1 / r
                    c = c * s
                else:
                    s = f / g
                    r = ctx.hypot(s, 1)
                    e[i + 1] = g * r
                    c = 1 / r
                    s = s * c
                g = d[i + 1] - p
                r = (d[i] - g) * s + 2 * c * b
                p = s * r
                d[i + 1] = g + p
                g = c * r - b

                if not isinstance(z, bool):
                    # calculate eigenvectors
                    for w in xrange(z.rows):
                        f = z[w,i+1]
                        z[w,i+1] = s * z[w,i] + c * f
                        z[w,i  ] = c * z[w,i] - s * f

            d[l] = d[l] - p
            e[l] = g
            e[m] = 0

    for ii in xrange(1, n):
        # sort eigenvalues and eigenvectors (bubble-sort)
        i = ii - 1
        k = i
        p = d[i]
        for j in xrange(ii, n):
            if d[j] >= p:
                continue
            k = j
            p = d[k]
        if k == i:
            continue
        d[k] = d[i]
        d[i] = p

        if not isinstance(z, bool):
            for w in xrange(z.rows):
                p = z[w,i]
                z[w,i] = z[w,k]
                z[w,k] = p

########################################################################################
Remove and return a random member from a set
dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/matrices/eigen_symmetric.py/tridiag_eigen--crossover with--dataset/ETHPy150Open deldotdr/txRedis/txredis/client.py/RedisClient.spop