%%
x = load('gausspoints.dat');
w = load('gaussweights.dat');

%%
[qx,qw] = qs(3, x, w);

poly = [3 2 1];
ep = polyval(poly, qx);
integ = qw'*ep;

ref = integral(@(arg)(polyval(poly, arg)), -1, 1);
fprintf('%g %g %g\n', integ, ref, abs(integ-ref)/ref);

%%
poly = [1 1 1 1 3 2 1];
ep = polyval(poly, qx);
integ = qw'*ep;

ref = integral(@(arg)(polyval(poly, arg)), -1, 1);
fprintf('%g %g %g\n', integ, ref, abs(integ-ref)/ref);

%%
poly = [1 1 1 1 1 1 3 2 1];
ep = polyval(poly, qx);
integ = qw'*ep;

ref = integral(@(arg)(polyval(poly, arg)), -1, 1);
fprintf('%g %g %g\n', integ, ref, abs(integ-ref)/ref);

%%
[qx,qw] = qs(4, x, w);
poly = [1 1 1 1 1 1 3 2 1];
ep = polyval(poly, qx);
integ = qw'*ep;

ref = integral(@(arg)(polyval(poly, arg)), -1, 1);
fprintf('%g %g %g\n', integ, ref, abs(integ-ref)/ref);