Função para Gauss-Legendre para Gretl

function [x,w]=gauss_legendre(n)

    // [x,w]=gauss_legendre(n)

    // Calcula as raízes "x" e os pesos "w" dos polinômios de Gauss-Legendre de

    // ordem "n".

    // (C) Rudnei Dias da Cunha 2008

    tol = number_properties('eps')

    a = -1.0;

    b = 1.0;

    pi = acos(-1.0);

    xl = (b-a)/2;

    xm = (b+a)/2;

    m = floor((n+1)/2);

    for i=1:m

        z1 = 0

        z = cos(pi*(i-0.25)/(n+0.25));

        while abs(z-z1)>tol

            p1 = 1.0;

            p2 = 0.0;

            for j=1:n

                p3 = p2;

                p2 = p1;

                p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/real(j);

            end

            pp = n*(z*p1-p2)/(z*z-1.0);

            z1 = z;

            z = z1-p1/pp;

        end

        x(i) = xm+xl*z;

        x(n+1-i) = xm-xl*z;

        w(i) = 2.0*xl/((1.0-z^2)*pp^2);

        w(n+1-i) = w(i);

    end

endfunction