procedure GAULEG(X1, X2:real;var X, W:array of real; N:integer);
label 1;
const
    EPS = 0.3e-13;
var
    M,I,J:integer;
    XM,XL,P1,P2,P3,PP,Z,Z1:real;
begin
    M:=Trunc((N + 1) / 2);
    XM:=0.5 * (X2 + X1);
    XL:=0.5 * (X2 - X1);
    For I:=1 To M do
    begin
        Z:=Cos(3.141592654 * (I - 0.25) / (N + 0.5));
1:      P1:=1;
        P2:=0;
        For J:=1 To N do
        begin
            P3:=P2;
            P2:=P1;
            P1:=((2 * J - 1) * Z * P2 - (J - 1) * P3) / J;
        end;
        PP:=N * (Z * P1 - P2) / (Z * Z - 1);
        Z1:=Z;
        Z:=Z1 - P1 / PP;
        If Abs(Z - Z1) > EPS Then GoTo 1;
        X[I]:=XM - XL * Z;
        X[N + 1 - I]:=XM + XL * Z;
        W[I]:=2 * XL / ((1 - Z * Z) * PP * PP);
        W[N + 1 - I]:=W[I];
    end;
end;