/*
 *  Trpex.c: Trapezoidal 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 Trapez(double a, double b, double (*f)(double x), int mxhlf)
/*
    Integration of function "f(x)" from "a" to "b"
    by using Trapezoidal formula.
    Result is stored in "s".
    "eps" is relative error of "s".
 */
{
    double eps = 1e-10;
    int j, jmax, nhlf;
    double h, sum, ss, ds, s;

/*  Initialization */
    jmax = 1; nhlf = 0;
    h = b-a;
    s = (h/2)*((*f)(a)+(*f)(b));
/*  Start of the trapezoidal 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*s+h*sum;
        ds = fabs(s-ss);
        printf(" nhlf = %3d integral = %e ds = %e\n",nhlf,s,ds);
        if (nhlf == mxhlf)
            printf("*** Simpson method does not converge. MXHLF = %d\n", mxhlf);
    } while ((ds >= eps*fabs(s)) && (nhlf < mxhlf));
    return s;
}
/*  Main  */
int main()
{
    double a, b, s;

    a = 0;
    b = 1;
    s = Trapez(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>Trpex
 nhlf =   1 integral =  3.100000e+00 ds =  1.000000e-01
 nhlf =   2 integral =  3.131176e+00 ds =  3.117647e-02
 nhlf =   3 integral =  3.138988e+00 ds =  7.812024e-03
 nhlf =   4 integral =  3.140942e+00 ds =  1.953118e-03
 nhlf =   5 integral =  3.141430e+00 ds =  4.882811e-04
 nhlf =   6 integral =  3.141552e+00 ds =  1.220703e-04
 nhlf =   7 integral =  3.141582e+00 ds =  3.051758e-05
 nhlf =   8 integral =  3.141590e+00 ds =  7.629395e-06
 nhlf =   9 integral =  3.141592e+00 ds =  1.907349e-06
 nhlf =  10 integral =  3.141592e+00 ds =  4.768372e-07
 nhlf =  11 integral =  3.141593e+00 ds =  1.192093e-07
 nhlf =  12 integral =  3.141593e+00 ds =  2.980232e-08
 nhlf =  13 integral =  3.141593e+00 ds =  7.450580e-09
 nhlf =  14 integral =  3.141593e+00 ds =  1.862646e-09
 nhlf =  15 integral =  3.141593e+00 ds =  4.656560e-10
 nhlf =  16 integral =  3.141593e+00 ds =  1.164198e-10
-Integral of 4/(1+x*x) from 0 to 1 =  3.141592653551
*/