算法——龙贝格算法

#include<stdio.h>
#include<math.h>
#define eps 1e-5

double sum=0.0;
double T[100][100];

double f(double x)
{
    double y;
    if(x==0)
      y=1;
    else
      y=sin(x)/x;
    return y;
}
double ss(double x,int n)
{
    int i;
    double s=1.0;
    for(i=1;i<=n;i++)
        s*=x;
    return s;
}
double Romberg(double a,double b)
{
    int i,k,j,n=1;
    double h,s=0;
    h=b-a;
    T[0][0]=h/2*(f(a)+f(b));
    i=0;
    j=0;
    do
    {
         i++;
         j=0;
         h=(b-a)/n;
         s=0.0;
         for(k=0;k<n;k++)
         {
           s+=f(a+(k+0.5)*h);
         }
         T[i][j]=T[i-1][j]/2+h/2*s;
         n=2*n;
         for(j=1;j<=i;j++)
         {
             T[i][j]=(ss(4,j)*T[i][j-1]-T[i-1][j-1])/(ss(4,j)-1);
         }
    }while(fabs(T[i][i]-T[i][i-1])>eps);
     sum=T[i][i];
    return sum;
}

int main()
{
    Romberg(0.0,1.0);
    printf("%.8lf\n",sum);
    return 0;
}