Romberg积分算法的C语言实现

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

#define f(x) log(x)
#define N_H 1
#define MAXREPT 10
#define a 1.0
#define b 2.0
#define epsilon pow(10, -4)

double calculateT(double aa, double bb, long int n)
{
    double h = (bb-aa)/n;
    double sum = 0;
    for (int i = 1; i < n; i++)
    {
        sum += f(aa+i*h);
    }
    sum += (f(aa) + f(bb))/2;
    return (h*sum);
}

void printR(int m, double R[MAXREPT][MAXREPT])
{
    for (int i = 0; i < MAXREPT; i++) //calculate R[0]
    {
        printf("%lf\n", R[m][i]);
    }
}

int main()
{
    int i;
    long int n = N_H;
    double R[MAXREPT][MAXREPT] = {0}; //R(k, j) is R[j-1][k-1]
    
    for (int m = 0; m < MAXREPT; m++)
    {       
        if (m == 0)//calculate R[0]
        {
            for (int i = 0; i < MAXREPT; i++) 
            {
                R[0][i] = calculateT(a, b, n*pow(2, i));
            }
            printf("j = %d, R = \n", m+1);
            printR(0, R);
            continue;
        }
        
        for (int i = m; i < MAXREPT; i++)
        {
            R[m][i] = R[m-1][i] + (R[m-1][i] - R[m-1][i-1])/(pow(4, m) - 1);
        }
        printf("j = %d, R = \n", m+1);
        printR(m, R);

        double limit = fabs(R[m][m] - R[m-1][m-1]);
        if (limit < epsilon)
        {
            printf("The Integrate is %lf \n", R[m][m]);
            break;
        }
        
    }
    
    
    return 0;
}