#include "stdio.h"
#include "conio.h"
#include "math.h"
#define EPS 1e-6

double f(double x) { return 4/(1+x*x); }

double Romberg(double a,double b,double(*f)(double),double eps)
{
    int i , m , n=1 , j=1 ;
    double h , sum , p ;
    double T[10][10];

    h = (b-a)/2;
    T[0][0] = h*( f(a) + f(b) );
    printf("T(0,0)=%f\n",T[0][0]);

    do                                //Romberg算法
    {

         sum=0.0;
        for( i=1 ; i <= n ; i++ )
          sum = sum + f( a + ((2*i-1)*h) );

        T[j][0] = ( T[j-1][0]/2) + ( h*sum );
        printf("T(0,%d)=%f\n",j,T[j][0]);

        for( m=1 ; m <= j ; m++ )
        {
            p = pow( 4 , (double)(m) );

            T[j-m][m] = ( p * T[j-m+1][m-1] - T[j-m][m-1] ) / (p-1);
            printf("T(%d,%d)=%f",m,j-m,T[j-m][m]);
        }

        m--;
        h = h/2;
        n = n*2;
        j++;
       
    } while ( fabs( T[0][m] - T[0][m-1] ) > eps );
    return (T[0][m]);
}

main()
{
    printf("\n\n积分结果 =  %f\n", Romberg(0 , 1 , f , EPS));
}