/* 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 sum = f(MIN)*0.5+f(MIN+h)+f(MAX)*0.5;
	 double a0,a1 = sum*h;
	 for(int step=1;step<=20;step++) {
         a0 = a1;
         h *= 0.5;
         double x = MIN+h;
         while(x<MAX-h*0.5) {
             sum += f(x);
             x   += (h+h);
        }
         a1 = sum*h;
         cout << "Step: " << step << " h: " << setprecision(18) << h << " Numerical: " << \
         setprecision(18) << a1 << " diff: " << fabs(fint_exact-a1) << endl;

//         printf("Step %2d => h: %.16f, Numerical: %.16f, diff: %.16f\n", step,h,a1,fabs(fint_exact-a1));
         if (fabs(a1-a0)<EPS) break;
	 }
}