#include "vs.h"
int main() {
	// A. Bode's Integration Formula--exact up to 5-th order
	double x = 0.0,
   	w[17] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                              64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                              64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                              64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};

   for(int i = 0; i < 11; i++) {
   	Quadrature q1(w, 0.0, x, 17), q2(w, x, 1.0, 17);
      double c = 1.0/2.0, B = -1.0/6.0;
      H0 z1(q1), z2(q2);
      H0 A = (1.0 + z2) / 2.0,
      	a = 1.0/2.0+1.0/2.0*z1.pow(2)-1.0/6.0*z1.pow(3),
      	b = -1.0 - 1.0/2.0*z1.pow(2),
      	integrand1 = a + b*(1.0-x)+ c*pow((1.0-x), 2),
         integrand2 = A*pow(x, 2) + B*pow(x, 3);
      C0 integal_1, integal_2;
      if(i != 0) integal_1 &= integrand1 | J(x / 16.0); else integal_1 &= C0(0.0);
      if(i != 10) integal_2 &= integrand2 | J((1.0-x) / 16.0); else integal_2 &= C0(0.0);
      double u = (double) (integal_1 + integal_2);
   	cout << "u(" << x << "): " << u << endl;
   	if(i != 10) x += 0.1;
   }
	return 0;
}