龙贝格求积法

#include<iostream>
using namespace std;
#include"math.h"
#define e  1e-5
double  f(double x)
{
    double y;
    if (x==0)
        return y=1.0;//f(0)=1
    else y=sin(x)/x;
    return y;
}
void romberg(double a,double b)
{
    int n=1,k=0,i;
    double h,T2,S2=0,C2=0,R2=0,T1,C1,S1,R1;
    double x;
    h=b-a;
    T2=h/2*(f(a)+f(b));
    do//在要求精度之下
    {
        x=a+h/2;
        R1=R2;
        T1=T2;
        S1=S2;
        C1=C2;
        double sum=0;
        for(i=0;i<=pow(2,k)-1;i++)
           {
               sum=sum+f(x);///都是公式
               x+=h;
           }
        T2=(T1+sum*h)/2;
        S2=(4*T2-T1)/3;
        C2=(16*S2-S1)/15;
        R2=(64*C2-C1)/63;
        k++;
        h=h/2;
    }while (fabs((R2-R1))>e);
    cout<<R2<<endl;
}
int main()
{
    double a,b;
    cout<<"输入上下限a,b的值用空格隔开"<<endl;
    cin>>a>>b;
    cout<<"下限a="<<a<<endl;
    cout<<"上限b="<<b<<endl;
    cout<<"被积函数为：y=sin(x)/x"<<endl;
    romberg(a,b);
}