C语言龙贝格求积算法

#include <stdio.h>
#include <math.h>
#define EPS 1e-6
void main(){
double f(double);
double Romberg(double,double,double(*)(double),double);
double r=Romberg(0,1,f,EPS);
printf("积分结果=%.6f\n",r);
}
double f(double x){return 4/(1+x*x);}
double Romberg(double a,double b,double(*f)(double),double eps){
double T[40][40];
int l=0,k,m;
T[0][0]=(b-a)/2*(f(a)+f(b));
printf("T(0,0)=%.6f\n",T[0][0]);
do{
l++;
double fvalue=0;
for(int i=1; i<=(int)pow(2,l-1); i++){
fvalue+=(b-a)/pow(2,l)*f(a+(2*i-1)*(b-a)/pow(2,l));
}
T[0][l]=0.5*T[0][l-1]+fvalue;
printf("T(%d,%d)=%.6f ",0,l,T[0][l]);
for(k=l-1; k>=0; k--){
m=l-k;
T[m][k]=(pow(4,m)*T[m-1][k+1]-T[m-1][k])/(pow(4,m)-1);
printf("T(%d,%d)=%.6f ",m,k,T[m][k]);
}
printf("\n");
}while(fabs(T[m][0]-T[m-1][0])>eps);