/*
//////////////////////////////////////////////////////////////////////////
	This code is share for anybody. If someone want to use this free code. 
Don't erase these copyright descriptions. Otherwise, it is invalid.   
	This code was written by Lai ShengJian, cem@uestc.edu.cn
	If one has some questions about the code, please contact me.
//////////////////////////////////////////////////////////////////////////	
*/


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

#define PI 3.1415926535897932384626433832795

double GaussPoints2[]={-0.57735026918963,
	0.57735026918963};
double GaussPtCeof2[]={1,
	1}; 
double GaussPoints3[]={-0.77459666924148,
	0,
	0.77459666924148};
double GaussPtCeof3[]={0.55555555555556,
	0.88888888888889,
	0.55555555555556}; 
double GaussPoints4[]={-0.86113631159405,
	-0.33998104358486,
	0.33998104358486,
	0.86113631159405};
double GaussPtCeof4[]={0.34785484513745,
	0.65214515486255,
	0.65214515486255,
	0.34785484513745}; 
double GaussPoints5[]={-0.90617984593866,
	-0.53846931010568,
	0,
	0.53846931010568,
	0.90617984593866};
double GaussPtCeof5[]={0.23692688505619,
	0.47862867049937,
	0.56888888888889,
	0.47862867049937,
	0.23692688505619}; 
double GaussPoints6[]={-0.93246951420315,
  -0.66120938646626,
  -0.23861918608320,
   0.23861918608320,
   0.66120938646626,
   0.93246951420315};
double GaussPtCeof6[]={0.17132449237917,
   0.36076157304814,
   0.46791393457269,
   0.46791393457269,
   0.36076157304814,
   0.17132449237917}; 
double GaussPoints7[]={-0.94910791234276,
  -0.74153118559940,
  -0.40584515137740,
                  0,
   0.40584515137740,
   0.74153118559940,
   0.94910791234276};
double GaussPtCeof7[]={0.12948496616887,
   0.27970539148928,
   0.38183005050512,
   0.41795918367347,
   0.38183005050512,
   0.27970539148928,
   0.12948496616887}; 
double GaussPoints8[]={-0.96028985649754,
  -0.79666647741362,
  -0.52553240991633,
  -0.18343464249565,
   0.18343464249565,
   0.52553240991633,
   0.79666647741362,
   0.96028985649754};
double GaussPtCeof8[]={0.10122853629036,
   0.22238103445338,
   0.31370664587789,
   0.36268378337836,
   0.36268378337836,
   0.31370664587789,
   0.22238103445338,
   0.10122853629037}; 
double GaussPoints9[]={-0.96816023950763,
  -0.83603110732663,
  -0.61337143270059,
  -0.32425342340381,
                  0,
   0.32425342340381,
   0.61337143270059,
   0.83603110732664,
   0.96816023950763};
double GaussPtCeof9[]={0.08127438836157,
   0.18064816069486,
   0.26061069640294,
   0.31234707704000,
   0.33023935500126,
   0.31234707704000,
   0.26061069640293,
   0.18064816069486,
   0.08127438836157}; 
double GaussPoints10[]={-0.97390652851717,
  -0.86506336668899,
  -0.67940956829903,
  -0.43339539412925,
  -0.14887433898163,
   0.14887433898163,
   0.43339539412925,
   0.67940956829903,
   0.86506336668898,
   0.97390652851717};
double GaussPtCeof10[]={0.06667134430870,
   0.14945134915058,
   0.21908636251598,
   0.26926671931000,
   0.29552422471475,
   0.29552422471475,
   0.26926671931000,
   0.21908636251598,
   0.14945134915058,
   0.06667134430869}; 
double GaussPoints11[]={-0.97822865814607,
  -0.88706259976809,
  -0.73015200557405,
  -0.51909612920681,
  -0.26954315595234,
                  0,
   0.26954315595234,
   0.51909612920681,
   0.73015200557405,
   0.88706259976808,
   0.97822865814607};
double GaussPtCeof11[]={0.05566856711615,
   0.12558036946493,
   0.18629021092774,
   0.23319376459199,
   0.26280454451025,
   0.27292508677790,
   0.26280454451025,
   0.23319376459199,
   0.18629021092774,
   0.12558036946492,
   0.05566856711614}; 
double GaussPoints12[]={-0.98156063424672,
  -0.90411725637047,
  -0.76990267419431,
  -0.58731795428662,
  -0.36783149899818,
  -0.12523340851147,
   0.12523340851147,
   0.36783149899818,
   0.58731795428661,
   0.76990267419432,
   0.90411725637046,
   0.98156063424673};
double GaussPtCeof12[]={0.04717533638650,
   0.10693932599532,
   0.16007832854335,
   0.20316742672307,
   0.23349253653836,
   0.24914704581340,
   0.24914704581340,
   0.23349253653836,
   0.20316742672307,
   0.16007832854333,
   0.10693932599534,
   0.04717533638649}; 
double GaussPoints13[]={-0.98418305471861,
  -0.91759839922294,
  -0.80157809073333,
  -0.64234933944034,
  -0.44849275103644,
  -0.23045831595513,
                  0,
   0.23045831595514,
   0.44849275103644,
   0.64234933944035,
   0.80157809073329,
   0.91759839922300,
   0.98418305471858};
double GaussPtCeof13[]={0.04048400476525,
   0.09212149983780,
   0.13887351021975,
   0.17814598076195,
   0.20781604753689,
   0.22628318026290,
   0.23255155323087,
   0.22628318026290,
   0.20781604753689,
   0.17814598076194,
   0.13887351021982,
   0.09212149983769,
   0.04048400476534}; 
double GaussPoints14[]={-0.98628380869672,
  -0.92843488366379,
  -0.82720131506959,
  -0.68729290481175,
  -0.51524863635815,
  -0.31911236892789,
  -0.10805494870734,
   0.10805494870734,
   0.31911236892789,
   0.51524863635815,
   0.68729290481175,
   0.82720131506960,
   0.92843488366375,
   0.98628380869674};
double GaussPtCeof14[]={0.03511946033198,
   0.08015808715954,
   0.12151857068806,
   0.15720316715817,
   0.18553839747794,
   0.20519846372130,
   0.21526385346316,
   0.21526385346316,
   0.20519846372130,
   0.18553839747794,
   0.15720316715816,
   0.12151857068800,
   0.08015808715954,
   0.03511946033191}; 
double GaussPoints15[]={-0.98799251802067,
  -0.93727339240030,
  -0.84820658341075,
  -0.72441773136004,
  -0.57097217260858,
  -0.39415134707756,
  -0.20119409399743,
                  0,
   0.20119409399743,
   0.39415134707756,
   0.57097217260857,
   0.72441773136009,
   0.84820658341064,
   0.93727339240040,
   0.98799251802063};
double GaussPtCeof15[]={0.03075324199563,
   0.07036604748857,
   0.10715922046708,
   0.13957067792619,
   0.16626920581700,
   0.18616100001556,
   0.19843148532711,
   0.20257824192556,
   0.19843148532711,
   0.18616100001557,
   0.16626920581698,
   0.13957067792617,
   0.10715922046709,
   0.07036604748844,
   0.03075324199577}; 


double IntegralGaussLegendre(double(*fun)(double),double a,double b,double eps=1e-7,int nGauss=10) 
{ 
	double preVal,val,valtmp,normxt;
	double *xt,*ceo;
	int i,j,m;
	double a0,b0,dx;
	const int mIter=10;

	switch(nGauss) {
	case 2:
		xt=&GaussPoints2[0];
		ceo=&GaussPtCeof2[0];
		break;
	case 3:
		xt=&GaussPoints3[0];
		ceo=&GaussPtCeof3[0];
		break;
	case 4:
		xt=&GaussPoints4[0];
		ceo=&GaussPtCeof4[0];
		break;
	case 5:
		xt=&GaussPoints5[0];
		ceo=&GaussPtCeof5[0];
		break;
	case 6:
		xt=&GaussPoints6[0];
		ceo=&GaussPtCeof6[0];
		break;
	case 7:
		xt=&GaussPoints7[0];
		ceo=&GaussPtCeof7[0];
		break;
	case 8:
		xt=&GaussPoints8[0];
		ceo=&GaussPtCeof8[0];
		break;
	case 9:
		xt=&GaussPoints9[0];
		ceo=&GaussPtCeof9[0];
		break;
	case 10:
		xt=&GaussPoints10[0];
		ceo=&GaussPtCeof10[0];
		break;
	case 11:
		xt=&GaussPoints11[0];
		ceo=&GaussPtCeof11[0];
		break;
	case 12:
		xt=&GaussPoints12[0];
		ceo=&GaussPtCeof12[0];
		break;
	case 13:
		xt=&GaussPoints13[0];
		ceo=&GaussPtCeof13[0];
		break;
	case 14:
		xt=&GaussPoints14[0];
		ceo=&GaussPtCeof14[0];
		break;
	case 15:
		xt=&GaussPoints15[0];
		ceo=&GaussPtCeof15[0];
		break;
	default:
		std::cout<<" number of Gauss points must be between 2 and 15!! "<<std::endl;
		return 0;
		break;
	}
	preVal=0;
	m=1;
	while (true) {
		val=0.0; 
		dx=(b-a)/m;
		for(j=0;j<m;j++){
			a0=a+j*dx;
			b0=a+(j+1)*dx;
			valtmp=0;
			for(i=0;i<nGauss;i++) 
			{ 
				normxt=(xt[i]*(b0-a0)+a0+b0)/2.0; 
				valtmp+=ceo[i]*fun(normxt); 
			}
			val+=valtmp;
		}
		val=val*dx/2;
		if (fabs(val-preVal)/fabs(val)<eps) break;
		preVal=val;
		++m;
		if(m>mIter) {
			std::cout<<" !!! number of iter have already more than"<<mIter<<std::endl;
			break;
		}
	}
	return val;

} 

double IntegralGaussLegendre2D_1(double(*fun)(double,double),double x0,double a,double b,double eps=1e-7,int nGauss=10)
{
	double preVal,val,valtmp,normxt;
	double *xt,*ceo;
	int i,j,m;
	double a0,b0,dx;
	const int mIter=10;

	switch(nGauss) {
	case 2:
		xt=&GaussPoints2[0];
		ceo=&GaussPtCeof2[0];
		break;
	case 3:
		xt=&GaussPoints3[0];
		ceo=&GaussPtCeof3[0];
		break;
	case 4:
		xt=&GaussPoints4[0];
		ceo=&GaussPtCeof4[0];
		break;
	case 5:
		xt=&GaussPoints5[0];
		ceo=&GaussPtCeof5[0];
		break;
	case 6:
		xt=&GaussPoints6[0];
		ceo=&GaussPtCeof6[0];
		break;
	case 7:
		xt=&GaussPoints7[0];
		ceo=&GaussPtCeof7[0];
		break;
	case 8:
		xt=&GaussPoints8[0];
		ceo=&GaussPtCeof8[0];
		break;
	case 9:
		xt=&GaussPoints9[0];
		ceo=&GaussPtCeof9[0];
		break;
	case 10:
		xt=&GaussPoints10[0];
		ceo=&GaussPtCeof10[0];
		break;
	case 11:
		xt=&GaussPoints11[0];
		ceo=&GaussPtCeof11[0];
		break;
	case 12: