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

// f(x) = 1 + x^2 + x^3 + ... + x^power
double f(double x, int power) {
    double val = 0.;
    for (int i=0; i<=power; i++) val += pow(x,i);
    return val;
}
double fint(double x, int power) {
    double val = 0.;
    for (int i=0; i<=power; i++) val += pow(x,i+1)/(double)(i+1);
    return val;
}

#define MIN	0.0
#define MAX	1.0

#define NPOINTS 5
const double GQ_weights[NPOINTS][2] = {
    0.5688888888888889,	+0.0000000000000000,
    0.4786286704993665,	-0.5384693101056831,
    0.4786286704993665,	+0.5384693101056831,
    0.2369268850561891,	-0.9061798459386640,
    0.2369268850561891,	+0.9061798459386640};

int main()
{
    for (int test_power = 1; test_power<=15; test_power++) {

        double fint_exact = fint(MAX,test_power)-fint(MIN,test_power);
        double area = 0.;

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

        cout << "Power: " << test_power << " Exact: " << setprecision(18) << fint_exact << " Numerical: " << \
        setprecision(18) << area << " diff: " << fabs(fint_exact-area) << endl;

//        printf("Power: %2d, Exact: %.16f, Numerical: %.16f, diff: %.16f\n",test_power,fint_exact,area,fabs(fint_exact-area));
    }
}