#include <iostream>
#include <cmath>

using namespace std;

// integrated function
double func(double x, double param_k)
{
	return sqrt(1.0 - param_k * param_k * x * x) / sqrt(1.0 - x * x);
}

int main(int argc, char *argv[])
{
	long int div_num, i;
	double interval_min, interval_max, h, k;
	double xi, xip1;
	double delta_min, delta_max, delta_ave;

	// set integration interval and a parameter
	k = 1.0 / 3.0;
	interval_min = 0;
	interval_max = k;

	// set the number of division of integration interval
	div_num = 1000;

	// set stepsize
	h = (interval_max - interval_min) / (double)div_num;

	// Numerical integration
	delta_min = 0.0;
	delta_max = 0.0;
	delta_ave = 0.0;
	for(i = 0; i < div_num; i++)
	{
		xi    = interval_min + h * (double)i;
		xip1  = xi + h;
		delta_min += min(func(xi, k), func(xip1, k));
		delta_max += max(func(xi, k), func(xip1, k));
		delta_ave += (func(xi, k) + func(xip1, k)) / 2.0;
	}

	delta_min *= h;
	delta_max *= h;
	delta_ave *= h;

	// print 3 results
	cout << "Delta_min = " << delta_min << endl;
	cout << "Delta_max = " << delta_max << endl;
	cout << "Delta_ave = " << delta_ave << endl;

	return 0;
}