/* calculation of integral trapezoidal method*/
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

#define MIN 0.0
#define MAX 1.2
#define EPS 1E-10

double f(double x) {return x - x*x + x*x*x - x*x*x*x + sin(x*13.)/13.;}
double fint(double x) {return x*x/2. - x*x*x/3. + x*x*x*x/4. - x*x*x*x*x/5. - cos(x*13.)/169.;}


int main() {
    double fint_exact = fint(MAX)-fint(MIN);
 	cout << "Exact: " << setprecision(18) << fint_exact << endl;

	 double h = (MAX-MIN)*0.5;
     double sum1  = f(MIN)+f(MAX);
     double sum2  = 0.;
     double sum4  = f(MIN+h);
//	 double sum = f(MIN)*0.5+f(MIN+h)+f(MAX)*0.5;
	 double a0,a1 = (sum1+sum2*2.+sum4*4.)*h/3.;
	 for(int step=1;step<=20;step++) {
         a0 = a1;
         h *= 0.5;
         sum2 += sum4;
         sum4  = 0.;
         double x = MIN+h;
         while(x<MAX-h*0.5) {
             sum4 += f(x);
             x    += (h+h);
//             sum += f(x);
//             x   += (h+h);
        }
         a1 = (sum1+sum2*2.+sum4*4.)*h/3.; //sum*h;
         cout << "Step: " << step << " h: " << setprecision(18) << h << " Numerical: " << \
         setprecision(18) << a1 << " diff: " << fabs(fint_exact-a1) << endl;

         if (fabs(a1-a0)<EPS) break;
	 }
}