(* 
Notebook calculates abcissas and weights for the Gaussian quadrature
for the integral Integrate[F[x]/(1+Exp[x]),{x,1,Infinity}].
Pre-calculated related orthogonal polynomials are used.
*)

num=32 (* Ngauss+1 *)

poly=Table[Table[ToExpression["a"<>ToString[i]<>ToString[j]]*x^i,{i,0,j}]/.List->Plus,{j,0,num}];

ortoDAT=Import["https://ribes.if.uj.edu.pl/~va00/Orto"<>ToString[87]<>".m"];
(* 
coefficients of the highest order ortogonal polynomial found until now; 87, 2010-07-20
*)

orto=poly/.ortoDAT; (* list of the orthogonal polynomials *)

Xtbl=x/.FindInstance[orto[[-2]]==0,x,num+2]//Sort//Reverse; (* abcissas *)

Wtbl1=Table[1/Sum[(orto[[k]]/.x->Xtbl[[ii]])^2,{k,1,num-1}],{ii,1,num-1}]; 
(* weights,formula 1 *)

Wtbl2=Table[-Last[CoefficientList[orto[[-1]],x]]/Last[CoefficientList[orto[[-2]],x]]/(orto[[-1]]/.x->Xtbl[[ii]])/D[orto[[-2]],x]/.x->Xtbl[[ii]],{ii,1,num-1}];
(* weights,formula 2 *)

Wtbl1-Wtbl2

quadrature=Sum[Wtbl1[[k]]F[Xtbl[[k]]],{k,1,num-1}];

testFUN1=#^7-#^4+3  #^3-12& (* polynomial *)

testFUN2=#^3 Sqrt[1+#]& (* non-polynomial *)

result1=quadrature /. F->testFUN1

directresult1=NIntegrate[testFUN1[x]/(1+Exp[x]),{x,1,Infinity},WorkingPrecision->128]

abserr1=result1-directresult1

relerr1=1-result1/directresult1

result2=quadrature /. F->testFUN2

directresult2=NIntegrate[testFUN2[x]/(1+Exp[x]),{x,1,Infinity},WorkingPrecision->128]

abserr2=result2-directresult2//N

relerr2=1-result2/directresult2//N