/*  Romberg求定積分
    輸入：積分區間[a,b]，被積函數f(x,y,z)
    輸出：積分結果

    f(x,y,z)示例：
    double f0( double x, double l, double t )
    {
	    return sqrt(1.0+l*l*t*t*cos(t*x)*cos(t*x));
    }
*/

double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, 
				double l, double t)

double Romberg (double a, double b, double (*f)(double x, double y, double z), double eps, 
				double l, double t)
{
#define MAX_N  1000

	int i, j, temp2, min;
	double h, R[2][MAX_N], temp4;

	for (i=0; i<MAX_N; i++) {
		R[0][i] = 0.0;
		R[1][i] = 0.0;
	}
	h = b-a;
	min = (int)(log(h*10.0)/log(2.0)); //h should be at most 0.1
	R[0][0] = ((*f)(a, l, t)+(*f)(b, l, t))*h*0.50;
	i = 1;
	temp2 = 1;
	while (i<MAX_N){
		i++;
		R[1][0] = 0.0;
		for (j=1; j<=temp2; j++)
			R[1][0] += (*f)(a+h*((double)j-0.50), l, t);
		R[1][0] = (R[0][0] + h*R[1][0])*0.50;
		temp4 = 4.0;
		for (j=1; j<i; j++) {
			R[1][j] = R[1][j-1] + (R[1][j-1]-R[0][j-1])/(temp4-1.0);
			temp4 *= 4.0;
		}
		if ((fabs(R[1][i-1]-R[0][i-2])<eps)&&(i>min))
			return R[1][i-1];
		h *= 0.50;
		temp2 *= 2;
		for (j=0; j<i; j++)
			R[0][j] = R[1][j];
	}
	return R[1][MAX_N-1];
}

double Integral(double a, double b, double (*f)(double x, double y, double z), double eps, 
				double l, double t)
{
#define pi 3.1415926535897932

	int n;
	double R, p, res;

	n = (int)(floor)(b * t * 0.50 / pi);
	p = 2.0 * pi / t;
	res = b - (double)n * p;
	if (n) 
		R = Romberg (a, p, f0, eps/(double)n, l, t);
	R = R * (double)n + Romberg( 0.0, res, f0, eps, l, t );

	return R/100.0;
}