龍貝格求積公式，是一種梯形法的遞推化，也是將求積區間逐次二分，將前一次二分的積分和Tn與下一次二分的積分和T2n有個遞推關係，通過梯形公式的餘項泰勒展開成級數形式，變數代換，進行整體的加減組合，能夠得到關於T2n和Tn的線性組合S。同理，對S的線性組合可以推出C，對C的線性組合，可以推出R。精度在逐次提高，每一次遞推提高2階。

下面是實現程式碼：

/*龍貝格求積分
*double romb(double a,double b,double eps,double(*f)(double))
*double a:給出積分割槽間下限
*double b:積分割槽間上限
*double eps:誤差精度
*double (*f)(double) :被積函式
*函式返回double
*/
 



#include<iostream>
#include<cmath>

double romb(double a, double b, double eps, double(*f)(double))
{
    int   i, k;
    double y[10], p, x,  q;
    double h = b - a;
    y[0] = h*((*f)(a) + (*f)(b)) / 2.0;
    int m = 1,n=1;
    double ep = eps + 1.0;

    while ((ep >= eps) && (m <= 9
 
))
    {
        p = 0.0;
        for (i = 0; i < n; i++)
        {
            x = a + (i + 0.5)*h;
            p = p + (*f)(x);
        }
        p = (y[0] + h*p) / 2.0;
        for (k = 1; k <= m; k++)
        {
            q = (4.0*p - y[k - 1]) / 3.0;
            y[k - 1] = p;
            p = q;
        }
        ep = fabs
 
(q - y[m - 1]);
        m = m + 1;
        y[m - 1] = q;
        n = n * 2;
        h = h / 2.0;
    }
    return q;
}


//積分函式:f(x)=x/(4+x*x)

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


int main()
{
    double a = 0.0;
    double b = 1.0;
    double eps = 0.0000001;
    double t = romb(a, b, eps, func);
    std::cout << "the value is " << t << std::endl;
    return 0;
}