/*
 *  Gausin.c: Gauss Integral Method ガウスの積分公式
 */
#include <stdio.h>
#include <math.h>

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

double Gauss(double a, double b, double f(double x), double *error, int *nmax)
/*
    Integration of function "f(x)" from "a" to "b"
    by using Gauss's method.
    Result is stored in "S".
    "eps" is relative error of "S".
 */
{
    double eps = 1e-10;
    double axis[16] = {0, 
          0.57735026918962576451, 0.77459666924148337703, 0.0,
          0.86113631159405257522, 0.33998104358485626480,
          0.90617984593866399280, 0.53846931010586309104, 0.0,
          0.93246951420315202781, 0.66120938646626451336,
          0.23861918608319690863, 0.94910791234275852453,
          0.74153118559939443986, 0.40584515137739716691, 0.0 };
    double weight[16] = {0, 1.0,
          0.55555555555555555556, 0.88888888888888888889,
          0.34785484513745385737, 0.65214515486254614263,
          0.23692688505618908751, 0.47862867049936646804,
          0.56888888888888888889, 0.17132449237917034504,
          0.36076157304813860757, 0.46791393457269104739,
          0.12948496616886969327, 0.27970539148927666790,
          0.38183005050511894495, 0.41795918367346938776 };
    int start[8] = {0,1,2,4,6,9,12,16};
    int i, j, k, n;
    double c1, c2, S, S0, sum, xi;
    double x[8], w[8];

//  Initialization
    c1 = (b-a)/2; c2 = (b+a)/2;
    S = 2*c1*(*f)(c2);
//  Start of the Gaussian integral method
    for (n=2; n<=7; n++)
    {
        S0 = S;
    //  n次の積分値 dS
    //  座標と重みの設定
        k = 0;
        for (i=start[n-1]; i<=start[n]-1; i++)
        {
            k++; j = n-k+1;
            x[k] = -axis[i]; w[k] = weight[i];
            x[j] =  axis[i]; w[j] = weight[i];
        }
    //  積分値の計算
        sum = 0;
        for (i=1; i<=n; i++)
        {
            xi = c1*x[i]+c2;
            sum += w[i]*(*f)(xi);
        }
        S = c1*sum;
    //  収束判定
        *error = fabs(S-S0);
        if (*error < eps*fabs(S)) break;
        if (n == 7)
            printf("*** Gaussian method does not converge. nmax = 7\n");
    }
    *nmax = n;
    return S;
}
/*   -----   Main   -----   */
int main()
{
    int i, nmax, N;
    double a, b, d, S, dS, xmin, xmax, error;

    N = 8;
    a = 0; b = 1;
    d = (b-a)/N; S = 0;
    for (i=1; i<=N; i++)
    {
        xmin = a+(i-1)*d; xmax = xmin+d;
        dS = Gauss(xmin, xmax, f, &error, &nmax);
        S += dS;
    }
    printf("-Integral of 4/(1+x*x) from 0 to 1 = %17.14f\n",S);
    printf(" error =  %e     n =  %3d\n", error, nmax);
    return 0;
}
/*   ---------   実行例: 教科書 p.130 例題5.6   ---------   */
/*
Z:\src\C>Gausin
-Integral of 4/(1+x*x) from 0 to 1 =  3.14159265358972
 error =  3.491651e-14     n =  5
*/