/*
 *  Simps.c: Simpson Method  シンプソン(1/3)則
 */
#include <stdio.h>
#include <math.h>
#define MXHLF  50

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

double Simpsn(double a, double b, double (*f)(double x), int mxhlf)
/*
    Integration of function "f(x)" from "a" to "b"
    by using Simpson (1/3) rule.
    Result is stored in "s".
    "eps" is relative error of "s".
 */
{
    double eps = 1e-10;
    int j, jmax, nhlf;
    double h, s1, s2, s4, ss, ds, s;

/*  Initialization */
    jmax = 1; nhlf = 0;
    h = b-a;
    s1 = (*f)(a)+(*f)(b);
    s2 = 0; s4 = 0;
    s = (h/2)*s1;
/*  Start of the Simpson method */
    do {
        nhlf++;
        h /= 2.0;
        s2 += s4;
        s4 = 0;
        for (j=1; j<=jmax; j++) s4 += (*f)(a+(2*j-1)*h);
        jmax *= 2;
        ss = s;
        s = (h/3)*(s1+2*s2+4*s4);
        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 = \n", mxhlf);
    } while ((ds >= eps*fabs(s)) && (nhlf < mxhlf));
    return s;
}
/*  Main  */
int main()
{
    double a, b, s;
    a = 0;
    b = 1;
    s = Simpsn(a, b, f, MXHLF);
    printf("-Integral of 4/(1+x*x) from 0 to 1 = %17.14f\n", s);
    return 0;
}
/*  ---------   実行例: 教科書 p.119 例題5.4   ---------  */
/*
Z:\src\C>Simps
 nhlf =   1 integral =  3.133333e+00 ds =  1.333333e-01
 nhlf =   2 integral =  3.141569e+00 ds =  8.235294e-03
 nhlf =   3 integral =  3.141593e+00 ds =  2.387501e-05
 nhlf =   4 integral =  3.141593e+00 ds =  1.487661e-07
 nhlf =   5 integral =  3.141593e+00 ds =  2.328014e-09
 nhlf =   6 integral =  3.141593e+00 ds =  3.637979e-11
-Integral of 4/(1+x*x) from 0 to 1 =  3.14159265358922
*/