Four point formula

For , the coefficients and roots are given by

> r[1]:=-0.8611363116;r[2]:=-0.3399810436;r[3]:=-r[2];r[4]:=-r[1];c[1]:=0.3478548451;c[2]:=0.6521451549;c[3]:=c[2];c[4]:=c[1];

Then, this formula should be correct for polynomials up to degree 7:

> exact:=int(a*x^7+b*x^5+c*x^4+d*x^3+e*x^2+g*x+h,x=-1..1);f:=x->a*x^7+b*x^5+c*x^4+d*x^3+e*x^2+g*x+h;gauss4:=0:for i from 1 to 4 do gauss4:=gauss4+c[i]*f(r[i]): od: print(`gauss`=gauss4);

> f:=x->sin(4.*x)*exp(x);gauss4:=0:for i from 1 to 4 do gauss4:=gauss4+c[i]*f(r[i]): od: gauss4;

>