@defun
def gauss_quadrature(ctx, n, qtype = "legendre", alpha = 0, beta = 0):

    d = ctx.zeros(n, 1)
    e = ctx.zeros(n, 1)
    z = ctx.zeros(1, n)

    z[0,0] = 1

    if qtype == "legendre":
        # legendre on the range -1 +1 , abramowitz, table 25.4, p.916
        w = 2
        for i in xrange(n):
            j = i + 1
            e[i] = ctx.sqrt(j * j / (4 * j * j - ctx.mpf(1)))
    elif qtype == "legendre01":
        # legendre shifted to 0 1        , abramowitz, table 25.8, p.921
        w = 1
        for i in xrange(n):
            d[i] = 1 / ctx.mpf(2)
            j = i + 1
            e[i] = ctx.sqrt(j * j / (16 * j * j - ctx.mpf(4)))
    elif qtype == "hermite":
        # hermite on the range -inf +inf , abramowitz, table 25.10,p.924
        w = ctx.sqrt(ctx.pi)
        for i in xrange(n):
            j = i + 1
            e[i] = ctx.sqrt(j / ctx.mpf(2))
    elif qtype == "laguerre":
        # laguerre on the range 0 +inf , abramowitz, table 25.9, p. 923
        w = 1
        for i in xrange(n):
            j = i + 1
            d[i] = 2 * j - 1
            e[i] = j
    elif qtype=="chebyshev1":
        # chebyshev polynimials of the first kind
        w = ctx.pi
        for i in xrange(n):
            e[i] = 1 / ctx.mpf(2)
        e[0] = ctx.sqrt(1 / ctx.mpf(2))
    elif qtype == "chebyshev2":
        # chebyshev polynimials of the second kind
        w = ctx.pi / 2
        for i in xrange(n):
            e[i] = 1 / ctx.mpf(2)
    elif qtype == "glaguerre":
        # generalized laguerre on the range 0 +inf
        w = ctx.gamma(1 + alpha)
        for i in xrange(n):
            j = i + 1
            d[i] = 2 * j - 1 + alpha
            e[i] = ctx.sqrt(j * (j + alpha))
    elif qtype == "jacobi":
        # jacobi polynomials
        alpha = ctx.mpf(alpha)
        beta = ctx.mpf(beta)
        ab = alpha + beta
        abi = ab + 2
        w = (2**(ab+1)) * ctx.gamma(alpha + 1) * ctx.gamma(beta + 1) / ctx.gamma(abi)
        d[0] = (beta - alpha) / abi
        e[0] = ctx.sqrt(4 * (1 + alpha) * (1 + beta) / ((abi + 1) * (abi * abi)))
        a2b2 = beta * beta - alpha * alpha
        for i in xrange(1, n):
            j = i + 1
            abi = 2 * j + ab
            d[i] = a2b2 / ((abi - 2) * abi)
            e[i] = ctx.sqrt(4 * j * (j + alpha) * (j + beta) * (j + ab) / ((abi * abi - 1) * abi * abi))
    elif isinstance(qtype, str):
        raise ValueError("unknown quadrature rule \"%s\"" % qtype)
    elif not isinstance(qtype, str):
        w = qtype(d, e)
    else:
        assert 0

    tridiag_eigen(ctx, d, e, z)

    for i in xrange(len(z)):
        z[i] *= z[i]

    z = z.transpose()
    return (d, w * z)

##################################################################################################
##################################################################################################
##################################################################################################
This routine calulates gaussian quadrature rules for different
dataset/ETHPy150Open fredrik-johansson/mpmath/mpmath/matrices/eigen_symmetric.py/gauss_quadrature