/*
 *  Simps.c: Simpson Method シンプソン(3/8)則
 */
#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 (3/8) rule.
    Result is stored in "s".
    "eps" is relative error of "s".
 */
{
    double eps = 1e-10;
    int j, jmax, nhlf;
    double h, s1, s2, s3, ss, ds, s;

/*  Initialization */
    jmax = 1; nhlf = 0;
    h = b-a;
    s1 = (*f)(a)+(*f)(b);
    s2 = 0; s3 = 0;
    s = (h/2)*s1;
/*  Start of the Simpson (3/8) method */
    do {
        nhlf++;
        h /= 3;
        s2 += s3; s3 = 0;
        for (j=0; j<=jmax-1; j++)
            s3 += (*f)(a+(3*j+1)*h)+(*f)(a+(3*j+2)*h);
        jmax *= 3;
        ss = s;
        s = (3*h/8)*(s1+2*s2+3*s3);
        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>Simps38
 nhlf =   1 integral =  3.138462e+00 ds =  1.384615e-01
 nhlf =   2 integral =  3.141592e+00 ds =  3.130771e-03
 nhlf =   3 integral =  3.141593e+00 ds =  3.438400e-07
 nhlf =   4 integral =  3.141593e+00 ds =  4.602576e-10
 nhlf =   5 integral =  3.141593e+00 ds =  6.310508e-13
-Integral of 4/(1+x*x) from 0 to 1 =  3.14159265358979
*/