#include <functional>
#include <iostream>
#include <math.h>
#include <vector>
using namespace std;

const double pi = acos(-1);

double xm2[2] = {5.857864376269050e-1, 3.414213562373095};
double wm2[2] = {8.535533905932738e-1,1.464466094067262e-1};
double xm4[4] = {3.225476896193923e-1,1.745761101158347,  4.536620296921128,  9.395070912301133};
double wm4[4] = {6.031541043416336e-1,  3.574186924377997e-1,  3.888790851500538e-2,  5.392947055613275e-4 };
double xm6[6] = {2.228466041792607e-1 , 1.188932101672623    , 2.992736326059314     , 5.775143569104511   , 9.837467418382590   , 1.598287398060170e-1};
double wm6[6] = {4.589646739499636e-1  ,  4.170008307721210e-1   , 1.133733820740450e-1  , 1.039919745314907e-2   , 2.610172028149321e-4   , 8.985479064296212e-7};
double xmwm8[8][3] = 
{ 
	0.170279632305, 0.369188589342, 0.437723410493 ,	
	0.903701776799, 0.418786780814, 1.03386934767 ,
	2.25108662987, 0.175794986637 ,1.66970976566 ,
	4.26670017029, 0.0333434922612 ,2.37692470176, 
	7.04590540239, 0.00279453623523 ,3.20854091335 ,
	10.7585160102, 9.07650877338E-005, 4.26857551083 ,
	15.7406786413, 8.48574671626E-007 ,5.81808336867 ,
	22.8631317369, 1.04800117487E-009 ,8.90622621529 
};

class Function : public unary_function<double,double>
{
public:
	Function() {};
	virtual result_type operator()(argument_type x){return x;};
};


class BlackBody : public Function
{
public:
	BlackBody() {};
	double operator()(double x)
	{
		return x*x*x/(1-exp(-x));
	}
};

class Approximation
{
public:
	Approximation(double lowerBound, double upperBound, Function* function)
		:a(lowerBound), b(upperBound), f(function)
	{};

	virtual double calculate(double precision=0, int maxTries = 3000)
	{
		cerr << "Default approximation not implemented" << endl;
		return -1;
	}
protected:
	double a;
	double b;
	Function* f;
	//TBD: set precision & maxTries?
};

class GaussLaguerre : public Approximation
{
public:
	GaussLaguerre(int numberOfPoints, Function* f)	:Approximation(1, numberOfPoints, f), weights(8), abscissas(8)
	{
		int i;
		switch(numberOfPoints)
		{
		case 2:
			for (i=0;i<2;i++)
			{
				weights[i] = wm2[i];
				abscissas[i] = xm2[i];
			}
			break;
		case 4:
			for (i=0;i<4;i++)
			{
				weights[i] = wm4[i];
				abscissas[i] = xm4[i];
			}
			break;
		case 6:
			for (i=0;i<6;i++)
			{
				weights[i] = wm6[i];
				abscissas[i] = xm6[i];
			}
			break;
		case 8:
			for (i=0;i<8;i++)
			{
				weights[i] = xmwm8[i][1];
				abscissas[i] = xmwm8[i][0];
				//cout << abscissas[i] << "\t" << weights[i] << endl;
			}
			break;
		};	
	}

	virtual double calculate()
	{
		double sum=0;
		for (int i =a;i<=b;i++)
			sum += (*f)(abscissas[i-1])*weights[i-1];
		return sum;
	}

protected:
	vector<double> weights;
	vector<double> abscissas;

};

int main(int argc, char* argv)
{
	cout.precision(15);
	cout << "Here are the Gauss Laguerre Approximations of Plank's Black Body equation where pi = " << pi << endl;
	cout << "N\tvalue\terror" << endl;
	BlackBody bb;
	for (int i =2;i<=8;i+=2)
	{
		GaussLaguerre gl(i,&bb);
		cout << i << "\t" << gl.calculate() << "\t" << fabs(gl.calculate() - pi*pi*pi*pi/15) <<endl;; 
	}
return 0;
}

/*5.85786 43762 69050(-1)  3.41421 35623 73095
3.22547 68961 93923(-1)  1.74576 11011 58347  4.53662 02969 21128  9.39507 09123 01133
8.53553 39059 32738(-1)  1.46446 60940 67262(-1)
6.03154 10434 16336(-1)  3.57418 69243 77997(-1)  3.88879 08515 00538(-2)  5.39294 70556 13275(-4)
2.22846 60417 92607(-1)  1.18893 21016 72623     2.99273 63260 59314      5.77514 35691 04511    9.83746 74183 82590    1.59828 73980 60170(1)
4.58964 67394 99636(-1)    4.17000 83077 21210(-1)    1.13373 38207 40450(-1)   1.03991 97453 14907(-2)    2.61017 20281 49321(-4)    8.98547 90642 96212(-7)


1 0.322547689619 0.603154104342 0.832739123838 
2 1.74576110116 0.357418692438 2.04810243845 
3 4.53662029692 0.038887908515 3.63114630582 
4 9.3950709123 0.000539294705561 6.48714508441 


1 0.170279632305 0.369188589342 0.437723410493 
2 0.903701776799 0.418786780814 1.03386934767 
3 2.25108662987 0.175794986637 1.66970976566 
4 4.26670017029 0.0333434922612 2.37692470176 
5 7.04590540239 0.00279453623523 3.20854091335 
6 10.7585160102 9.07650877338E-005 4.26857551083 
7 15.7406786413 8.48574671626E-007 5.81808336867 
8 22.8631317369 1.04800117487E-009 8.90622621529 
*/