C语言实现龙贝格求积

#include <stdio.h>
#include <stdlib.h>
#include<conio.h>
#include<math.h>
#include<string.h>


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


double Romberg(double top,double bottom,double precision)
{
        int k=1;
        double S,x,T1,T2,S1,S2,C1,C2,R1,R2,h=bottom-top;


/*
转载请注明出处：去转盘网www.quzhuanpan.com
*/


        //S梯形公式，T梯型变步长，S梯形加速，C幸普森加速，R龙贝格求积
        T1=h*(f(top)+f(bottom))/2;//梯形公式
        while(1)
        {
                S=0;
                x=top+h/2;
                do
                {
                   S+=f(x);
                   x+=h;


                }while(x<bottom);
                 T2=(T1+h*S)/2.0;
                 if(fabs(T2-T1)<precision)
                 {
                         return T2;
                 }
                 S2=T2+(T2-T1)/3;//梯形加速
                 if(k==1)
                 {
                         T1=T2;
                         S1=S2;//几下原来的S2
                         h/=2;
                         k+=1;
                         continue;
                 }
                 C2=S2+(S2-S1)/15;//新的S2减原来的S2，即使S1,幸普森加速
                 if(k==2)
                 {
                     C1=C2;
                     T1=T2;
                     S1=S2;
                     h/=2;
                     k+=1;
                     continue;
                 }
                 R2=C2+(C2-C1)/63.0;
                 if(k==3)
                 {
                         R1=R2;
                         C1=C2;
                         T1=T2;
                         S1=S2;
                         h/=2;
                         k+=1;
                         continue;
                 }
                 if(fabs(S2-S1)<precision)
                 {
                         return S2;
                 }
                 C1=C2;
                 T1=T2;
                 S1=S2;
                 h/=2;
                 k+=1;
                 if(fabs(R2-R1)<precision);
                 {
                         return R2;
                 }
        }
}


int main()
{
        double top,bottom,precision,S;
        printf("Pletopse input the begin: ");
        scanf("%lf",&top);
        printf("Pletopse input the end: ");
        scanf("%lf",&bottom);
        printf("Pletopse input the precision:");
        scanf("%lf",&precision);
        S=Romberg(top,bottom,precision);
        printf("The result is:%lf",S);
        getch();
        return 0;
}