def get_points_and_weights(w_func=lambda x : np.ones(x.shape), 
    left=-1.0, right=1.0, num_points=5, n=4096):
    dx = (float(right)-left)/n
    z = np.hstack(np.linspace(left+0.5*dx, right-0.5*dx, n))
    w = dx*w_func(z)  
    (a, b) = discrete_gautschi(z, w, num_points)
    alpha = a
    beta = np.sqrt(b)
    J = np.diag(alpha)
    J += np.diag(beta, k=-1)
    J += np.diag(beta, k=1)    
    (points,v) = np.linalg.eigh(J)
    ind = points.argsort()
    points = points[ind]
    weights = v[0,:]**2 * w.sum()
    weights = weights[ind]
    return (points, weights)