#include<iostream.h>
#include<math.h>
#define pa 3.1415926
double f1(double x);
double f2(double x);
double f3(double x);
double cT(double a,double b,int m,int k);//組合梯形求積公式,積分限a,b，對第k個函數積分
double cSimpson(double a,double b,int m,int k)//組合simpson求積公式
{
	double h=(b-a)/(2*m);
	double F1=0;//存放f2j-1
	double F2=0;//存放f2j
	double F0;
	int j;
	switch(k)
	{
	case 1:
	F0=f1(a)+f1(b);
	for(j=1;j<2*m;j++)
	{
		double Xj=a+j*h;
		if(j%2==0)
			F2=F2+f1(Xj);
		else
			F1=F1+f1(Xj);
	}
	break;
	case 2:
    F0=f2(a)+f2(b);
	for(j=1;j<2*m;j++)
	{
		double Xj=a+j*h;
		if(j%2==0)
			F2=F2+f2(Xj);
		else
			F1=F1+f2(Xj);
	}
	break;
    case 3:
    F0=f3(a)+f3(b);
	for(j=1;j<2*m;j++)
	{
		double Xj=a+j*h;
		if(j%2==0)
			F2=F2+f3(Xj);
		else
			F1=F1+f3(Xj);
	}
	break;
	default:
		cout<<"error!";
	}
	
	double Xf=h*(F0+2*F2+4*F1)/3;
	return Xf;

}

void main()
{
/*	double sf=cSimpson(0,1,10);
	double tf=cT(0,1,100);
	double re=pa*sqrt(2)/4;
	cout<<sf;
	cout<<endl;
	cout<<fabs(re-sf);
	cout<<endl;
	cout<<tf;
	cout<<endl;
	cout<<fabs(tf-re);*/
	double re1=pa*sqrt(2)/4;
	double re2=0.5;
	double re3=pa/4;
	///f1//
	double xfs=cSimpson(0,1,10,1);
	double xft=cT(0,1,300,1);
    cout<<"(1): "<<endl;
	cout<<"simpson積分值：  "<<xfs;
	cout<<"    誤差:  "<<fabs(xfs-re1)<<endl;
	cout<<"梯形積分值：  "<<xft;
	cout<<"    誤差:  "<<fabs(xft-re1)<<endl;
	//f2//
	xfs=cSimpson(0,1,1000,2);
	xft=cT(0,1,4000,2);
    cout<<"(2): "<<endl;
	cout<<"simpson積分值：  "<<xfs;
	cout<<"    誤差:  "<<fabs(xfs-re2)<<endl;
	cout<<"梯形積分值：  "<<xft;
	cout<<"    誤差:  "<<fabs(xft-re2)<<endl;
	//f3//
	xfs=cSimpson(0,pa/2,10,3);
	xft=cT(0,pa/2,10,3);
    cout<<"(3): "<<endl;
	cout<<"simpson積分值：  "<<xfs;
	cout<<"    誤差:  "<<fabs(xfs-re3)<<endl;
	cout<<"梯形積分值：  "<<xft;
	cout<<"    誤差:  "<<fabs(xft-re3)<<endl;


	
}
double cT(double a,double b,int m,int k)//組合梯形求積公式,積分限a,b，對第k個函數積分
{
    double h=(b-a)/m;
	double Fj=0;//
	double F0;
	int j;
	switch(k)
	{
	case 1:
	F0=f1(a)+f1(b);
	for(j=1;j<m;j++)
	{
		double Xj=a+j*h;
		Fj+=f1(Xj);
	}
	break;
	case 2:
	F0=f2(a)+f2(b);
	for(j=1;j<m;j++)
	{
		double Xj=a+j*h;
		Fj+=f2(Xj);
	}
	break;
	case 3:
	F0=f3(a)+f3(b);
	for(j=1;j<m;j++)
	{
		double Xj=a+j*h;
		Fj+=f3(Xj);
	}
	break;
    default:
		cout<<"error!";
	}
	
	double Xf=h*(2*Fj+F0)/2;
	return Xf;
}
double f1(double x)
{
	return (1+x*x)/(1+x*x*x*x);
}
double f2(double x)
{
	return log(1+sqrt(x));
}
double f3(double x)
{
	return 1/(pow(cos(x),2)+4*pow(sin(x),2));
}