龙贝格算法实现

#include<stdio.h>
#include<math.h>

double f(double x)
{
 double h;
 h=sin(x)/x;
 return(h);
}

void main()
{
 double T1,T2;
 double S,S1,S2;
 double C1,C2;
 double R1,R2;
 double e=1e-6;
 double a,b,h,x;
 int k;
 printf("Input a and b:/n");
 scanf("%lf %lf",&a,&b);
 if(a==0)
  a=0.0000001;
 h=b-a;
 T1=(f(a)+f(b))*h/2;
    for(k=1;k<=3;k++)
 {
  S=0;
  x=a+h/2;
  while(x<b)
  {
    S=S+f(x);
    x=x+h;
  }
       T2=T1/2+S*(h/2);
    S2=T2+1/3*(T2-T1);
    if(k==1){h=h/2;T1=T2;S1=S2; continue;}
    C2=S2+1/15*(S2-S1);
    if(k==2){C1=C2;h=h/2;T1=T2;S1=S2; continue;}
    R2=C2+1/63*(C2-C1);
    if(k==3){R1=R2;C1=C2;h=h/2;T1=T2;S1=S2; continue;}
 }
 do{
  if(fabs(R2-R1)<e)
   break;
  R1=R2;
  C1=C2;
  h=h/2;
  T1=T2;
  S1=S2;
  S=0;
  x=a+h/2;
  while(x<b)
  {
   S=S+f(x);
   x=x+h;
  }
  T2=T1/2+S*(h/2);
  S2=T2+1/3*(T2-T1);
  C2=S2+1/15*(S2-S1);
  R2=C2+1/63*(C2-C1);
 }while(1);
 printf("result:%lf ",R2);
 printf("/n");
}