#include <math.h>
#include <stdio.h>
#include <stdlib.h>

// weights for Simpson's Rule
const double weightA = 0.333333333333333;
const double weightB = 1.333333333333333;

#define Nmax 40
#define Nmax2 Nmax*2

typedef double (*PFI) (double);

typedef struct {
    double x[Nmax2], f[Nmax2];      // f is cached value of f(x)
    double E[Nmax], Int[Nmax];
    double oldI, finI, newI;
    int n, np;
    PFI pn;
} Simpson;

#define _x      this->x
#define _f      this->f
#define _n      this->n
#define _np     this->np
#define _E      this->E
#define _Int    this->Int
#define _finI   this->finI
#define _oldI   this->oldI
#define _newI   this->newI
#define _pn     this->pn


double simpsonStep(Simpson* this, int a, int b, int c) {
    double dx = _x[c] - _x[b];
    return dx * (weightA * (_f[a] + _f[c]) + weightB * _f[c]);
}


void initialize(Simpson* this, double a, double b, double eps, PFI pn) {
    _n = 2;
    _E[0] = eps;
    _x[0] = a;
    _x[2] = b;
    _x[1] = 0.5 * (a + b);
    _f[0] = pn(a);
    _f[2] = pn(b);
    _f[1] = pn(_x[1]);
    _Int[0] = simpsonStep(this, 0, 1, 2);
    _finI = 0.0;
    _pn = pn;
}


Simpson* subdivide(Simpson* this) {
    _n += 2;

    if ((_n) > Nmax - 1) {	/* check n */
        printf(" Too many subdivisions");
        abort();
    }
    _np = _n / 2 - 2;
    _oldI = _Int[_np];
    _E[_np + 1] = 0.5 * _E[_np];

    return this;
}


Simpson* move(Simpson* this, int to, int from) {
    _x[to] = _x[from];
    _f[to] = _f[from];
    return this;
}


Simpson* moveDown(Simpson* this) {
           move(this,   _n, _n-2);
    return move(this, _n-2, _n-3);
}


Simpson* averageSurrounding(Simpson* this, int point) {
    _x[point] = 0.5 * (_x[point-1] + _x[point+1]);
    _f[point] = _pn(_x[point]);
    return this;
}


Simpson* newPoints(Simpson* this) {
    averageSurrounding(this, _n-1);
    averageSurrounding(this, _n-3);
    _Int[_np + 1] = simpsonStep(this, _n-2, _n-1, _n);
    _Int[_np]     = simpsonStep(this, _n-4, _n-3, _n-2);
    return this;
}


int converged(Simpson* this) {
    double deltaI;

    _newI = _Int[_np] + _Int[_np + 1];
    deltaI = _newI - _oldI;
    if (fabs(deltaI) < _E[_np]) {
        _finI = deltaI / 15 + _newI + _finI;
        _n -= 4;
    }
    return _n <= 0;
}


double integral(double a, double b, double eps, PFI pn) {
    Simpson i;

    initialize(&i, a, b, eps, pn);
    while (!converged(newPoints(moveDown(subdivide(&i)))));
    return i.finI;
}


int main(int argc, char** argv)
{
    double a, b, eps;

    printf("What are a, b, and epsilon? ");
    scanf(" %lf%lf%lf", &a, &b, &eps);
    printf(" %.10le %.10le %.10le\n", a, b, eps);

    printf(" Integral with sqrt is %.10lf\n", integral(a, b, eps, sqrt));
    printf(" Integral with exp is %.10lf\n",  integral(a, b, eps, exp));

    return 0;
}