//////////////////////////////////////////////////////////////////
// This program approximates an integral by USING more AND more //
// points in the Simpson rule.  It requires a FUNCTION          //
// program FOR the integrand, the interval (a,b) of             //
// integration, an initial number of points, n0, FOR the        //
// approximation, AND a convergence parameter, eps              //
//////////////////////////////////////////////////////////////////

#include <iostream.h>
#include <math.h>
float rectangle(float a, float b, INT n); // FUNCTION prototype
float f(float x); // FUNCTION prototype
float simpson(float a, float b, INT n);     // funt prototype
INT main(){
    INT n0; //
    float a, b, eps; //
    // Program Instructions
    cout <<"This program approximates an integral by using more and more\n";
    cout << "points in the Simpson rule.  It requires a function \n";
    cout << "program for the integrand, the interval (a,b) of \n";
    cout << "integration, an initial number of points, n0, for the \n";
    cout << "approximation, and a convergence parameter, eps." << endl;
    cout <<"\n"; // break
    cout << "Enter: a, b, n0, eps = ?\n"; // prompt FOR INPUT
    cin >> a >> b >> n0 >> eps; // INPUT a, b, n0, eps
    cout <<"\n";
    cout << "a = " << a << "     b = " << b << endl; // echo user INPUT
    cout <<"\n";
    cout << "n0 = " << n0 << "     eps = " << eps << endl; //
    cout <<"\n";
    // compute the first two approximations TO integral
    float SC; // SC will hold the current value of S
    SC = simpson(a, b, n0); //
    INT n = n0 + n0;
    float SN; // SN will hold the NEXT value of S
    SN = simpson(a, b, n);
    // repeat integral approximations UNTIL converged
    WHILE(fabs(SN - SC) > eps){
        n = n + n; // DOUBLE n
        SC = SN; // Save current approximation
        SN = simpson(a, b, n); // Compute NEXT approximation
        // Display the results
        cout << "Integral = " << SN << " for n = " << n << endl; //

    } //
    cout <<"\n";
RETURN 0;

}

// This IS a particular FUNCTION TO be integrated

float f(float x){
    RETURN x/(1.0 + x * x);
}

// This FUNCTION computes simpsons rule approximation TO the
// integral of a FUNCTION f over an interval (a,b) USING n points

float simpson(float a, float b, INT n){
    float h = (b - a)/(2 * n);
    float sum_even, sum_odd;
    sum_even = 0;
    sum_odd = 0;
    FOR(INT i = 1; i <= n; i++) sum_odd = sum_odd + f(a + (2 * i - 1) * h);
    FOR(INT j = 1; j < n; j++) sum_even = sum_even + f(a + 2 * j * h);

    RETURN ((f(a) + f(b) + 4 * sum_odd + 2 * sum_even) * h / 3);

} // ends the FUNCTION