function [x, w] = lp_gauss(p)
% LP_GAUSS  Gauss-Legendre nodes and weights on [-1, 1].
%
% Newton's method on the Legendre polynomial P_p, starting from the Chebyshev
% estimate of the roots; three or four iterations reach roundoff for the small
% p used here. The weights are 2 / ((1 - x^2) P_p'(x)^2).
x = cos(pi * ((1:p)' - 0.25) / (p + 0.5));
dp = ones(p, 1);
for it = 1:50
    p0 = ones(p, 1);
    p1 = x;
    for k = 2:p
        p2 = ((2*k - 1) * x .* p1 - (k - 1) * p0) / k;
        p0 = p1;
        p1 = p2;
    end
    dp = p * (x .* p1 - p0) ./ (x.^2 - 1);
    dx = p1 ./ dp;
    x = x - dx;
    if max(abs(dx)) < 1e-15
        break;
    end
end
w = 2 ./ ((1 - x.^2) .* dp.^2);
end