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

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.;}

#define MIN	0.0
#define MAX	1.2

#define NPOINTS 21
const double GQ_weights[NPOINTS][2] = {
    0.1460811336496904, +0.0000000000000000,
    0.1445244039899700, -0.1455618541608951,
    0.1445244039899700, +0.1455618541608951,
    0.1398873947910731, -0.2880213168024011,
    0.1398873947910731, +0.2880213168024011,
    0.1322689386333375, -0.4243421202074388,
    0.1322689386333375, +0.4243421202074388,
    0.1218314160537285, -0.5516188358872198,
    0.1218314160537285, +0.5516188358872198,
    0.1087972991671484, -0.6671388041974123,
    0.1087972991671484, +0.6671388041974123,
    0.0934444234560339, -0.7684399634756779,
    0.0934444234560339, +0.7684399634756779,
    0.0761001136283793, -0.8533633645833173,
    0.0761001136283793, +0.8533633645833173,
    0.0571344254268572, -0.9200993341504008,
    0.0571344254268572, +0.9200993341504008,
    0.0369537897708525, -0.9672268385663063,
    0.0369537897708525, +0.9672268385663063,
    0.0160172282577743, -0.9937521706203895,
    0.0160172282577743, +0.9937521706203895};

int main()
{
    double fint_exact = fint(MAX)-fint(MIN);

    double area = 0.;

    for (int i=0; i<NPOINTS; i++) {
        double x = ((MAX-MIN)*GQ_weights[i][1] + (MAX+MIN))/2.;
        area += f(x)*GQ_weights[i][0];
    }
    area *= (MAX-MIN)/2.;

    cout << "Exact: " << setprecision(18) <<fint_exact << " Numerical: " << \
    setprecision(18) << area << " diff: " << fabs(fint_exact-area) << endl;
//    printf("Exact: %.16f, Numerical: %.16f, diff: %.16f\n",fint_exact,area,fabs(fint_exact-area));
}