用c语言数值求积的Romberg方法计算圆周率π

罗马堡算法是一种高精度数值积分方法，常用于计算难以精确解析的函数的定积分。在C语言中，我们可以利用该算法来逼近圆周率π的值，因为π可以表示为某个特定函数在其区间上的面积。这里是一个简单的例子，我们将使用矩形规则逐步改进到 Romberg 近似：

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

// 计算n阶Romberg公式的近似值
double romberg(int n, double a, double b, double (*f)(double)) {
    double r = (b - a) / 4.0, h = r;
    double T[5][5]; // 用于存储中间结果

    // 初始化第一行和第一列
    for (int i = 0; i <= n; ++i)
        T[i][0] = f(a + (i - 1) * h);
        
    // 使用Romberg公式计算其余的T[i][j]
    for (int j = 1; j <= n; ++j) {
        for (int i = 1; i <= j; ++i) {
            if (i % 2 == 1)
                T[i][j] = (4.0 * pow(h, 2j)) * (T[i][j - 1] - T[i - 1][j - 1]);
            else
                T[i][j] = (h * (T[i + 1][j - 1] + 4.0 * T[i][j - 1]) + T[i - 1][j - 1]) / 3.0;
        }
    }

    return T[n][n];
}

// 示例函数，仅作积分示例，不是真正计算π
double func(double x) {
    return 1.0 / sqrt(1.0 - x * x); // 用以计算sin(x)^(-1)
}

int main() {
    const int N = 100; // 精度等级
    double pi_approx = romberg(N, 0, 1, func); // 积分区间从0到1，积分函数为反正切
    printf("Using Romberg method with %d levels of refinement, pi approximation is %.9f\n", N, pi_approx);

    return 0;
}