code: 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)