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

// Kahan's compenseated summation
// compile this w/o optinization ( g++ -O0 test.cc - o test )
// program takes an integer parameter : the grid step .

template <class MyFloat>
MyFloat func(MyFloat x,MyFloat y )
{
	MyFloat r2=x*x+y*y;
	return r2;
}

template <class MyFloat>
MyFloat area(MyFloat xmin,MyFloat xmax,MyFloat ymin,MyFloat ymax,int nbx,int nby)
{
	MyFloat dx=(xmax-xmin)/(nbx);
	MyFloat dy=(ymax-ymin)/(nby);
	MyFloat detj=dx*dy;
	MyFloat sum=0.;
	for (int i=0;i<nbx;++i)
	{
		MyFloat x=xmin+dx*i+dx/2;
		for (int j=0;j<nby;++j)
		{
			MyFloat y=ymin+dy*j+dy/2;
			MyFloat res=func<MyFloat>(x,y);
			sum+=res*detj;
		}
	}
	return sum;
}


template <class MyFloat>
MyFloat area_kahan(MyFloat xmin,MyFloat xmax,MyFloat ymin,MyFloat ymax,int nbx,int nby)
{
	MyFloat dx=(xmax-xmin)/(nbx);
	MyFloat dy=(ymax-ymin)/(nby);
	MyFloat detj=dx*dy;
	MyFloat sum=0.;
	MyFloat C=0;
	for (int i=0;i<nbx;++i)
	{
		MyFloat x=xmin+dx*i+dx/2;
		for (int j=0;j<nby;++j)
		{
			MyFloat y=ymin+dy*j+dy/2;
			MyFloat res=func<MyFloat>(x,y);
			MyFloat Y=res*detj-C; 
			MyFloat T=sum+Y; 
			C=(T-sum)-Y;
			sum=T;
		}
	}
	return sum;
}



int main(int argc,char ** argv)
{
	float integ; double integ2;long double integ3; float integ4;
	if ((argc)!= 2) { std::cout << "usage : test n with  n integer >=1" << std::endl; return -1;}
	int n=atoi(argv[1]);
	std::cout << std::setprecision(20);
	integ=area<float>(0.,1.,0.,1,n,n);
	std::cout << "1) sum (float  )=" << integ << std::endl;
	integ2=area<double>(0.,1.,0.,1,n,n);
	std::cout << "2) sum (double )=" << integ2 << std::endl;
	integ3=area<long double>(0.,1.,0.,1,n,n);
	std::cout << "3) sum (ldouble)=" << integ3 << std::endl;
	integ4=area_kahan<float>(0.,1.,0.,1,n,n);
	std::cout << "4) sum (kahan  )=" << integ4 << std::endl;
	return 0;
}