c++实现龙贝格积分法

`
double 龙贝格(double a, double b, double esp)//输入范围a~b,误差限度esp
{
double T[2] = { 0.0,0.0 };//复化梯形
double S[2] = { 0.0,0.0 };//复化辛浦生
double C[2] = { 0.0,0.0 };//复化牛顿-柯特斯
double R[2] = { 0.0,0.0 };//龙贝格

int k = 1;
double h = b - a;
T[0] = (h / 2.0) * (func(a) + func(b));//自定义函数func()
cout << "k=" << k << " ";
cout <<"T[0]："<< T[0]<<" ";
if (CT) {
    cout << endl;
    CT = false;
}
while (true) {
    cout << "k=" << k << " ";
    double Sum = 0;
    double x = a + h / 2.0;
    while (true) {
        Sum += func(x);
        x += h;
        if (x >= b)
            break;
    }
    T[1] = T[0] / 2.0 + (h * Sum) / 2.0;
    cout << "T[1]：" << T[1] << " ";
    S[1] = T[1] + (1.0 / 3.0) * (T[1] - T[0]);
    cout << "S[1]：" << S[1]<<" ";
    if (CS) {
        cout << endl;
        CS = false;
    }
    if (k == 1) {
        k += 1;
        h = h / 2.0;
        T[0] = T[1];
        S[0] = S[1];
        continue;
    }
    else {
        C[1] = S[1] + (1.0 / 15.0) * (S[1] - S[0]);
        cout << "C[1]：" << C[1]<<" ";
        if (CC) {
            cout << endl;
            CC = false;
        }
        if (k == 2) {
            C[0] = C[1];
            k += 1;
            h = h / 2.0;
            T[0] = T[1];
            S[0] = S[1];
            continue;
        }
        else {
            R[1] = C[1] + (1.0 / 63.0) * (C[1] - C[0]);
            cout << "R[1]：" << R[1] << endl;
            if (k == 3) {
                R[0] = R[1];
                C[0] = C[1];
                k += 1;
                h = h / 2.0;
                T[0] = T[1];
                S[0] = S[1];
                continue;
            }
            else {
                if (std::abs(R[1] - R[0]) < esp)
                    return R[1];
                else {
                    R[0] = R[1];
                    C[0] = C[1];
                    k += 1;
                    h = h / 2.0;
                    T[0] = T[1];
                    S[0] = S[1];
                    continue;
                }
            }
        }
    }
}

}
`
例如如下函数：

double func(double x) { if (x == 0) return 1; return sin(x) / x; }

求0~1的函数积分，结果可得