/*
 *  Rombrg.c: Romberg Integral Method ロンバーグ積分
 */
#include <stdio.h>
#include <math.h>
#define MXHLF 30

double f(double x)
/*  f(x) = 4/(1+x^2) */
{
    return 4.0/(1+x*x);
}

double Romberg(double a, double b, double (*f)(double x), int mxhlf)
/*
    Integration of function "f(x)" from "a" to "b"
    by using Romberg integral method.
    Result is stored in "s".
    "eps" is relative error of "s".
 */
{
    double eps = 1e-10;
    int j, jmax, nhlf;
    double h, sum, ss, s0, r, s, ds;
    double T[MXHLF+1];

/*  Initialization */
    jmax = 1; nhlf = 0;
    h = b-a;
    s = (h/2.0)*((*f)(a)+(*f)(b));
    T[0] = s;
/*  Start of the Romberg integral method */
    do {
    /*  台形則 */
        nhlf++;
        h /= 2.0; sum = 0;
        for (j=1; j<=jmax; j++) sum += (*f)(a+(2*j-1)*h);
        jmax *= 2;
        ss = s;
        s = 0.5*T[0]+h*sum;
    /*  リチャードソンの補外法 */
        r = 1;
        for (j=1; j<=nhlf; j++) {
            s0 = T[j-1];
            T[j-1]= s;
            r = 4*r;
            s = (r*s-s0)/(r-1);
        }
        T[nhlf] = s;
    /*  収束判定 */
        ds = fabs(s-ss);
        printf(" nhlf = %3d integral = %15.12f ds =  %e\n",nhlf,s,ds);
        if (nhlf == mxhlf)
            printf("*** Simpson method does not converge. MXHLF = %d", mxhlf);
    } while ((ds >= eps*fabs(s)) && (nhlf < mxhlf));
    return s;
}

/*  Main  */
int main()
{
   double a, b, s;

    a = 0;
    b = 1;
    s = Romberg(a, b, f, MXHLF);
    printf("-Integral of 4/(1+x*x) from 0 to 1 = %15.12f\n", s);
    return 0;
}
/*  ---------   実行例: 教科書 p.114 例題5.2   ---------  */
/*
Z:\src\C>Rombrg
 nhlf =   1 integral =  3.133333333333 ds =  1.333333e-01
 nhlf =   2 integral =  3.142117647059 ds =  8.784314e-03
 nhlf =   3 integral =  3.141585783762 ds =  5.318633e-04
 nhlf =   4 integral =  3.141592665278 ds =  6.881516e-06
 nhlf =   5 integral =  3.141592653638 ds =  1.163947e-08
 nhlf =   6 integral =  3.141592653590 ds =  4.852119e-11
-Integral of 4/(1+x*x) from 0 to 1 =  3.141592653590
*/