//=======================================================================
// ch07_composite_trapezoidal_simpson.cpp
//
// Composite trapezoidal and Simpson rules for numerical integration.
//
// The program evaluates
//
//       integral sin(x) dx,   0 <= x <= pi,
//
// for successively finer uniform grids. The exact value is 2.
//
// The example illustrates the different convergence rates of the
// composite trapezoidal and Simpson rules as the number of subintervals
// is doubled.
//
// The Simpson routine is based on code written by Alexander Godunov,
// October 2009. The companion trapezoidal routine and C++ presentation
// were prepared for the book companion website, 2026.
//=======================================================================

#include <cmath>
#include <iomanip>
#include <iostream>

using namespace std;

double f(double x);
double composite_trapezoidal(double (*fun)(double),
                             double a, double b, int n);
double composite_simpson(double (*fun)(double),
                         double a, double b, int n);

int main()
{
    const double pi = 3.14159265358979323846;
    const double a = 0.0;
    const double b = pi;
    const double exact = 2.0;

    cout << "Composite trapezoidal and Simpson rules\n";
    cout << "Integral of sin(x) from 0 to pi\n";
    cout << scientific << setprecision(7);
    cout << "Exact value          = " << setw(15) << exact << "\n\n";

    cout << "        n       Trapezoidal          Simpson\n";
    cout << "------------------------------------------------\n";

    int n = 2;

    for (int i = 0; i < 16; ++i) {
        const double trap = composite_trapezoidal(f, a, b, n);
        const double simp = composite_simpson(f, a, b, n);

        cout << setw(9) << defaultfloat << n
             << scientific << setprecision(7)
             << setw(19) << trap
             << setw(19) << simp << '\n';

        n *= 2;
    }

    return 0;
}


double f(double x)
{
    //-------------------------------------------------------------------
    // Test function for the integration example.
    //
    // To use another integrand, replace only the line below and change
    // the integration limits in main() if needed.
    //-------------------------------------------------------------------
    return sin(x);
}


double composite_trapezoidal(double (*fun)(double),
                             double a, double b, int n)
{
    //-------------------------------------------------------------------
    // Composite trapezoidal rule on [a,b].
    //
    // Input:
    //   fun - function to integrate
    //   a,b - integration limits
    //   n   - number of equal subintervals
    //
    // Output:
    //   return value - numerical approximation to the integral
    //-------------------------------------------------------------------

    if (n < 1) {
        cerr << "composite_trapezoidal: n must be positive.\n";
        return 0.0;
    }

    const double h = (b-a)/static_cast<double>(n);

    double sum = 0.5*(fun(a) + fun(b));

    for (int i = 1; i < n; ++i) {
        const double x = a + static_cast<double>(i)*h;
        sum += fun(x);
    }

    return h*sum;
}


double composite_simpson(double (*fun)(double),
                         double a, double b, int n)
{
    //-------------------------------------------------------------------
    // Composite Simpson rule on [a,b].
    //
    // Input:
    //   fun - function to integrate
    //   a,b - integration limits
    //   n   - number of equal subintervals; n must be even
    //
    // Output:
    //   return value - numerical approximation to the integral
    //-------------------------------------------------------------------

    if (n < 2 || n % 2 != 0) {
        cerr << "composite_simpson: n must be a positive even integer.\n";
        return 0.0;
    }

    const double h = (b-a)/static_cast<double>(n);

    double sum_odd = 0.0;
    for (int i = 1; i < n; i += 2) {
        const double x = a + static_cast<double>(i)*h;
        sum_odd += fun(x);
    }

    double sum_even = 0.0;
    for (int i = 2; i <= n-2; i += 2) {
        const double x = a + static_cast<double>(i)*h;
        sum_even += fun(x);
    }

    return h*(fun(a) + fun(b) + 4.0*sum_odd + 2.0*sum_even)/3.0;
}