#include <cassert>
#include <cmath>

// Approximate the following integral using the composite Simpson's rule
//  - this finds the ell-two error between the two functions
//
//                ---------------------
//               /  b
//              /  /              2
//             /   | (f(x) - g(x))  dx
//         \  /    / 
//          \/    a
//                                    4
// The composite Simpson's rule is O(h ), therefore, if we have two
// successive approximations
//      prev_approx          // using h
//      curr_approx          // using h/2
// then the error will be reduced by a factor of 16
//
//                curr_approx - prev_approx
//   - Therefore, ------------------------- will be a good
//                            15
//     approximation of the error of 'curr_approx'
//
//   - Therefore, local extrapolation allows us to find an
//     even better approximation of the integral using the
//     formula
//              16 curr_approx - prev_approx
//              ----------------------------
//                           15
//
// We will iteratively find better and better approximations,
// each time halving 'h' until 1.5 times the approximation of the
// error of the most current approximation is less than the
// required amount. That is, we are building in a safety factor
// of 150%.
//
// Algorithm
// =========
//
// Recall that if we use the composite Simpson's rule, once with 'h'
// and then again with 'h/2', then the second evaluation will require
// a recalculation of the difference squared at each second point:
//
//
//         n = 3                   ****o*******o
//                       ******o***             ***
//               ******o*                          ****o
//         ****o*
//     o***
//
//     1       4       2       4       2       4       1
//     |-----------------------------------------------|
//     a                                               b
//
//         n = 6                   o***o***o***o
//                       **o***o***             ***
//               **o***o*                          o***o
//         o***o*
//     o***
//
//     1   4   2   4   2   4   2   4   2   4   2   4   1
//     |-----------------------------------------------|
//     a                                               b
//
//
// Note that those points with coefficient '4' at one step
// are absorbed into those points that have coefficent '2' at
// the next step. Thus, we will calculate separately the sum
// of coefficient-'4' terms and the coefficient-'2' terms.
// Then, for the next step, we will absorb the previous 
// coefficient-'4' terms into the coefficient-'2' terms
// and calculate the new sum of coefficieint-'4' terms.

template <typename T>
static T l2sum( T f( T x ), T g( T x ), T a, T h, unsigned int n, T offset );

template <typename T>
T l2error(
    T f( T x ),
    T g( T x ),
    T a,
    T b,
    unsigned int n,
    T eps_abs,
    unsigned int max_iterations
) {
    assert( b > a );

    T const SAFETY_FACTOR{ 1.5 };     // 150% safety factor

    T width{ b - a };
    T h{ width/n };

    // End-points
    T delta_a{ f( a ) - g( a ) };
    T delta_b{ f( b ) - g( b ) };
    T end_points{ delta_a*delta_a + delta_b*delta_b };

    T coeff_2_sum{ l2sum( f, g, a, h, n - 1, 1.0 ) };
    T coeff_4_sum{ l2sum( f, g, a, h, n,     0.5 ) };

    T approx_current{ (end_points + 2.0*coeff_4_sum + coeff_2_sum)*h/6.0 };

    for ( unsigned int k{1}; k <= max_iterations; ++k ) {
        // Store the most recent approximation
        T approx_previous{ approx_current };
        coeff_2_sum += coeff_4_sum;

        // Calculate the next interior points
	n *= 2;
	h /= 2.0;
        coeff_4_sum = l2sum( f, g, a, h, n, 0.5 );
        approx_current = (end_points + 2.0*coeff_4_sum + coeff_2_sum)*h/6.0;;

        std::clog << " " << k << ": " << approx_current << std::endl;

	T error{ std::abs( (approx_current - approx_previous)/15.0 ) };

	if ( SAFETY_FACTOR*error < eps_abs ) {
            // Use local extrapolation
	    return std::sqrt( (16.0*approx_current - approx_previous)/15.0 );
	}
    }

    assert( false );
    return 0.0;
}

// Calculate
//
//        n - 1
//        -----
//       \                                                      2
//  2.0   \      ( f(a + (k + offset)h) - g(a + (k + offset)h) )  dx
//        /    
//       /_____
//        k = 0
//
// where 'offset' is either 0.5 or 1.0.

template <typename T>
static T l2sum( T f( T x ), T g( T x ), T a, T h, unsigned int n, T offset ) {
    T sum{ 0.0 };

    for ( unsigned int k{0}; k < n; ++k ) {
        T x{ a + (k + offset)*h };
        T delta{ f( x ) - g( x ) };
	sum += delta*delta;
    }

    return 2.0*sum;
}