#include <math.h>
#include <stdio.h>
#include <iostream>
using namespace std;

double MyFunc(double x) // Ôóíêö³ÿ, â³ä ÿêî¿ òðåáà îá÷èñëèòè ³íòåãðàë
{
    return exp(x) * log(1 + x * x);
}

double Simpson(double (*func)(double x), double a, double b, double eps, double *h_)
{
    double h, x, I, I1, I2, I3;
    I2 = 1;
    h = b - a;
    I = (*func)(a) + (*func)(b);

    do
    {
            I3 = I2;
            h /= 2;
            I1 = 0;
            x = a +h;
            do
            {
                I1 += 2 * (*func)(x);
                x += 2*h;
            }
            while(x < b);
            I += I1;
            I2 = (I + I1) * h / 3;
            x = fabs(I3 - I2) / 15;
    }
    while(x > eps);
    *h_ = h;
    return I2;
}
double Rect(double (*func)(double x), double a, double b, double eps, double *h_)
{
    double x, h, I, I1;
    h = b - a;
    I = (*func)((a + b) / 2) * h;
    x = a;
    do
    {
        x = a;
        I1 = I;
        I = 0;
        h /= 2;
        do
        {
            x += 2 * h;
            I += (*func)(x) * (2 * h);
        }
        while(x < b);
        x = fabs(I1 - I) / 3;
    }
    while(x > eps);
    *h_ = h;
    return I;
}
double Trap(double (*func)(double x), double a, double b, double eps, double *h_)
{
    double x, h, I, I1, I3;
    h = b - a;
    I = 1;
    do
    {

        I3 = I;
        I = ((*func)(a) + (*func)(b)) / 2;
        h /= 2;
        x = a + h;
        do
        {
            I += (*func)(x);
            x += h;
        }
        while(x < b);
        I *= h;
        x = fabs(I3 - I) / 3;
    }
    while(x > eps);
    *h_ = h;
    return I;
}

int main(int argc, char* argv[])
{

    double x0 = 0;
    double x1 = 1;
    int k = 3;
    double int_S, int_R, int_T, eps = 0.1;
    double h_S = 1, h_R = 1, h_T = 1;
    printf("Simpson     h_S        Rect        h_R        Trap        h_T        eps\n");
    for(int i = 0; i < k; i++)
    {
        int_S = Simpson(&MyFunc, x0, x1, eps, &h_S);
        int_R = Rect(&MyFunc, x0, x1, eps, &h_R);
        int_T = Trap(&MyFunc, x0, x1, eps, &h_T);
        printf("%+.6f   %.6f   %+.6f   %.6f   %+.6f   %.6f   %.3f\n", int_S, h_S, int_R, h_R, int_T, h_T, eps);
        eps /= 10;
    }
    
    system("PAUSE");
    return 0;
}